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 -mat_solver_type <>\n\n";

  4: /*
  5: -mat_solver_type:
  6:   superlu
  7:   superlu_dist
  8:   mumps
  9:   mkl_pardiso
 10:   cusparse
 11:   petsc
 12: */

 14: #include <petscmat.h>

 16: PetscErrorCode CreateRandom(PetscInt n, PetscInt m, Mat *A)
 17: {
 18:   PetscFunctionBeginUser;
 19:   PetscCall(MatCreate(PETSC_COMM_WORLD, A));
 20:   PetscCall(MatSetType(*A, MATAIJ));
 21:   PetscCall(MatSetFromOptions(*A));
 22:   PetscCall(MatSetSizes(*A, PETSC_DECIDE, PETSC_DECIDE, n, m));
 23:   PetscCall(MatSeqAIJSetPreallocation(*A, 5, NULL));
 24:   PetscCall(MatMPIAIJSetPreallocation(*A, 5, NULL, 5, NULL));
 25:   PetscCall(MatSetRandom(*A, NULL));
 26:   PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
 27:   PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
 28:   PetscFunctionReturn(PETSC_SUCCESS);
 29: }

 31: PetscErrorCode CreateIdentity(PetscInt n, Mat *A)
 32: {
 33:   PetscFunctionBeginUser;
 34:   PetscCall(MatCreate(PETSC_COMM_WORLD, A));
 35:   PetscCall(MatSetType(*A, MATAIJ));
 36:   PetscCall(MatSetFromOptions(*A));
 37:   PetscCall(MatSetSizes(*A, PETSC_DECIDE, PETSC_DECIDE, n, n));
 38:   PetscCall(MatSetUp(*A));
 39:   PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
 40:   PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
 41:   PetscCall(MatShift(*A, 1.0));
 42:   PetscFunctionReturn(PETSC_SUCCESS);
 43: }

 45: int main(int argc, char **args)
 46: {
 47:   Mat           A, Ae, RHS = NULL, RHS1 = NULL, C, F, X;
 48:   Vec           u, x, b;
 49:   PetscMPIInt   size;
 50:   PetscInt      m, n, nfact, nsolve, nrhs, ipack = 5;
 51:   PetscReal     norm, tol = 10 * PETSC_SQRT_MACHINE_EPSILON;
 52:   IS            perm = NULL, iperm = NULL;
 53:   MatFactorInfo info;
 54:   PetscRandom   rand;
 55:   PetscBool     flg, symm, testMatSolve = PETSC_TRUE, testMatMatSolve = PETSC_TRUE, testMatMatSolveTranspose = PETSC_TRUE, testMatSolveTranspose = PETSC_TRUE, match = PETSC_FALSE;
 56:   PetscBool     chol = PETSC_FALSE, view = PETSC_FALSE, matsolvexx = PETSC_FALSE;
 57: #if defined(PETSC_HAVE_MUMPS)
 58:   PetscBool test_mumps_opts = PETSC_FALSE;
 59: #endif
 60:   PetscViewer fd;                       /* viewer */
 61:   char        file[PETSC_MAX_PATH_LEN]; /* input file name */
 62:   char        pack[PETSC_MAX_PATH_LEN];

 64:   PetscFunctionBeginUser;
 65:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 66:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));

 68:   /* Determine file from which we read the matrix A */
 69:   PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg));
 70:   if (flg) { /* Load matrix A */
 71:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
 72:     PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 73:     PetscCall(MatSetFromOptions(A));
 74:     PetscCall(MatLoad(A, fd));
 75:     PetscCall(PetscViewerDestroy(&fd));
 76:   } else {
 77:     n = 13;
 78:     PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
 79:     PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 80:     PetscCall(MatSetType(A, MATAIJ));
 81:     PetscCall(MatSetFromOptions(A));
 82:     PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
 83:     PetscCall(MatSetUp(A));
 84:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 85:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 86:     PetscCall(MatShift(A, 1.0));
 87:   }

 89:   /* if A is symmetric, set its flag -- required by MatGetInertia() */
 90:   PetscCall(MatIsSymmetric(A, 0.0, &symm));
 91:   PetscCall(MatSetOption(A, MAT_SYMMETRIC, symm));

 93:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-cholesky", &chol, NULL));

 95:   /* test MATNEST support */
 96:   flg = PETSC_FALSE;
 97:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_nest", &flg, NULL));
 98:   if (flg) {
 99:     Mat B;

101:     flg = PETSC_FALSE;
102:     PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_nest_bordered", &flg, NULL));
103:     if (!flg) {
104:       Mat mats[9] = {NULL, NULL, A, NULL, A, NULL, A, NULL, NULL};

106:       /* Create a nested matrix representing
107:               | 0 0 A |
108:               | 0 A 0 |
109:               | A 0 0 |
110:       */
111:       PetscCall(MatCreateNest(PETSC_COMM_WORLD, 3, NULL, 3, NULL, mats, &B));
112:     } else {
113:       Mat mats[4];

115:       /* Create a nested matrix representing
116:               | Id  R |
117:               | R^t A |
118:       */
119:       PetscCall(MatGetSize(A, NULL, &n));
120:       m = n + 12;
121:       PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
122:       PetscCall(CreateIdentity(m, &mats[0]));
123:       PetscCall(CreateRandom(m, n, &mats[1]));
124:       mats[3] = A;

126:       /* use CreateTranspose/CreateHermitianTranspose or explicit matrix for debugging purposes */
127:       flg = PETSC_FALSE;
128:       PetscCall(PetscOptionsGetBool(NULL, NULL, "-expl", &flg, NULL));
129:       if (PetscDefined(USE_COMPLEX)) {
130:         if (chol) { /* Hermitian transpose not supported by MUMPS Cholesky factor */
131:           if (!flg) PetscCall(MatCreateTranspose(mats[1], &mats[2]));
132:           else PetscCall(MatTranspose(mats[1], MAT_INITIAL_MATRIX, &mats[2]));
133:         } else {
134:           if (!flg) PetscCall(MatCreateHermitianTranspose(mats[1], &mats[2]));
135:           else PetscCall(MatHermitianTranspose(mats[1], MAT_INITIAL_MATRIX, &mats[2]));
136:         }
137:       } else {
138:         if (!flg) PetscCall(MatCreateTranspose(mats[1], &mats[2]));
139:         else PetscCall(MatTranspose(mats[1], MAT_INITIAL_MATRIX, &mats[2]));
140:       }
141:       PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, mats, &B));
142:       PetscCall(MatDestroy(&mats[0]));
143:       PetscCall(MatDestroy(&mats[1]));
144:       PetscCall(MatDestroy(&mats[2]));
145:     }
146:     PetscCall(MatDestroy(&A));
147:     A = B;
148:     PetscCall(MatSetOption(A, MAT_SYMMETRIC, symm));

