Actual source code: ex27.c
petsc-3.13.4 2020-08-01
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;
23: PetscInt M,N,i,its;
24: PetscReal norm;
25: PetscScalar val=1.0;
26: PetscMPIInt size;
27: PC pc;
29: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
30: MPI_Comm_size(PETSC_COMM_WORLD,&size);
31: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
33: /* Read matrix and right-hand-side vector */
34: PetscOptionsGetString(NULL,NULL,"-f",file[0],PETSC_MAX_PATH_LEN,&flg);
35: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option");
37: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&fd);
38: MatCreate(PETSC_COMM_WORLD,&A);
39: MatSetType(A,MATAIJ);
40: MatLoad(A,fd);
41: VecCreate(PETSC_COMM_WORLD,&b);
42: VecLoad(b,fd);
43: PetscViewerDestroy(&fd);
45: /*
46: If the loaded matrix is larger than the vector (due to being padded
47: to match the block size of the system), then create a new padded vector.
48: */
49: {
50: PetscInt m,n,j,mvec,start,end,indx;
51: Vec tmp;
52: PetscScalar *bold;
54: /* Create a new vector b by padding the old one */
55: MatGetLocalSize(A,&m,&n);
56: if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%D, %D)", m, n);
57: VecCreate(PETSC_COMM_WORLD,&tmp);
58: VecSetSizes(tmp,m,PETSC_DECIDE);
59: VecSetFromOptions(tmp);
60: VecGetOwnershipRange(b,&start,&end);
61: VecGetLocalSize(b,&mvec);
62: VecGetArray(b,&bold);
63: for (j=0; j<mvec; j++) {
64: indx = start+j;
65: VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);
66: }
67: VecRestoreArray(b,&bold);
68: VecDestroy(&b);
69: VecAssemblyBegin(tmp);
70: VecAssemblyEnd(tmp);
71: b = tmp;
72: }
73: VecDuplicate(b,&x);
74: VecDuplicate(b,&u);
75: VecSet(x,0.0);
77: /* Create dense matric B and X. Set B as an identity matrix */
78: MatGetSize(A,&M,&N);
79: MatCreate(MPI_COMM_SELF,&B);
80: MatSetSizes(B,M,N,M,N);
81: MatSetType(B,MATSEQDENSE);
82: MatSeqDenseSetPreallocation(B,NULL);
83: for (i=0; i<M; i++) {
84: MatSetValues(B,1,&i,1,&i,&val,INSERT_VALUES);
85: }
86: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
87: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
89: MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&X);
91: /* Compute X=inv(A) by MatMatSolve() */
92: KSPCreate(PETSC_COMM_WORLD,&ksp);
93: KSPSetOperators(ksp,A,A);
94: KSPGetPC(ksp,&pc);
95: PCSetType(pc,PCLU);
96: KSPSetFromOptions(ksp);
97: KSPSetUp(ksp);
98: PCFactorGetMatrix(pc,&F);
99: MatMatSolve(F,B,X);
100: MatDestroy(&B);
102: /* Now, set X=inv(A) as a preconditioner */
103: PCSetType(pc,PCSHELL);
104: PCShellSetContext(pc,(void*)X);
105: PCShellSetApply(pc,PCShellApply_Matinv);
106: KSPSetFromOptions(ksp);
108: /* Solve preconditioned system A*x = b */
109: KSPSolve(ksp,b,x);
110: KSPGetIterationNumber(ksp,&its);
112: /* Check error */
113: MatMult(A,x,u);
114: VecAXPY(u,-1.0,b);
115: VecNorm(u,NORM_2,&norm);
116: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
117: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);
119: /* Free work space. */
120: MatDestroy(&X);
121: MatDestroy(&A); VecDestroy(&b);
122: VecDestroy(&u); VecDestroy(&x);
123: KSPDestroy(&ksp);
124: PetscFinalize();
125: return ierr;
126: }
128: PetscErrorCode PCShellApply_Matinv(PC pc,Vec xin,Vec xout)
129: {
131: Mat X;
134: PCShellGetContext(pc,(void**)&X);
135: MatMult(X,xin,xout);
136: return(0);
137: }
139: /*TEST
141: test:
142: args: -f ${DATAFILESPATH}/matrices/small
143: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
144: output_file: output/ex27.out
146: TEST*/