Actual source code: ex192.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  2: static char help[] = "Tests MatSolve() and MatMatSolve() with MUMPS or MKL_PARDISO sequential solvers in Schur complement mode.\n\
  3: Example: mpiexec -n 1 ./ex192 -f <matrix binary file> -nrhs 4 -symmetric_solve -hermitian_solve -schur_ratio 0.3\n\n";

  5: #include <petscmat.h>

  9: int main(int argc,char **args)
 10: {
 11:   Mat            A,RHS,C,F,X,S;
 12:   Vec            u,x,b;
 13:   Vec            xschur,bschur,uschur;
 14:   IS             is_schur;
 16:   PetscMPIInt    size;
 17:   PetscInt       isolver=0,size_schur,m,n,nfact,nsolve,nrhs;
 18:   PetscReal      norm,tol=PETSC_SQRT_MACHINE_EPSILON;
 19:   PetscRandom    rand;
 20:   PetscBool      data_provided,herm,symm,use_lu;
 21:   PetscReal      sratio = 5.1/12.;
 22:   PetscViewer    fd;              /* viewer */
 23:   char           solver[256];
 24:   char           file[PETSC_MAX_PATH_LEN]; /* input file name */

 26:   PetscInitialize(&argc,&args,(char*)0,help);
 27:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 28:   if (size > 1) SETERRQ(PETSC_COMM_WORLD,1,"This is a uniprocessor test");
 29:   /* Determine which type of solver we want to test for */
 30:   herm = PETSC_FALSE;
 31:   symm = PETSC_FALSE;
 32:   PetscOptionsGetBool(NULL,NULL,"-symmetric_solve",&symm,NULL);
 33:   PetscOptionsGetBool(NULL,NULL,"-hermitian_solve",&herm,NULL);
 34:   if (herm) symm = PETSC_TRUE;

 36:   /* Determine file from which we read the matrix A */
 37:   PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,&data_provided);
 38:   if (!data_provided) { /* get matrices from PETSc distribution */
 39:     sprintf(file,PETSC_DIR);
 40:     PetscStrcat(file,"/share/petsc/datafiles/matrices/");
 41:     if (symm) {
 42: #if defined (PETSC_USE_COMPLEX)
 43:       PetscStrcat(file,"hpd-complex-");
 44: #else
 45:       PetscStrcat(file,"spd-real-");
 46: #endif
 47:     } else {
 48: #if defined (PETSC_USE_COMPLEX)
 49:       PetscStrcat(file,"nh-complex-");
 50: #else
 51:       PetscStrcat(file,"ns-real-");
 52: #endif
 53:     }
 54: #if defined(PETSC_USE_64BIT_INDICES)
 55:     PetscStrcat(file,"int64-");
 56: #else
 57:     PetscStrcat(file,"int32-");
 58: #endif
 59: #if defined (PETSC_USE_REAL_SINGLE)
 60:     PetscStrcat(file,"float32");
 61: #else
 62:     PetscStrcat(file,"float64");
 63: #endif
 64:   }
 65:   /* Load matrix A */
 66:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
 67:   MatCreate(PETSC_COMM_WORLD,&A);
 68:   MatLoad(A,fd);
 69:   PetscViewerDestroy(&fd);
 70:   MatGetSize(A,&m,&n);
 71:   if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);

 73:   /* Create dense matrix C and X; C holds true solution with identical colums */
 74:   nrhs = 2;
 75:   PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);
 76:   MatCreate(PETSC_COMM_WORLD,&C);
 77:   MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
 78:   MatSetType(C,MATDENSE);
 79:   MatSetFromOptions(C);
 80:   MatSetUp(C);

 82:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
 83:   PetscRandomSetFromOptions(rand);
 84:   MatSetRandom(C,rand);
 85:   MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);

 87:   /* Create vectors */
 88:   VecCreate(PETSC_COMM_WORLD,&x);
 89:   VecSetSizes(x,n,PETSC_DECIDE);
 90:   VecSetFromOptions(x);
 91:   VecDuplicate(x,&b);
 92:   VecDuplicate(x,&u); /* save the true solution */

 94:   PetscOptionsGetInt(NULL,NULL,"-solver",&isolver,NULL);
 95:   switch (isolver) {
 96: #if defined(PETSC_HAVE_MUMPS)
 97:     case 0:
 98:       PetscStrcpy(solver,MATSOLVERMUMPS);
 99:       break;
100: #endif
101: #if defined(PETSC_HAVE_MKL_PARDISO)
102:     case 1:
103:       PetscStrcpy(solver,MATSOLVERMKL_PARDISO);
104:       break;
105: #endif
106:     default:
107:       PetscStrcpy(solver,MATSOLVERPETSC);
108:       break;
109:   }