150:     /* not all the combinations of MatMat operations are supported by MATNEST. */
151:     PetscCall(MatComputeOperator(A, MATAIJ, &Ae));
152:   } else {
153:     PetscCall(PetscObjectReference((PetscObject)A));
154:     Ae = A;
155:   }
156:   PetscCall(MatGetLocalSize(A, &m, &n));
157:   PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")", m, n);

159:   PetscCall(MatViewFromOptions(A, NULL, "-A_view"));
160:   PetscCall(MatViewFromOptions(Ae, NULL, "-A_view_expl"));

162:   /* Create dense matrix C and X; C holds true solution with identical columns */
163:   nrhs = 2;
164:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL));
165:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ex125: nrhs %" PetscInt_FMT "\n", nrhs));
166:   PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
167:   PetscCall(MatSetOptionsPrefix(C, "rhs_"));
168:   PetscCall(MatSetSizes(C, m, PETSC_DECIDE, PETSC_DECIDE, nrhs));
169:   PetscCall(MatSetType(C, MATDENSE));
170:   PetscCall(MatSetFromOptions(C));
171:   PetscCall(MatSetUp(C));

173:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_factor", &view, NULL));
174:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matmatsolve", &testMatMatSolve, NULL));
175:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matmatsolvetranspose", &testMatMatSolveTranspose, NULL));
176:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matsolvetranspose", &testMatSolveTranspose, NULL));
177: #if defined(PETSC_HAVE_MUMPS)
178:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_mumps_opts", &test_mumps_opts, NULL));
179: #endif

181:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
182:   PetscCall(PetscRandomSetFromOptions(rand));
183:   PetscCall(MatSetRandom(C, rand));
184:   PetscCall(MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &X));

186:   /* Create vectors */
187:   PetscCall(MatCreateVecs(A, &x, &b));
188:   PetscCall(VecDuplicate(x, &u)); /* save the true solution */

190:   /* Test Factorization */
191:   PetscCall(MatGetOrdering(A, MATORDERINGND, &perm, &iperm));

193:   PetscCall(PetscOptionsGetString(NULL, NULL, "-mat_solver_type", pack, sizeof(pack), NULL));
194: #if defined(PETSC_HAVE_SUPERLU)
195:   PetscCall(PetscStrcmp(MATSOLVERSUPERLU, pack, &match));
196:   if (match) {
197:     PetscCheck(!chol, PETSC_COMM_WORLD, PETSC_ERR_SUP, "SuperLU does not provide Cholesky!");
198:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, " SUPERLU LU:\n"));
199:     PetscCall(MatGetFactor(A, MATSOLVERSUPERLU, MAT_FACTOR_LU, &F));
200:     matsolvexx = PETSC_FALSE; /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix, need further work */
201:     ipack      = 0;
202:     goto skipoptions;
203:   }
204: #endif
205: #if defined(PETSC_HAVE_SUPERLU_DIST)
206:   PetscCall(PetscStrcmp(MATSOLVERSUPERLU_DIST, pack, &match));
207:   if (match) {
208:     PetscCheck(!chol, PETSC_COMM_WORLD, PETSC_ERR_SUP, "SuperLU does not provide Cholesky!");
209:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, " SUPERLU_DIST LU:\n"));
210:     PetscCall(MatGetFactor(A, MATSOLVERSUPERLU_DIST, MAT_FACTOR_LU, &F));
211:     matsolvexx = PETSC_TRUE;
212:     if (symm) { /* A is symmetric */
213:       testMatMatSolveTranspose = PETSC_TRUE;
214:       testMatSolveTranspose    = PETSC_TRUE;
215:     } else { /* superlu_dist does not support solving A^t x = rhs yet */
216:       testMatMatSolveTranspose = PETSC_FALSE;
217:       testMatSolveTranspose    = PETSC_FALSE;
218:     }
219:     ipack = 1;
220:     goto skipoptions;
221:   }
222: #endif
223: #if defined(PETSC_HAVE_MUMPS)
224:   PetscCall(PetscStrcmp(MATSOLVERMUMPS, pack, &match));
225:   if (match) {
226:     if (chol) {
227:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MUMPS CHOLESKY:\n"));
228:       PetscCall(MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_CHOLESKY, &F));
229:     } else {
230:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MUMPS LU:\n"));
231:       PetscCall(MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_LU, &F));
232:     }
233:     matsolvexx = PETSC_TRUE;
234:     if (test_mumps_opts) {
235:       /* test mumps options */
236:       PetscInt  icntl;
237:       PetscReal cntl;

239:       icntl = 2; /* sequential matrix ordering */
240:       PetscCall(MatMumpsSetIcntl(F, 7, icntl));

242:       cntl = 1.e-6; /* threshold for row pivot detection */
243:       PetscCall(MatMumpsSetIcntl(F, 24, 1));
244:       PetscCall(MatMumpsSetCntl(F, 3, cntl));
245:     }
246:     ipack = 2;
247:     goto skipoptions;
248:   }
249: #endif
250: #if defined(PETSC_HAVE_MKL_PARDISO)
251:   PetscCall(PetscStrcmp(MATSOLVERMKL_PARDISO, pack, &match));
252:   if (match) {
253:     if (chol) {
254:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MKL_PARDISO CHOLESKY:\n"));
255:       PetscCall(MatGetFactor(A, MATSOLVERMKL_PARDISO, MAT_FACTOR_CHOLESKY, &F));
256:     } else {
257:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MKL_PARDISO LU:\n"));
258:       PetscCall(MatGetFactor(A, MATSOLVERMKL_PARDISO, MAT_FACTOR_LU, &F));
259:     }
260:     ipack = 3;
261:     goto skipoptions;
262:   }
263: #endif
264: #if defined(PETSC_HAVE_CUDA)
265:   PetscCall(PetscStrcmp(MATSOLVERCUSPARSE, pack, &match));
266:   if (match) {
267:     if (chol) {
268:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, " CUSPARSE CHOLESKY:\n"));
269:       PetscCall(MatGetFactor(A, MATSOLVERCUSPARSE, MAT_FACTOR_CHOLESKY, &F));
270:     } else {
271:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, " CUSPARSE LU:\n"));
272:       PetscCall(MatGetFactor(A, MATSOLVERCUSPARSE, MAT_FACTOR_LU, &F));
273:     }
274:     testMatSolveTranspose    = PETSC_FALSE;
275:     testMatMatSolveTranspose = PETSC_FALSE;
276:     ipack                    = 4;
277:     goto skipoptions;
278:   }
279: #endif
280:   /* PETSc */
281:   match = PETSC_TRUE;
282:   if (match) {
283:     if (chol) {
284:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, " PETSC CHOLESKY:\n"));
285:       PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &F));
286:     } else {
287:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, " PETSC LU:\n"));
288:       PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &F));
289:     }
290:     matsolvexx = PETSC_TRUE;
291:     ipack      = 5;
292:     goto skipoptions;
293:   }

