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 = NULL, RHS1 = NULL, 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, testMatMatSolveTranspose = 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 */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

231:         /* Check the error */
232:         MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN);
233:         MatNorm(X, NORM_FROBENIUS, &norm);
234:         if (norm > tol) PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the MatMatSolve: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve);
235:       }

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) PetscPrintf(PETSC_COMM_WORLD, "MatMatSolve(F,RHS,RHS): Norm of error %g\n", (double)norm);
245:       }

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

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

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

268:     /* Test testMatMatSolveTranspose(), A^T X = B, where B can be dense or sparse */
269:     if (testMatMatSolveTranspose) {
270:       if (!nfact) {
271:         MatTransposeMatMult(A, C, MAT_INITIAL_MATRIX, 2.0, &RHS1);
272:       } else {
273:         MatTransposeMatMult(A, C, MAT_REUSE_MATRIX, 2.0, &RHS1);
274:       }

276:       for (nsolve = 0; nsolve < 2; nsolve++) {
277:         MatMatSolveTranspose(F, RHS1, X);

279:         /* Check the error */
280:         MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN);
281:         MatNorm(X, NORM_FROBENIUS, &norm);
282:         if (norm > tol) PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the MatMatSolveTranspose: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve);
283:       }

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

288:         MatTranspose(RHS1, MAT_INITIAL_MATRIX, &RHST);
289:         MatConvert(RHST, MATAIJ, MAT_INITIAL_MATRIX, &spRHST);
290:         MatCreateTranspose(spRHST, &spRHS);
291:         for (nsolve = 0; nsolve < 2; nsolve++) {
292:           MatMatSolveTranspose(F, spRHS, X);

294:           /* Check the error */
295:           MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN);
296:           MatNorm(X, NORM_FROBENIUS, &norm);
297:           if (norm > tol) PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the sparse MatMatSolveTranspose: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve);
298:         }
299:         MatDestroy(&spRHST);
300:         MatDestroy(&spRHS);
301:         MatDestroy(&RHST);
302:       }
303:     }

305:     /* Test MatSolve() */
306:     if (testMatSolve) {
307:       for (nsolve = 0; nsolve < 2; nsolve++) {
308:         VecSetRandom(x, rand);
309:         VecCopy(x, u);
310:         MatMult(A, x, b);

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

315:         /* Check the error */
316:         VecAXPY(u, -1.0, x); /* u <- (-1.0)x + u */
317:         VecNorm(u, NORM_2, &norm);
318:         if (norm > tol) {
319:           PetscReal resi;
320:           MatMult(A, x, u);    /* u = A*x */
321:           VecAXPY(u, -1.0, b); /* u <- (-1.0)b + u */
322:           VecNorm(u, NORM_2, &resi);
323:           PetscPrintf(PETSC_COMM_WORLD, "MatSolve: Norm of error %g, resi %g, numfact %" PetscInt_FMT "\n", (double)norm, (double)resi, nfact);
324:         }
325:       }
326:     }
327:   }

329:   /* Free data structures */
330:   MatDestroy(&A);
331:   MatDestroy(&C);
332:   MatDestroy(&F);
333:   MatDestroy(&X);
334:   MatDestroy(&RHS);
335:   MatDestroy(&RHS1);

337:   PetscRandomDestroy(&rand);
338:   ISDestroy(&perm);
339:   ISDestroy(&iperm);
340:   VecDestroy(&x);
341:   VecDestroy(&b);
342:   VecDestroy(&u);
343:   PetscFinalize();
344:   return 0;
345: }

347: /*TEST

349:    test:
350:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
351:       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 10
352:       output_file: output/ex125.out

354:    test:
355:       suffix: 2
356:       args: -mat_solver_type 10
357:       output_file: output/ex125.out

359:    test:
360:       suffix: mkl_pardiso
361:       requires: mkl_pardiso datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
362:       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 3

364:    test:
365:       suffix: mkl_pardiso_2
366:       requires: mkl_pardiso
367:       args: -mat_solver_type 3
368:       output_file: output/ex125_mkl_pardiso.out

370:    test:
371:       suffix: mumps
372:       requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
373:       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2
374:       output_file: output/ex125_mumps_seq.out

376:    test:
377:       suffix: mumps_2
378:       nsize: 3
379:       requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
380:       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2
381:       output_file: output/ex125_mumps_par.out

383:    test:
384:       suffix: mumps_3
385:       requires: mumps
386:       args: -mat_solver_type 2
387:       output_file: output/ex125_mumps_seq.out

389:    test:
390:       suffix: mumps_4
391:       nsize: 3
392:       requires: mumps
393:       args: -mat_solver_type 2
394:       output_file: output/ex125_mumps_par.out

396:    test:
397:       suffix: mumps_5
398:       nsize: 3
399:       requires: mumps
400:       args: -mat_solver_type 2 -cholesky
401:       output_file: output/ex125_mumps_par_cholesky.out

403:    test:
404:       suffix: superlu_dist
405:       nsize: {{1 3}}
406:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu_dist
407:       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 1 -mat_superlu_dist_rowperm NOROWPERM

409:    test:
410:       suffix: superlu_dist_2
411:       nsize: {{1 3}}
412:       requires: superlu_dist !complex
413:       args: -n 36 -mat_solver_type 1 -mat_superlu_dist_rowperm NOROWPERM
414:       output_file: output/ex125_superlu_dist.out

416:    test:
417:       suffix: superlu_dist_complex
418:       nsize: 3
419:       requires: datafilespath double superlu_dist complex !defined(PETSC_USE_64BIT_INDICES)
420:       args: -f ${DATAFILESPATH}/matrices/farzad_B_rhs -mat_solver_type 1
421:       output_file: output/ex125_superlu_dist_complex.out

423:    test:
424:       suffix: superlu_dist_complex_2
425:       nsize: 3
426:       requires: superlu_dist complex
427:       args: -mat_solver_type 1
428:       output_file: output/ex125_superlu_dist_complex.out

430:    test:
431:       suffix: cusparse
432:       requires: cuda datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
433:       #todo: fix the bug with cholesky
434:       #args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type 4 -cholesky {{0 1}separate output}
435:       args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type 4 -cholesky {{0}separate output}

437:    test:
438:       suffix: cusparse_2
439:       requires: cuda
440:       args: -mat_type aijcusparse -mat_solver_type 4 -cholesky {{0 1}separate output}

442: TEST*/