Actual source code: ex125.c

  1: static char help[] = "Tests MatSolve() and MatMatSolve() (interface to superlu_dist, mumps and mkl_pardiso).\n\
  2: Example: mpiexec -n <np> ./ex125 -f <matrix binary file> -nrhs 4 \n\n";

  4: #include <petscmat.h>

  6: int main(int argc,char **args)
  7: {
  8:   Mat            A,RHS,C,F,X;
  9:   Vec            u,x,b;
 10:   PetscMPIInt    size;
 11:   PetscInt       m,n,nfact,nsolve,nrhs,ipack=0;
 12:   PetscReal      norm,tol=1.e-10;
 13:   IS             perm,iperm;
 14:   MatFactorInfo  info;
 15:   PetscRandom    rand;
 16:   PetscBool      flg,testMatSolve=PETSC_TRUE,testMatMatSolve=PETSC_TRUE;
 17:   PetscBool      chol=PETSC_FALSE,view=PETSC_FALSE,matsolvexx = PETSC_FALSE;
 18: #if defined(PETSC_HAVE_MUMPS)
 19:   PetscBool      test_mumps_opts=PETSC_FALSE;
 20: #endif
 21:   PetscViewer    fd;              /* viewer */
 22:   char           file[PETSC_MAX_PATH_LEN]; /* input file name */

 24:   PetscInitialize(&argc,&args,(char*)0,help);
 25:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

 27:   /* Determine file from which we read the matrix A */
 28:   PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg);
 29:   if (flg) { /* Load matrix A */
 30:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
 31:     MatCreate(PETSC_COMM_WORLD,&A);
 32:     MatSetFromOptions(A);
 33:     MatLoad(A,fd);
 34:     PetscViewerDestroy(&fd);
 35:   } else {
 36:     n = 13;
 37:     PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 38:     MatCreate(PETSC_COMM_WORLD,&A);
 39:     MatSetType(A,MATAIJ);
 40:     MatSetFromOptions(A);
 41:     MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
 42:     MatSetUp(A);
 43:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 44:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 45:     MatShift(A,1.0);
 46:   }
 47:   MatGetLocalSize(A,&m,&n);

 50:   /* if A is symmetric, set its flag -- required by MatGetInertia() */
 51:   MatIsSymmetric(A,0.0,&flg);

 53:   MatViewFromOptions(A,NULL,"-A_view");

 55:   /* Create dense matrix C and X; C holds true solution with identical columns */
 56:   nrhs = 2;
 57:   PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);
 58:   PetscPrintf(PETSC_COMM_WORLD,"ex125: nrhs %" PetscInt_FMT "\n",nrhs);
 59:   MatCreate(PETSC_COMM_WORLD,&C);
 60:   MatSetOptionsPrefix(C,"rhs_");
 61:   MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
 62:   MatSetType(C,MATDENSE);
 63:   MatSetFromOptions(C);
 64:   MatSetUp(C);

 66:   PetscOptionsGetBool(NULL,NULL,"-view_factor",&view,NULL);
 67:   PetscOptionsGetBool(NULL,NULL,"-test_matmatsolve",&testMatMatSolve,NULL);
 68:   PetscOptionsGetBool(NULL,NULL,"-cholesky",&chol,NULL);
 69: #if defined(PETSC_HAVE_MUMPS)
 70:   PetscOptionsGetBool(NULL,NULL,"-test_mumps_opts",&test_mumps_opts,NULL);
 71: #endif

 73:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
 74:   PetscRandomSetFromOptions(rand);
 75:   MatSetRandom(C,rand);
 76:   MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);

 78:   /* Create vectors */
 79:   MatCreateVecs(A,&x,&b);
 80:   VecDuplicate(x,&u); /* save the true solution */

 82:   /* Test Factorization */
 83:   MatGetOrdering(A,MATORDERINGND,&perm,&iperm);

 85:   PetscOptionsGetInt(NULL,NULL,"-mat_solver_type",&ipack,NULL);
 86:   switch (ipack) {
 87: #if defined(PETSC_HAVE_SUPERLU)
 88:   case 0:
 90:     PetscPrintf(PETSC_COMM_WORLD," SUPERLU LU:\n");
 91:     MatGetFactor(A,MATSOLVERSUPERLU,MAT_FACTOR_LU,&F);
 92:     matsolvexx = PETSC_TRUE;
 93:     break;
 94: #endif
 95: #if defined(PETSC_HAVE_SUPERLU_DIST)
 96:   case 1:
 98:     PetscPrintf(PETSC_COMM_WORLD," SUPERLU_DIST LU:\n");
 99:     MatGetFactor(A,MATSOLVERSUPERLU_DIST,MAT_FACTOR_LU,&F);
100:     matsolvexx = PETSC_TRUE;
101:     break;
102: #endif
103: #if defined(PETSC_HAVE_MUMPS)
104:   case 2:
105:     if (chol) {
106:       PetscPrintf(PETSC_COMM_WORLD," MUMPS CHOLESKY:\n");
107:       MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_CHOLESKY,&F);
108:     } else {
109:       PetscPrintf(PETSC_COMM_WORLD," MUMPS LU:\n");
110:       MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,&F);
111:     }
112:     matsolvexx = PETSC_TRUE;
113:     if (test_mumps_opts) {
114:       /* test mumps options */
115:       PetscInt  icntl;
116:       PetscReal cntl;

118:       icntl = 2;        /* sequential matrix ordering */
119:       MatMumpsSetIcntl(F,7,icntl);

121:       cntl = 1.e-6; /* threshold for row pivot detection */
122:       MatMumpsSetIcntl(F,24,1);
123:       MatMumpsSetCntl(F,3,cntl);
124:     }
125:     break;
126: #endif
127: #if defined(PETSC_HAVE_MKL_PARDISO)
128:   case 3:
129:     if (chol) {
130:       PetscPrintf(PETSC_COMM_WORLD," MKL_PARDISO CHOLESKY:\n");
131:       MatGetFactor(A,MATSOLVERMKL_PARDISO,MAT_FACTOR_CHOLESKY,&F);
132:     } else {
133:       PetscPrintf(PETSC_COMM_WORLD," MKL_PARDISO LU:\n");
134:       MatGetFactor(A,MATSOLVERMKL_PARDISO,MAT_FACTOR_LU,&F);
135:     }
136:     break;
137: #endif
138: #if defined(PETSC_HAVE_CUDA)
139:   case 4:
140:     if (chol) {
141:       PetscPrintf(PETSC_COMM_WORLD," CUSPARSE CHOLESKY:\n");
142:       MatGetFactor(A,MATSOLVERCUSPARSE,MAT_FACTOR_CHOLESKY,&F);
143:     } else {
144:       PetscPrintf(PETSC_COMM_WORLD," CUSPARSE LU:\n");
145:       MatGetFactor(A,MATSOLVERCUSPARSE,MAT_FACTOR_LU,&F);
146:     }
147:     break;
148: #endif
149:   default:
150:     if (chol) {
151:       PetscPrintf(PETSC_COMM_WORLD," PETSC CHOLESKY:\n");
152:       MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&F);
153:     } else {
154:       PetscPrintf(PETSC_COMM_WORLD," PETSC LU:\n");
155:       MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
156:     }
157:     matsolvexx = PETSC_TRUE;
158:   }