295: skipoptions:
296:   PetscCall(MatFactorInfoInitialize(&info));
297:   info.fill      = 5.0;
298:   info.shifttype = (PetscReal)MAT_SHIFT_NONE;
299:   if (chol) {
300:     PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info));
301:   } else {
302:     PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, &info));
303:   }

305:   for (nfact = 0; nfact < 2; nfact++) {
306:     if (chol) {
307:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the CHOLESKY numfactorization \n", nfact));
308:       PetscCall(MatCholeskyFactorNumeric(F, A, &info));
309:     } else {
310:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the LU numfactorization \n", nfact));
311:       PetscCall(MatLUFactorNumeric(F, A, &info));
312:     }
313:     if (view) {
314:       PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO));
315:       PetscCall(MatView(F, PETSC_VIEWER_STDOUT_WORLD));
316:       PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
317:       view = PETSC_FALSE;
318:     }

320: #if defined(PETSC_HAVE_SUPERLU_DIST)
321:     if (ipack == 1) { /* Test MatSuperluDistGetDiagU()
322:        -- input: matrix factor F; output: main diagonal of matrix U on all processes */
323:       PetscInt     M;
324:       PetscScalar *diag;
325:   #if !defined(PETSC_USE_COMPLEX)
326:       PetscInt nneg, nzero, npos;
327:   #endif

329:       PetscCall(MatGetSize(F, &M, NULL));
330:       PetscCall(PetscMalloc1(M, &diag));
331:       PetscCall(MatSuperluDistGetDiagU(F, diag));
332:       PetscCall(PetscFree(diag));

334:   #if !defined(PETSC_USE_COMPLEX)
335:       /* Test MatGetInertia() */
336:       if (symm) { /* A is symmetric */
337:         PetscCall(MatGetInertia(F, &nneg, &nzero, &npos));
338:         PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, " MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n", nneg, nzero, npos));
339:       }
340:   #endif
341:     }
342: #endif