111: #if defined (PETSC_USE_COMPLEX)
112:   if (isolver == 0 && symm && !data_provided) { /* MUMPS (5.0.0) does not have support for hermitian matrices, so make them symmetric */
113:     PetscScalar im = PetscSqrtScalar((PetscScalar)-1.);
114:     PetscScalar val = -1.0;
115:     val = val + im;
116:     MatSetValue(A,1,0,val,INSERT_VALUES);
117:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
118:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
119:   }
120: #endif

122:   PetscOptionsGetReal(NULL,NULL,"-schur_ratio",&sratio,NULL);
123:   if (sratio < 0. || sratio > 1.) {
124:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Invalid ratio for schur degrees of freedom %f", sratio);
125:   }
126:   size_schur = (PetscInt)(sratio*m);

128:   PetscPrintf(PETSC_COMM_SELF,"Solving with %s: nrhs %d, sym %d, herm %d, size schur %d, size mat %d\n",solver,nrhs,symm,herm,size_schur,m);

130:   /* Test LU/Cholesky Factorization */
131:   use_lu = PETSC_FALSE;
132:   if (!symm) use_lu = PETSC_TRUE;
133: #if defined (PETSC_USE_COMPLEX)
134:   if (isolver == 1) use_lu = PETSC_TRUE;
135: #endif

137:   if (herm && !use_lu) { /* test also conversion routines inside the solver packages */
138:     MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
139:     MatConvert(A,MATSEQSBAIJ,MAT_INPLACE_MATRIX,&A);
140:   }


143:   if (use_lu) {
144:     MatGetFactor(A,solver,MAT_FACTOR_LU,&F);
145:   } else {
146:     if (herm) {
147:       MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
148:       MatSetOption(A,MAT_SPD,PETSC_TRUE);
149:     } else {
150:       MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
151:       MatSetOption(A,MAT_SPD,PETSC_FALSE);
152:     }
153:     MatGetFactor(A,solver,MAT_FACTOR_CHOLESKY,&F);
154:   }
155:   ISCreateStride(PETSC_COMM_SELF,size_schur,m-size_schur,1,&is_schur);
156:   MatFactorSetSchurIS(F,is_schur);
157:   ISDestroy(&is_schur);
158:   if (use_lu) {
159:     MatLUFactorSymbolic(F,A,NULL,NULL,NULL);
160:   } else {
161:     MatCholeskyFactorSymbolic(F,A,NULL,NULL);
162:   }