160:   MatFactorInfoInitialize(&info);
161:   info.fill      = 5.0;
162:   info.shifttype = (PetscReal) MAT_SHIFT_NONE;
163:   if (chol) {
164:     MatCholeskyFactorSymbolic(F,A,perm,&info);
165:   } else {
166:     MatLUFactorSymbolic(F,A,perm,iperm,&info);
167:   }

169:   for (nfact = 0; nfact < 2; nfact++) {
170:     if (chol) {
171:       PetscPrintf(PETSC_COMM_WORLD," %" PetscInt_FMT "-the CHOLESKY numfactorization \n",nfact);
172:       MatCholeskyFactorNumeric(F,A,&info);
173:     } else {
174:       PetscPrintf(PETSC_COMM_WORLD," %" PetscInt_FMT "-the LU numfactorization \n",nfact);
175:       MatLUFactorNumeric(F,A,&info);
176:     }
177:     if (view) {
178:       PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO);
179:       MatView(F,PETSC_VIEWER_STDOUT_WORLD);
180:       PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
181:       view = PETSC_FALSE;
182:     }

184: #if defined(PETSC_HAVE_SUPERLU_DIST)
185:     if (ipack == 1) { /* Test MatSuperluDistGetDiagU()
186:        -- input: matrix factor F; output: main diagonal of matrix U on all processes */
187:       PetscInt    M;
188:       PetscScalar *diag;
189: #if !defined(PETSC_USE_COMPLEX)
190:       PetscInt nneg,nzero,npos;
191: #endif

193:       MatGetSize(F,&M,NULL);
194:       PetscMalloc1(M,&diag);
195:       MatSuperluDistGetDiagU(F,diag);
196:       PetscFree(diag);

198: #if !defined(PETSC_USE_COMPLEX)
199:       /* Test MatGetInertia() */
200:       MatGetInertia(F,&nneg,&nzero,&npos);
201:       PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD," MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n",nneg,nzero,npos);
202: #endif
203:     }
204: #endif