344: #if defined(PETSC_HAVE_MUMPS)
345:     /* mumps interface allows repeated call of MatCholeskyFactorSymbolic(), while the succession calls do nothing */
346:     if (ipack == 2) {
347:       if (chol) {
348:         PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info));
349:         PetscCall(MatCholeskyFactorNumeric(F, A, &info));
350:       } else {
351:         PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, &info));
352:         PetscCall(MatLUFactorNumeric(F, A, &info));
353:       }
354:     }
355: #endif

357:     /* Test MatMatSolve(), A X = B, where B can be dense or sparse */
358:     if (testMatMatSolve) {
359:       if (!nfact) {
360:         PetscCall(MatMatMult(Ae, C, MAT_INITIAL_MATRIX, 2.0, &RHS));
361:       } else {
362:         PetscCall(MatMatMult(Ae, C, MAT_REUSE_MATRIX, 2.0, &RHS));
363:       }
364:       for (nsolve = 0; nsolve < 2; nsolve++) {
365:         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "   %" PetscInt_FMT "-the MatMatSolve \n", nsolve));
366:         PetscCall(MatMatSolve(F, RHS, X));

368:         /* Check the error */
369:         PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
370:         PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
371:         if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the MatMatSolve: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve));
372:       }

374:       if (matsolvexx) {
375:         /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix */
376:         PetscCall(MatCopy(RHS, X, SAME_NONZERO_PATTERN));
377:         PetscCall(MatMatSolve(F, X, X));
378:         /* Check the error */
379:         PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
380:         PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
381:         if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMatSolve(F,RHS,RHS): Norm of error %g\n", (double)norm));
382:       }

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

387:         PetscCall(MatTranspose(RHS, MAT_INITIAL_MATRIX, &RHST));
388:         PetscCall(MatConvert(RHST, MATAIJ, MAT_INITIAL_MATRIX, &spRHST));
389:         PetscCall(MatCreateTranspose(spRHST, &spRHS));
390:         for (nsolve = 0; nsolve < 2; nsolve++) {
391:           PetscCall(PetscPrintf(PETSC_COMM_WORLD, "   %" PetscInt_FMT "-the sparse MatMatSolve \n", nsolve));
392:           PetscCall(MatMatSolve(F, spRHS, X));

394:           /* Check the error */
395:           PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
396:           PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
397:           if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the sparse MatMatSolve: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve));
398:         }
399:         PetscCall(MatDestroy(&spRHST));
400:         PetscCall(MatDestroy(&spRHS));
401:         PetscCall(MatDestroy(&RHST));
402:       }
403:     }

405:     /* Test testMatMatSolveTranspose(), A^T X = B, where B can be dense or sparse */
406:     if (testMatMatSolveTranspose) {
407:       if (!nfact) {
408:         PetscCall(MatTransposeMatMult(Ae, C, MAT_INITIAL_MATRIX, 2.0, &RHS1));
409:       } else {
410:         PetscCall(MatTransposeMatMult(Ae, C, MAT_REUSE_MATRIX, 2.0, &RHS1));
411:       }

413:       for (nsolve = 0; nsolve < 2; nsolve++) {
414:         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "   %" PetscInt_FMT "-the MatMatSolveTranspose\n", nsolve));
415:         PetscCall(MatMatSolveTranspose(F, RHS1, X));

417:         /* Check the error */
418:         PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
419:         PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
420:         if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the MatMatSolveTranspose: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve));
421:       }

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

426:         PetscCall(MatTranspose(RHS1, MAT_INITIAL_MATRIX, &RHST));
427:         PetscCall(MatConvert(RHST, MATAIJ, MAT_INITIAL_MATRIX, &spRHST));
428:         PetscCall(MatCreateTranspose(spRHST, &spRHS));
429:         for (nsolve = 0; nsolve < 2; nsolve++) {
430:           PetscCall(MatMatSolveTranspose(F, spRHS, X));

432:           /* Check the error */
433:           PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
434:           PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
435:           if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the sparse MatMatSolveTranspose: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve));
436:         }
437:         PetscCall(MatDestroy(&spRHST));
438:         PetscCall(MatDestroy(&spRHS));
439:         PetscCall(MatDestroy(&RHST));
440:       }
441:     }

