Actual source code: ex27.c
2: static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
3: Test MatMatSolve(). Input parameters include\n\
4: -f <input_file> : file to load \n\n";
6: /*
7: Usage:
8: ex27 -f0 <mat_binaryfile>
9: */
11: #include petscksp.h
12: EXTERN PetscErrorCode PCShellApply_Matinv(PC,Vec,Vec);
16: int main(int argc,char **args)
17: {
18: KSP ksp;
19: Mat A,B,F,X;
20: Vec x,b,u; /* approx solution, RHS, exact solution */
21: PetscViewer fd; /* viewer */
22: char file[1][PETSC_MAX_PATH_LEN]; /* input file name */
23: PetscTruth flg;
25: PetscInt M,N,i,its;
26: PetscReal norm;
27: PetscScalar val=1.0;
28: PetscMPIInt size;
29: PC pc;
31: PetscInitialize(&argc,&args,(char *)0,help);
32: MPI_Comm_size(PETSC_COMM_WORLD,&size);
33: if (size != 1) SETERRQ(PETSC_ERR_SUP,"This is a uniprocessor example only!");
35: /* Read matrix and right-hand-side vector */
36: PetscOptionsGetString(PETSC_NULL,"-f",file[0],PETSC_MAX_PATH_LEN-1,&flg);
37: if (!flg) {
38: SETERRQ(PETSC_ERR_USER,"Must indicate binary file with the -f option");
39: }
41: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&fd);
42: MatLoad(fd,MATAIJ,&A);
43: PetscExceptionTry1(VecLoad(fd,PETSC_NULL,&b),PETSC_ERR_FILE_READ);
44: if (PetscExceptionCaught(ierr,PETSC_ERR_FILE_UNEXPECTED) || PetscExceptionCaught(ierr,PETSC_ERR_FILE_READ)) {
45: /* if file contains no RHS, then use a vector of all ones */
46: PetscInt m;
47: PetscScalar one = 1.0;
48: PetscInfo(0,"Using vector of ones for RHS\n");
49: MatGetLocalSize(A,&m,PETSC_NULL);
50: VecCreate(PETSC_COMM_WORLD,&b);
51: VecSetSizes(b,m,PETSC_DECIDE);
52: VecSetFromOptions(b);
53: VecSet(b,one);
54: } else
55: PetscViewerDestroy(fd);
57: /*
58: If the loaded matrix is larger than the vector (due to being padded
59: to match the block size of the system), then create a new padded vector.
60: */
61: {
62: PetscInt m,n,j,mvec,start,end,indx;
63: Vec tmp;
64: PetscScalar *bold;
66: /* Create a new vector b by padding the old one */
67: MatGetLocalSize(A,&m,&n);
68: if (m != n) {
69: SETERRQ2(PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
70: }
71: VecCreate(PETSC_COMM_WORLD,&tmp);
72: VecSetSizes(tmp,m,PETSC_DECIDE);
73: VecSetFromOptions(tmp);
74: VecGetOwnershipRange(b,&start,&end);
75: VecGetLocalSize(b,&mvec);
76: VecGetArray(b,&bold);
77: for (j=0; j<mvec; j++) {
78: indx = start+j;
79: VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);
80: }
81: VecRestoreArray(b,&bold);
82: VecDestroy(b);
83: VecAssemblyBegin(tmp);
84: VecAssemblyEnd(tmp);
85: b = tmp;
86: }
87: VecDuplicate(b,&x);
88: VecDuplicate(b,&u);
89: VecSet(x,0.0);
91: /* Create dense matric B and X. Set B as an identity matrix */
92: MatGetSize(A,&M,&N);
93: MatCreate(MPI_COMM_SELF,&B);
94: MatSetSizes(B,M,N,M,N);
95: MatSetType(B,MATSEQDENSE);
96: MatSeqDenseSetPreallocation(B,PETSC_NULL);
97: for (i=0; i<M; i++){
98: MatSetValues(B,1,&i,1,&i,&val,INSERT_VALUES);
99: }
100: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
101: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
103: MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&X);
104:
105: /* Compute X=inv(A) by MatMatSolve() */
106: KSPCreate(PETSC_COMM_WORLD,&ksp);
107: KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
108: KSPGetPC(ksp,&pc);
109: PCSetType(pc,PCLU);
110: KSPSetUp(ksp);
111: PCFactorGetMatrix(pc,&F);
112: MatMatSolve(F,B,X);
113: MatDestroy(B);
115: /* Now, set X=inv(A) as a preconditioner */
116: PCSetType(pc,PCSHELL);
117: PCShellSetContext(pc,(void *)X);
118: PCShellSetApply(pc,PCShellApply_Matinv);
119: KSPSetFromOptions(ksp);
121: /* Sove preconditioned system A*x = b */
122: KSPSolve(ksp,b,x);
123: KSPGetIterationNumber(ksp,&its);
125: /* Check error */
126: MatMult(A,x,u);
127: VecAXPY(u,-1.0,b);
128: VecNorm(u,NORM_2,&norm);
129: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
130: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %A\n",norm);
131:
132: /* Free work space. */
133: MatDestroy(X);
134: MatDestroy(A); VecDestroy(b);
135: VecDestroy(u); VecDestroy(x);
136: KSPDestroy(ksp);
137: PetscFinalize();
138: return 0;
139: }
141: PetscErrorCode PCShellApply_Matinv(PC pc,Vec xin,Vec xout)
142: {
144: Mat X;
147: PCShellGetContext(pc,(void**)&X);
148: MatMult(X,xin,xout);
149: return(0);
150: }