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);
14: int main(int argc,char **args)
15: {
16: KSP ksp;
17: Mat A,B,F,X;
18: Vec x,b,u; /* approx solution, RHS, exact solution */
19: PetscViewer fd; /* viewer */
20: char file[1][PETSC_MAX_PATH_LEN]; /* input file name */
21: PetscBool flg;
22: PetscInt M,N,i,its;
23: PetscReal norm;
24: PetscScalar val=1.0;
25: PetscMPIInt size;
26: PC pc;
28: PetscInitialize(&argc,&args,(char*)0,help);
29: MPI_Comm_size(PETSC_COMM_WORLD,&size);
32: /* Read matrix and right-hand-side vector */
33: PetscOptionsGetString(NULL,NULL,"-f",file[0],sizeof(file[0]),&flg);
36: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&fd);
37: MatCreate(PETSC_COMM_WORLD,&A);
38: MatSetType(A,MATAIJ);
39: MatLoad(A,fd);
40: VecCreate(PETSC_COMM_WORLD,&b);
41: VecLoad(b,fd);
42: PetscViewerDestroy(&fd);
44: /*
45: If the loaded matrix is larger than the vector (due to being padded
46: to match the block size of the system), then create a new padded vector.
47: */
48: {
49: PetscInt m,n,j,mvec,start,end,indx;
50: Vec tmp;
51: PetscScalar *bold;
53: /* Create a new vector b by padding the old one */
54: MatGetLocalSize(A,&m,&n);
56: VecCreate(PETSC_COMM_WORLD,&tmp);
57: VecSetSizes(tmp,m,PETSC_DECIDE);
58: VecSetFromOptions(tmp);
59: VecGetOwnershipRange(b,&start,&end);
60: VecGetLocalSize(b,&mvec);
61: VecGetArray(b,&bold);
62: for (j=0; j<mvec; j++) {
63: indx = start+j;
64: VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);
65: }
66: VecRestoreArray(b,&bold);
67: VecDestroy(&b);
68: VecAssemblyBegin(tmp);
69: VecAssemblyEnd(tmp);
70: b = tmp;
71: }
72: VecDuplicate(b,&x);
73: VecDuplicate(b,&u);
74: VecSet(x,0.0);
76: /* Create dense matric B and X. Set B as an identity matrix */
77: MatGetSize(A,&M,&N);
78: MatCreate(MPI_COMM_SELF,&B);
79: MatSetSizes(B,M,N,M,N);
80: MatSetType(B,MATSEQDENSE);
81: MatSeqDenseSetPreallocation(B,NULL);
82: for (i=0; i<M; i++) {
83: MatSetValues(B,1,&i,1,&i,&val,INSERT_VALUES);
84: }
85: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
86: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
88: MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&X);
90: /* Compute X=inv(A) by MatMatSolve() */
91: KSPCreate(PETSC_COMM_WORLD,&ksp);
92: KSPSetOperators(ksp,A,A);
93: KSPGetPC(ksp,&pc);
94: PCSetType(pc,PCLU);
95: KSPSetFromOptions(ksp);
96: KSPSetUp(ksp);
97: PCFactorGetMatrix(pc,&F);
98: MatMatSolve(F,B,X);
99: MatDestroy(&B);
101: /* Now, set X=inv(A) as a preconditioner */
102: PCSetType(pc,PCSHELL);
103: PCShellSetContext(pc,X);
104: PCShellSetApply(pc,PCShellApply_Matinv);
105: KSPSetFromOptions(ksp);
107: /* Solve preconditioned system A*x = b */
108: KSPSolve(ksp,b,x);
109: KSPGetIterationNumber(ksp,&its);
111: /* Check error */
112: MatMult(A,x,u);
113: VecAXPY(u,-1.0,b);
114: VecNorm(u,NORM_2,&norm);
115: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
116: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);
118: /* Free work space. */
119: MatDestroy(&X);
120: MatDestroy(&A)); PetscCall(VecDestroy(&b);
121: VecDestroy(&u)); PetscCall(VecDestroy(&x);
122: KSPDestroy(&ksp);
123: PetscFinalize();
124: return 0;
125: }
127: PetscErrorCode PCShellApply_Matinv(PC pc,Vec xin,Vec xout)
128: {
129: Mat X;
132: PCShellGetContext(pc,&X);
133: MatMult(X,xin,xout);
134: return 0;
135: }
137: /*TEST
139: test:
140: args: -f ${DATAFILESPATH}/matrices/small
141: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
142: output_file: output/ex27.out
144: TEST*/