443:     /* Test MatSolve() */
444:     if (testMatSolve) {
445:       for (nsolve = 0; nsolve < 2; nsolve++) {
446:         PetscCall(VecSetRandom(x, rand));
447:         PetscCall(VecCopy(x, u));
448:         PetscCall(MatMult(Ae, x, b));

450:         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "   %" PetscInt_FMT "-the MatSolve \n", nsolve));
451:         PetscCall(MatSolve(F, b, x));

453:         /* Check the error */
454:         PetscCall(VecAXPY(u, -1.0, x)); /* u <- (-1.0)x + u */
455:         PetscCall(VecNorm(u, NORM_2, &norm));
456:         if (norm > tol) {
457:           PetscReal resi;
458:           PetscCall(MatMult(Ae, x, u));   /* u = A*x */
459:           PetscCall(VecAXPY(u, -1.0, b)); /* u <- (-1.0)b + u */
460:           PetscCall(VecNorm(u, NORM_2, &resi));
461:           PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolve: Norm of error %g, resi %g, numfact %" PetscInt_FMT "\n", (double)norm, (double)resi, nfact));
462:         }
463:       }
464:     }

466:     /* Test MatSolveTranspose() */
467:     if (testMatSolveTranspose) {
468:       for (nsolve = 0; nsolve < 2; nsolve++) {
469:         PetscCall(VecSetRandom(x, rand));
470:         PetscCall(VecCopy(x, u));
471:         PetscCall(MatMultTranspose(Ae, x, b));

473:         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "   %" PetscInt_FMT "-the MatSolveTranspose\n", nsolve));
474:         PetscCall(MatSolveTranspose(F, b, x));

476:         /* Check the error */
477:         PetscCall(VecAXPY(u, -1.0, x)); /* u <- (-1.0)x + u */
478:         PetscCall(VecNorm(u, NORM_2, &norm));
479:         if (norm > tol) {
480:           PetscReal resi;
481:           PetscCall(MatMultTranspose(Ae, x, u)); /* u = A*x */
482:           PetscCall(VecAXPY(u, -1.0, b));        /* u <- (-1.0)b + u */
483:           PetscCall(VecNorm(u, NORM_2, &resi));
484:           PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolveTranspose: Norm of error %g, resi %g, numfact %" PetscInt_FMT "\n", (double)norm, (double)resi, nfact));
485:         }
486:       }
487:     }
488:   }

490:   /* Free data structures */
491:   PetscCall(MatDestroy(&Ae));
492:   PetscCall(MatDestroy(&A));
493:   PetscCall(MatDestroy(&C));
494:   PetscCall(MatDestroy(&F));
495:   PetscCall(MatDestroy(&X));
496:   PetscCall(MatDestroy(&RHS));
497:   PetscCall(MatDestroy(&RHS1));

499:   PetscCall(PetscRandomDestroy(&rand));
500:   PetscCall(ISDestroy(&perm));
501:   PetscCall(ISDestroy(&iperm));
502:   PetscCall(VecDestroy(&x));
503:   PetscCall(VecDestroy(&b));
504:   PetscCall(VecDestroy(&u));
505:   PetscCall(PetscFinalize());
506:   return 0;
507: }