164:   for (nfact = 0; nfact < 3; nfact++) {
165:     Mat AD;

167:     if (!nfact) {
168:       VecSetRandom(x,rand);
169:       if (symm && herm) {
170:         VecAbs(x);
171:       }
172:       MatDiagonalSet(A,x,ADD_VALUES);
173:     }
174:     if (use_lu) {
175:       MatLUFactorNumeric(F,A,NULL);
176:     } else {
177:       MatCholeskyFactorNumeric(F,A,NULL);
178:     }
179:     MatFactorCreateSchurComplement(F,&S);
180:     MatCreateVecs(S,&xschur,&bschur);
181:     VecDuplicate(xschur,&uschur);
182:     if (nfact == 1) {
183:       MatFactorInvertSchurComplement(F);
184:     }
185:     for (nsolve = 0; nsolve < 2; nsolve++) {
186:       VecSetRandom(x,rand);
187:       VecCopy(x,u);

189:       if (nsolve) {
190:         MatMult(A,x,b);
191:         MatSolve(F,b,x);
192:       } else {
193:         MatMultTranspose(A,x,b);
194:         MatSolveTranspose(F,b,x);
195:       }
196:       /* Check the error */
197:       VecAXPY(u,-1.0,x);  /* u <- (-1.0)x + u */
198:       VecNorm(u,NORM_2,&norm);
199:       if (norm > tol) {
200:         PetscReal resi;
201:         if (nsolve) {
202:           MatMult(A,x,u); /* u = A*x */
203:         } else {
204:           MatMultTranspose(A,x,u); /* u = A*x */
205:         }
206:         VecAXPY(u,-1.0,b);  /* u <- (-1.0)b + u */
207:         VecNorm(u,NORM_2,&resi);
208:         if (nsolve) {
209:           PetscPrintf(PETSC_COMM_SELF,"(f %d, s %d) MatSolve error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
210:         } else {
211:           PetscPrintf(PETSC_COMM_SELF,"(f %d, s %d) MatSolveTranspose error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
212:         }
213:       }
214:       VecSetRandom(xschur,rand);
215:       VecCopy(xschur,uschur);
216:       if (nsolve) {
217:         MatMult(S,xschur,bschur);
218:         MatFactorSolveSchurComplement(F,bschur,xschur);
219:       } else {
220:         MatMultTranspose(S,xschur,bschur);
221:         MatFactorSolveSchurComplementTranspose(F,bschur,xschur);
222:       }
223:       /* Check the error */
224:       VecAXPY(uschur,-1.0,xschur);  /* u <- (-1.0)x + u */
225:       VecNorm(uschur,NORM_2,&norm);
226:       if (norm > tol) {
227:         PetscReal resi;
228:         if (nsolve) {
229:           MatMult(S,xschur,uschur); /* u = A*x */
230:         } else {
231:           MatMultTranspose(S,xschur,uschur); /* u = A*x */
232:         }
233:         VecAXPY(uschur,-1.0,bschur);  /* u <- (-1.0)b + u */
234:         VecNorm(uschur,NORM_2,&resi);
235:         if (nsolve) {
236:           PetscPrintf(PETSC_COMM_SELF,"(f %d, s %d) MatFactorSolveSchurComplement error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
237:         } else {
238:           PetscPrintf(PETSC_COMM_SELF,"(f %d, s %d) MatFactorSolveSchurComplementTranspose error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
239:         }
240:       }
241:     }
242:     MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&AD);
243:     if (!nfact) {
244:       MatMatMult(AD,C,MAT_INITIAL_MATRIX,2.0,&RHS);
245:     } else {
246:       MatMatMult(AD,C,MAT_REUSE_MATRIX,2.0,&RHS);
247:     }
248:     MatDestroy(&AD);
249:     for (nsolve = 0; nsolve < 2; nsolve++) {
250:       MatMatSolve(F,RHS,X);

252:       /* Check the error */
253:       MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
254:       MatNorm(X,NORM_FROBENIUS,&norm);
255:       if (norm > tol) {
256:         PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatMatSolve: Norm of error %g\n",nfact,nsolve,norm);
257:       }
258:     }
259:     MatDestroy(&S);
260:     VecDestroy(&xschur);
261:     VecDestroy(&bschur);
262:     VecDestroy(&uschur);
263:   }
264:   /* Free data structures */
265:   MatDestroy(&A);
266:   MatDestroy(&C);
267:   MatDestroy(&F);
268:   MatDestroy(&X);
269:   MatDestroy(&RHS);
270:   PetscRandomDestroy(&rand);
271:   VecDestroy(&x);
272:   VecDestroy(&b);
273:   VecDestroy(&u);
274:   PetscFinalize();
275:   return 0;
276: }