206: #if defined(PETSC_HAVE_MUMPS)
207:     /* mumps interface allows repeated call of MatCholeskyFactorSymbolic(), while the succession calls do nothing */
208:     if (ipack == 2) {
209:       if (chol) {
210:         MatCholeskyFactorSymbolic(F,A,perm,&info);
211:         MatCholeskyFactorNumeric(F,A,&info);
212:       } else {
213:         MatLUFactorSymbolic(F,A,perm,iperm,&info);
214:         MatLUFactorNumeric(F,A,&info);
215:       }
216:     }
217: #endif

219:     /* Test MatMatSolve() */
220:     if (testMatMatSolve) {
221:       if (!nfact) {
222:         MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS);
223:       } else {
224:         MatMatMult(A,C,MAT_REUSE_MATRIX,2.0,&RHS);
225:       }
226:       for (nsolve = 0; nsolve < 2; nsolve++) {
227:         PetscPrintf(PETSC_COMM_WORLD,"   %" PetscInt_FMT "-the MatMatSolve \n",nsolve);
228:         MatMatSolve(F,RHS,X);

230:         /* Check the error */
231:         MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
232:         MatNorm(X,NORM_FROBENIUS,&norm);
233:         if (norm > tol) {
234:           PetscPrintf(PETSC_COMM_WORLD,"%" PetscInt_FMT "-the MatMatSolve: Norm of error %g, nsolve %" PetscInt_FMT "\n",nsolve,(double)norm,nsolve);
235:         }
236:       }
237:       if (matsolvexx) {
238:         /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix */
239:         MatCopy(RHS,X,SAME_NONZERO_PATTERN);
240:         MatMatSolve(F,X,X);
241:         /* Check the error */
242:         MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
243:         MatNorm(X,NORM_FROBENIUS,&norm);
244:         if (norm > tol) {
245:           PetscPrintf(PETSC_COMM_WORLD,"MatMatSolve(F,RHS,RHS): Norm of error %g\n",(double)norm);
246:         }
247:       }

249:       if (ipack == 2 && size == 1) {
250:         Mat spRHS,spRHST,RHST;

252:         MatTranspose(RHS,MAT_INITIAL_MATRIX,&RHST);
253:         MatConvert(RHST,MATAIJ,MAT_INITIAL_MATRIX,&spRHST);
254:         MatCreateTranspose(spRHST,&spRHS);
255:         for (nsolve = 0; nsolve < 2; nsolve++) {
256:           PetscPrintf(PETSC_COMM_WORLD,"   %" PetscInt_FMT "-the sparse MatMatSolve \n",nsolve);
257:           MatMatSolve(F,spRHS,X);

259:           /* Check the error */
260:           MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
261:           MatNorm(X,NORM_FROBENIUS,&norm);
262:           if (norm > tol) {
263:             PetscPrintf(PETSC_COMM_WORLD,"%" PetscInt_FMT "-the sparse MatMatSolve: Norm of error %g, nsolve %" PetscInt_FMT "\n",nsolve,(double)norm,nsolve);
264:           }
265:         }
266:         MatDestroy(&spRHST);
267:         MatDestroy(&spRHS);
268:         MatDestroy(&RHST);
269:       }
270:     }