509: /*TEST

511:    test:
512:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
513:       args: -f ${DATAFILESPATH}/matrices/medium -mat_solver_type petsc
514:       output_file: output/ex125.out

516:    test:
517:       suffix: 2
518:       args: -mat_solver_type petsc
519:       output_file: output/ex125.out

521:    test:
522:       suffix: mkl_pardiso
523:       requires: mkl_pardiso datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
524:       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mkl_pardiso

526:    test:
527:       suffix: mkl_pardiso_2
528:       requires: mkl_pardiso
529:       args: -mat_solver_type mkl_pardiso
530:       output_file: output/ex125_mkl_pardiso.out

532:    test:
533:       suffix: mumps
534:       requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
535:       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps
536:       output_file: output/ex125_mumps_seq.out

538:    test:
539:       suffix: mumps_nest
540:       requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
541:       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}}
542:       output_file: output/ex125_mumps_seq.out

544:    test:
545:       suffix: mumps_2
546:       nsize: 3
547:       requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
548:       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps
549:       output_file: output/ex125_mumps_par.out

551:    test:
552:       suffix: mumps_2_nest
553:       nsize: 3
554:       requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
555:       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}}
556:       output_file: output/ex125_mumps_par.out

558:    test:
559:       suffix: mumps_3
560:       requires: mumps
561:       args: -mat_solver_type mumps
562:       output_file: output/ex125_mumps_seq.out

564:    test:
565:       suffix: mumps_3_nest
566:       requires: mumps
567:       args: -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}}
568:       output_file: output/ex125_mumps_seq.out

570:    test:
571:       suffix: mumps_4
572:       nsize: 3
573:       requires: mumps
574:       args: -mat_solver_type mumps
575:       output_file: output/ex125_mumps_par.out

577:    test:
578:       suffix: mumps_4_nest
579:       nsize: 3
580:       requires: mumps
581:       args: -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}}
582:       output_file: output/ex125_mumps_par.out

584:    test:
585:       suffix: mumps_5
586:       nsize: 3
587:       requires: mumps
588:       args: -mat_solver_type mumps -cholesky
589:       output_file: output/ex125_mumps_par_cholesky.out

591:    test:
592:       suffix: mumps_5_nest
593:       nsize: 3
594:       requires: mumps
595:       args: -mat_solver_type mumps -cholesky -test_nest -test_nest_bordered {{0 1}}
596:       output_file: output/ex125_mumps_par_cholesky.out

598:    test:
599:       suffix: superlu
600:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu
601:       args: -f ${DATAFILESPATH}/matrices/medium -mat_solver_type superlu
602:       output_file: output/ex125_superlu.out

604:    test:
605:       suffix: superlu_dist
606:       nsize: {{1 3}}
607:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu_dist
608:       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
609:       output_file: output/ex125_superlu_dist.out

611:    test:
612:       suffix: superlu_dist_2
613:       nsize: {{1 3}}
614:       requires: superlu_dist !complex
615:       args: -n 36 -mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
616:       output_file: output/ex125_superlu_dist.out

618:    test:
619:       suffix: superlu_dist_3
620:       nsize: {{1 3}}
621:       requires: superlu_dist !complex
622:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu_dist
623:       args: -f ${DATAFILESPATH}/matrices/medium -mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
624:       output_file: output/ex125_superlu_dist_nonsymmetric.out

626:    test:
627:       suffix: superlu_dist_complex
628:       nsize: 3
629:       requires: datafilespath double superlu_dist complex !defined(PETSC_USE_64BIT_INDICES)
630:       args: -f ${DATAFILESPATH}/matrices/farzad_B_rhs -mat_solver_type superlu_dist
631:       output_file: output/ex125_superlu_dist_complex.out

633:    test:
634:       suffix: superlu_dist_complex_2
635:       nsize: 3
636:       requires: superlu_dist complex
637:       args: -mat_solver_type superlu_dist
638:       output_file: output/ex125_superlu_dist_complex_2.out

640:    test:
641:       suffix: cusparse
642:       requires: cuda datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
643:       #TODO: fix the bug with cholesky
644:       #args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type cusparse -cholesky {{0 1}separate output}
645:       args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type cusparse -cholesky {{0}separate output}

647:    test:
648:       suffix: cusparse_2
649:       requires: cuda
650:       args: -mat_type aijcusparse -mat_solver_type cusparse -cholesky {{0 1}separate output}

652: TEST*/