272:     /* Test MatSolve() */
273:     if (testMatSolve) {
274:       for (nsolve = 0; nsolve < 2; nsolve++) {
275:         VecSetRandom(x,rand);
276:         VecCopy(x,u);
277:         MatMult(A,x,b);

279:         PetscPrintf(PETSC_COMM_WORLD,"   %" PetscInt_FMT "-the MatSolve \n",nsolve);
280:         MatSolve(F,b,x);

282:         /* Check the error */
283:         VecAXPY(u,-1.0,x);  /* u <- (-1.0)x + u */
284:         VecNorm(u,NORM_2,&norm);
285:         if (norm > tol) {
286:           PetscReal resi;
287:           MatMult(A,x,u); /* u = A*x */
288:           VecAXPY(u,-1.0,b);  /* u <- (-1.0)b + u */
289:           VecNorm(u,NORM_2,&resi);
290:           PetscPrintf(PETSC_COMM_WORLD,"MatSolve: Norm of error %g, resi %g, numfact %" PetscInt_FMT "\n",(double)norm,(double)resi,nfact);
291:         }
292:       }
293:     }
294:   }

296:   /* Free data structures */
297:   MatDestroy(&A);
298:   MatDestroy(&C);
299:   MatDestroy(&F);
300:   MatDestroy(&X);
301:   if (testMatMatSolve) {
302:     MatDestroy(&RHS);
303:   }

305:   PetscRandomDestroy(&rand);
306:   ISDestroy(&perm);
307:   ISDestroy(&iperm);
308:   VecDestroy(&x);
309:   VecDestroy(&b);
310:   VecDestroy(&u);
311:   PetscFinalize();
312:   return 0;
313: }

315: /*TEST

317:    test:
318:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
319:       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 10
320:       output_file: output/ex125.out

322:    test:
323:       suffix: 2
324:       args: -mat_solver_type 10
325:       output_file: output/ex125.out

327:    test:
328:       suffix: mkl_pardiso
329:       requires: mkl_pardiso datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
330:       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 3

332:    test:
333:       suffix: mkl_pardiso_2
334:       requires: mkl_pardiso
335:       args: -mat_solver_type 3
336:       output_file: output/ex125_mkl_pardiso.out

338:    test:
339:       suffix: mumps
340:       requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
341:       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2
342:       output_file: output/ex125_mumps_seq.out

344:    test:
345:       suffix: mumps_2
346:       nsize: 3
347:       requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
348:       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2
349:       output_file: output/ex125_mumps_par.out

351:    test:
352:       suffix: mumps_3
353:       requires: mumps
354:       args: -mat_solver_type 2
355:       output_file: output/ex125_mumps_seq.out

357:    test:
358:       suffix: mumps_4
359:       nsize: 3
360:       requires: mumps
361:       args: -mat_solver_type 2
362:       output_file: output/ex125_mumps_par.out

364:    test:
365:       suffix: mumps_5
366:       nsize: 3
367:       requires: mumps
368:       args: -mat_solver_type 2 -cholesky
369:       output_file: output/ex125_mumps_par_cholesky.out

371:    test:
372:       suffix: superlu_dist
373:       nsize: {{1 3}}
374:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) superlu_dist
375:       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 1 -mat_superlu_dist_rowperm NOROWPERM

377:    test:
378:       suffix: superlu_dist_2
379:       nsize: {{1 3}}
380:       requires: superlu_dist !complex
381:       args: -n 36 -mat_solver_type 1 -mat_superlu_dist_rowperm NOROWPERM
382:       output_file: output/ex125_superlu_dist.out

384:    test:
385:       suffix: superlu_dist_complex
386:       nsize: 3
387:       requires: datafilespath superlu_dist complex double !defined(PETSC_USE_64BIT_INDICES)
388:       args: -f ${DATAFILESPATH}/matrices/farzad_B_rhs -mat_solver_type 1
389:       output_file: output/ex125_superlu_dist_complex.out

391:    test:
392:       suffix: superlu_dist_complex_2
393:       nsize: 3
394:       requires: superlu_dist complex
395:       args: -mat_solver_type 1
396:       output_file: output/ex125_superlu_dist_complex.out

398:    test:
399:       suffix: cusparse
400:       requires: cuda datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
401:       args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type 4 -cholesky {{0 1}separate output}

403:    test:
404:       suffix: cusparse_2
405:       requires: cuda
406:       args: -mat_type aijcusparse -mat_solver_type 4 -cholesky {{0 1}separate output}

408: TEST*/