Actual source code: ex52.c

      1: static char help[] = "Solves a linear system in parallel with KSP. Modified from ex2.c \n\
      2:                       Illustrate how to use external packages MUMPS, SUPERLU and STRUMPACK \n\
      3: Input parameters include:\n\
      4:   -random_exact_sol : use a random exact solution vector\n\
      5:   -view_exact_sol   : write exact solution vector to stdout\n\
      6:   -m <mesh_x>       : number of mesh points in x-direction\n\
      7:   -n <mesh_y>       : number of mesh points in y-direction\n\n";

      9: #include <petscksp.h>

     11: #if defined(PETSC_HAVE_MUMPS)
     12: /* Subroutine contributed by Varun Hiremath */
     13: PetscErrorCode printMumpsMemoryInfo(Mat F)
     14: {
     15:   PetscInt maxMem, sumMem;

     17:   PetscFunctionBeginUser;
     18:   PetscCall(MatMumpsGetInfog(F, 16, &maxMem));
     19:   PetscCall(MatMumpsGetInfog(F, 17, &sumMem));
     20:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n MUMPS INFOG(16) :: Max memory in MB = %" PetscInt_FMT, maxMem));
     21:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n MUMPS INFOG(17) :: Sum memory in MB = %" PetscInt_FMT "\n", sumMem));
     22:   PetscFunctionReturn(PETSC_SUCCESS);
     23: }
     24: #endif

     26: int main(int argc, char **args)
     27: {
     28:   Vec         x, b, u; /* approx solution, RHS, exact solution */
     29:   Mat         A, F;
     30:   KSP         ksp; /* linear solver context */
     31:   PC          pc;
     32:   PetscRandom rctx; /* random number generator context */
     33:   PetscReal   norm; /* norm of solution error */
     34:   PetscInt    i, j, Ii, J, Istart, Iend, m = 8, n = 7, its;
     35:   PetscBool   flg = PETSC_FALSE, flg_ilu = PETSC_FALSE, flg_ch = PETSC_FALSE;
     36: #if defined(PETSC_HAVE_MUMPS)
     37:   char      tmpdir[PETSC_MAX_PATH_LEN];
     38:   PetscBool flg_mumps = PETSC_FALSE, flg_mumps_ch = PETSC_FALSE, test_mumps_ooc_api = PETSC_FALSE;
     39: #endif
     40: #if defined(PETSC_HAVE_SUPERLU) || defined(PETSC_HAVE_SUPERLU_DIST)
     41:   PetscBool flg_superlu = PETSC_FALSE;
     42: #endif
     43: #if defined(PETSC_HAVE_STRUMPACK)
     44:   PetscBool flg_strumpack = PETSC_FALSE;
     45: #endif
     46:   PetscScalar   v;
     47:   PetscMPIInt   rank, size;
     48:   PetscLogStage stage;

     50: #if defined(PETSC_HAVE_STRUMPACK) && defined(PETSC_HAVE_SLATE)
     51:   PETSC_MPI_THREAD_REQUIRED = MPI_THREAD_MULTIPLE;
     52: #endif
     53:   PetscFunctionBeginUser;
     54:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
     55:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
     56:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
     57:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
     58:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
     59:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     60:          Compute the matrix and right-hand-side vector that define
     61:          the linear system, Ax = b.
     62:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     63:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
     64:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n));
     65:   PetscCall(MatSetFromOptions(A));
     66:   PetscCall(MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL));
     67:   PetscCall(MatSeqAIJSetPreallocation(A, 5, NULL));
     68:   PetscCall(MatSetUp(A));

     70:   /*
     71:      Currently, all PETSc parallel matrix formats are partitioned by
     72:      contiguous chunks of rows across the processors.  Determine which
     73:      rows of the matrix are locally owned.
     74:   */
     75:   PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));

     77:   /*
     78:      Set matrix elements for the 2-D, five-point stencil in parallel.
     79:       - Each processor needs to insert only elements that it owns
     80:         locally (but any non-local elements will be sent to the
     81:         appropriate processor during matrix assembly).
     82:       - Always specify global rows and columns of matrix entries.

     84:      Note: this uses the less common natural ordering that orders first
     85:      all the unknowns for x = h then for x = 2h etc; Hence you see J = Ii +- n
     86:      instead of J = I +- m as you might expect. The more standard ordering
     87:      would first do all variables for y = h, then y = 2h etc.

     89:    */
     90:   PetscCall(PetscLogStageRegister("Assembly", &stage));
     91:   PetscCall(PetscLogStagePush(stage));
     92:   for (Ii = Istart; Ii < Iend; Ii++) {
     93:     v = -1.0;
     94:     i = Ii / n;
     95:     j = Ii - i * n;
     96:     if (i > 0) {
     97:       J = Ii - n;
     98:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
     99:     }
    100:     if (i < m - 1) {
    101:       J = Ii + n;
    102:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
    103:     }
    104:     if (j > 0) {
    105:       J = Ii - 1;
    106:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
    107:     }
    108:     if (j < n - 1) {
    109:       J = Ii + 1;
    110:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
    111:     }
    112:     v = 4.0;
    113:     PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
    114:   }

    116:   /*
    117:      Assemble matrix, using the 2-step process:
    118:        MatAssemblyBegin(), MatAssemblyEnd()
    119:      Computations can be done while messages are in transition
    120:      by placing code between these two statements.
    121:   */
    122:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    123:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    124:   PetscCall(PetscLogStagePop());

    126:   /* A is symmetric. Set symmetric flag to enable ICC/Cholesky preconditioner */
    127:   PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));

    129:   /* Create parallel vectors */
    130:   PetscCall(MatCreateVecs(A, &u, &b));
    131:   PetscCall(VecDuplicate(u, &x));

    133:   /*
    134:      Set exact solution; then compute right-hand-side vector.
    135:      By default we use an exact solution of a vector with all
    136:      elements of 1.0;  Alternatively, using the runtime option
    137:      -random_sol forms a solution vector with random components.
    138:   */
    139:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-random_exact_sol", &flg, NULL));
    140:   if (flg) {
    141:     PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
    142:     PetscCall(PetscRandomSetFromOptions(rctx));
    143:     PetscCall(VecSetRandom(u, rctx));
    144:     PetscCall(PetscRandomDestroy(&rctx));
    145:   } else {
    146:     PetscCall(VecSet(u, 1.0));
    147:   }
    148:   PetscCall(MatMult(A, u, b));

    150:   /*
    151:      View the exact solution vector if desired
    152:   */
    153:   flg = PETSC_FALSE;
    154:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_exact_sol", &flg, NULL));
    155:   if (flg) PetscCall(VecView(u, PETSC_VIEWER_STDOUT_WORLD));

    157:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    158:                 Create the linear solver and set various options
    159:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    161:   /*
    162:      Create linear solver context
    163:   */
    164:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
    165:   PetscCall(KSPSetOperators(ksp, A, A));

    167:   /*
    168:     Example of how to use external package MUMPS
    169:     Note: runtime options
    170:           '-ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps -mat_mumps_icntl_7 3 -mat_mumps_icntl_1 0.0'
    171:           are equivalent to these procedural calls
    172:   */
    173: #if defined(PETSC_HAVE_MUMPS)
    174:   flg_mumps    = PETSC_FALSE;
    175:   flg_mumps_ch = PETSC_FALSE;
    176:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_mumps_lu", &flg_mumps, NULL));
    177:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_mumps_ch", &flg_mumps_ch, NULL));
    178:   if (flg_mumps || flg_mumps_ch) {
    179:     PetscCall(KSPSetType(ksp, KSPPREONLY));
    180:     PetscInt  ival, icntl;
    181:     PetscReal val;
    182:     PetscCall(KSPGetPC(ksp, &pc));
    183:     if (flg_mumps) PetscCall(PCSetType(pc, PCLU));
    184:     else if (flg_mumps_ch) {
    185:       PetscCall(MatSetOption(A, MAT_SPD, PETSC_TRUE)); /* set MUMPS id%SYM=1 */
    186:       PetscCall(PCSetType(pc, PCCHOLESKY));
    187:     }
    188:     PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS));
    189:     PetscCall(PCFactorSetUpMatSolverType(pc)); /* call MatGetFactor() to create F */
    190:     PetscCall(PCFactorGetMatrix(pc, &F));
    191:     PetscCall(MatMumpsSetIcntl(F, 24, 1));
    192:     PetscCall(MatMumpsGetIcntl(F, 24, &ival));
    193:     PetscCheck(ival == 1, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "ICNTL(24) = %" PetscInt_FMT " (!= 1)", ival);
    194:     PetscCall(MatMumpsSetCntl(F, 3, 1e-6));
    195:     PetscCall(MatMumpsGetCntl(F, 3, &val));
    196:     PetscCheck(PetscEqualReal(val, 1e-6), PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "CNTL(3) = %g (!= %g)", (double)val, 1e-6);
    197:     if (flg_mumps) {
    198:       /* Zero the first and last rows in the rank, they should then show up in corresponding null pivot rows output via
    199:          MatMumpsGetNullPivots */
    200:       flg = PETSC_FALSE;
    201:       PetscCall(PetscOptionsGetBool(NULL, NULL, "-zero_first_and_last_rows", &flg, NULL));
    202:       if (flg) {
    203:         PetscInt rows[2];
    204:         rows[0] = Istart;   /* first row of the rank */
    205:         rows[1] = Iend - 1; /* last row of the rank */
    206:         PetscCall(MatZeroRows(A, 2, rows, 0.0, NULL, NULL));
    207:       }
    208:       /* Get memory estimates from MUMPS' MatLUFactorSymbolic(), e.g. INFOG(16), INFOG(17).
    209:          KSPSetUp() below will do nothing inside MatLUFactorSymbolic() */
    210:       MatFactorInfo info;
    211:       PetscCall(MatFactorInfoInitialize(&info));
    212:       PetscCall(MatLUFactorSymbolic(F, A, NULL, NULL, &info));
    213:       flg = PETSC_FALSE;
    214:       PetscCall(PetscOptionsGetBool(NULL, NULL, "-print_mumps_memory", &flg, NULL));
    215:       if (flg) PetscCall(printMumpsMemoryInfo(F));
    216:     }

    218:     /* sequential ordering */
    219:     icntl = 7;
    220:     ival  = 2;
    221:     PetscCall(MatMumpsSetIcntl(F, icntl, ival));

    223:     /* threshold for row pivot detection */
    224:     PetscCall(MatMumpsGetIcntl(F, 24, &ival));
    225:     PetscCheck(ival == 1, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "ICNTL(24) = %" PetscInt_FMT " (!= 1)", ival);
    226:     icntl = 3;
    227:     PetscCall(MatMumpsGetCntl(F, icntl, &val));
    228:     PetscCheck(PetscEqualReal(val, 1e-6), PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "CNTL(3) = %g (!= %g)", (double)val, 1e-6);

    230:     /* compute determinant of A */
    231:     PetscCall(MatMumpsSetIcntl(F, 33, 1));
    232:   }
    233: #endif

    235:   /*
    236:     Example of how to use external package SuperLU
    237:     Note: runtime options
    238:           '-ksp_type preonly -pc_type ilu -pc_factor_mat_solver_type superlu -mat_superlu_ilu_droptol 1.e-8'
    239:           are equivalent to these procedual calls
    240:   */
    241: #if defined(PETSC_HAVE_SUPERLU) || defined(PETSC_HAVE_SUPERLU_DIST)
    242:   flg_ilu     = PETSC_FALSE;
    243:   flg_superlu = PETSC_FALSE;
    244:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_superlu_lu", &flg_superlu, NULL));
    245:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_superlu_ilu", &flg_ilu, NULL));
    246:   if (flg_superlu || flg_ilu) {
    247:     PetscCall(KSPSetType(ksp, KSPPREONLY));
    248:     PetscCall(KSPGetPC(ksp, &pc));
    249:     if (flg_superlu) PetscCall(PCSetType(pc, PCLU));
    250:     else if (flg_ilu) PetscCall(PCSetType(pc, PCILU));
    251:     if (size == 1) {
    252:   #if !defined(PETSC_HAVE_SUPERLU)
    253:       SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This test requires SUPERLU");
    254:   #else
    255:       PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERSUPERLU));
    256:   #endif
    257:     } else {
    258:   #if !defined(PETSC_HAVE_SUPERLU_DIST)
    259:       SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This test requires SUPERLU_DIST");
    260:   #else
    261:       PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERSUPERLU_DIST));
    262:   #endif
    263:     }
    264:     PetscCall(PCFactorSetUpMatSolverType(pc)); /* call MatGetFactor() to create F */
    265:     PetscCall(PCFactorGetMatrix(pc, &F));
    266:   #if defined(PETSC_HAVE_SUPERLU)
    267:     if (size == 1) PetscCall(MatSuperluSetILUDropTol(F, 1.e-8));
    268:   #endif
    269:   }
    270: #endif

    272:   /*
    273:     Example of how to use external package STRUMPACK
    274:     Note: runtime options
    275:           '-pc_type lu/ilu \
    276:            -pc_factor_mat_solver_type strumpack \
    277:            -mat_strumpack_reordering GEOMETRIC \
    278:            -mat_strumpack_geometric_xyz n,m \
    279:            -mat_strumpack_colperm 0 \
    280:            -mat_strumpack_compression_rel_tol 1.e-3 \
    281:            -mat_strumpack_compression_min_sep_size 15 \
    282:            -mat_strumpack_leaf_size 4'
    283:        are equivalent to these procedural calls

    285:     We refer to the STRUMPACK manual for more info on
    286:     how to tune the preconditioner, see for instance:
    287:      https://portal.nersc.gov/project/sparse/strumpack/master/prec.html
    288:   */
    289: #if defined(PETSC_HAVE_STRUMPACK)
    290:   flg_ilu       = PETSC_FALSE;
    291:   flg_strumpack = PETSC_FALSE;
    292:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_strumpack_lu", &flg_strumpack, NULL));
    293:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_strumpack_ilu", &flg_ilu, NULL));
    294:   if (flg_strumpack || flg_ilu) {
    295:     PetscCall(KSPSetType(ksp, KSPPREONLY));
    296:     PetscCall(KSPGetPC(ksp, &pc));
    297:     if (flg_strumpack) PetscCall(PCSetType(pc, PCLU));
    298:     else if (flg_ilu) PetscCall(PCSetType(pc, PCILU));
    299:   #if !defined(PETSC_HAVE_STRUMPACK)
    300:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This test requires STRUMPACK");
    301:   #endif
    302:     PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERSTRUMPACK));
    303:     PetscCall(PCFactorSetUpMatSolverType(pc)); /* call MatGetFactor() to create F */
    304:     PetscCall(PCFactorGetMatrix(pc, &F));

    306:     /* Set the fill-reducing reordering, MAT_STRUMPACK_METIS is       */
    307:     /* always supported, but is sequential. Parallel alternatives are */
    308:     /* MAT_STRUMPACK_PARMETIS and MAT_STRUMPACK_PTSCOTCH, but         */
    309:     /* strumpack needs to be configured with support for these.       */
    310:     /*PetscCall(MatSTRUMPACKSetReordering(F, MAT_STRUMPACK_METIS));   */
    311:     /* However, since this is a problem on a regular grid, we can use */
    312:     /* a simple geometric nested dissection implementation, which     */
    313:     /* requires passing the grid dimensions to strumpack.             */
    314:     PetscCall(MatSTRUMPACKSetReordering(F, MAT_STRUMPACK_GEOMETRIC));
    315:     PetscCall(MatSTRUMPACKSetGeometricNxyz(F, n, m, PETSC_DECIDE));
    316:     /* These are optional, defaults are 1                             */
    317:     PetscCall(MatSTRUMPACKSetGeometricComponents(F, 1));
    318:     PetscCall(MatSTRUMPACKSetGeometricWidth(F, 1));

    320:     /* Since this is a simple discretization, the diagonal is always  */
    321:     /* nonzero, and there is no need for the extra MC64 permutation.  */
    322:     PetscCall(MatSTRUMPACKSetColPerm(F, PETSC_FALSE));

    324:     if (flg_ilu) {
    325:       /* The compression tolerance used when doing low-rank compression */
    326:       /* in the preconditioner. This is problem specific!               */
    327:       PetscCall(MatSTRUMPACKSetCompRelTol(F, 1.e-3));

    329:       /* Set a small minimum (dense) matrix size for compression to     */
    330:       /* demonstrate the preconditioner on small problems.              */
    331:       /* For performance the default value should be better.            */
    332:       /* This size corresponds to the size of separators in the graph.  */
    333:       /* For instance on an m x n mesh, the top level separator is of   */
    334:       /* size m (if m <= n)                                             */
    335:       /*PetscCall(MatSTRUMPACKSetCompMinSepSize(F,15));*/

    337:       /* Set the size of the diagonal blocks (the leafs) in the HSS     */
    338:       /* approximation. The default value should be better for real     */
    339:       /* problems. This is mostly for illustration on a small problem.  */
    340:       /*PetscCall(MatSTRUMPACKSetCompLeafSize(F,4));*/
    341:     }
    342:   }
    343: #endif

    345:   /*
    346:     Example of how to use procedural calls that are equivalent to
    347:           '-ksp_type preonly -pc_type lu/ilu -pc_factor_mat_solver_type petsc'
    348:   */
    349:   flg     = PETSC_FALSE;
    350:   flg_ilu = PETSC_FALSE;
    351:   flg_ch  = PETSC_FALSE;
    352:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_petsc_lu", &flg, NULL));
    353:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_petsc_ilu", &flg_ilu, NULL));
    354:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_petsc_ch", &flg_ch, NULL));
    355:   if (flg || flg_ilu || flg_ch) {
    356:     Vec diag;

    358:     PetscCall(KSPSetType(ksp, KSPPREONLY));
    359:     PetscCall(KSPGetPC(ksp, &pc));
    360:     if (flg) PetscCall(PCSetType(pc, PCLU));
    361:     else if (flg_ilu) PetscCall(PCSetType(pc, PCILU));
    362:     else if (flg_ch) PetscCall(PCSetType(pc, PCCHOLESKY));
    363:     PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERPETSC));
    364:     PetscCall(PCFactorSetUpMatSolverType(pc)); /* call MatGetFactor() to create F */
    365:     PetscCall(PCFactorGetMatrix(pc, &F));

    367:     /* Test MatGetDiagonal() */
    368:     PetscCall(KSPSetUp(ksp));
    369:     PetscCall(VecDuplicate(x, &diag));
    370:     PetscCall(MatGetDiagonal(F, diag));
    371:     /* PetscCall(VecView(diag,PETSC_VIEWER_STDOUT_WORLD)); */
    372:     PetscCall(VecDestroy(&diag));
    373:   }

    375:   PetscCall(KSPSetFromOptions(ksp));

    377: #if defined(PETSC_HAVE_MUMPS)
    378:   // The OOC options must be set before PCSetUp()/KSPSetUp(), as MUMPS requires them be set after the initialization phase and before the (numeric) factorization phase.
    379:   if (flg_mumps || flg_mumps_ch) {
    380:     PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_mumps_ooc_api", &test_mumps_ooc_api, NULL));
    381:     if (test_mumps_ooc_api) {
    382:       PetscCall(PetscStrncpy(tmpdir, "OOC_XXXXXX", sizeof(tmpdir)));
    383:       if (rank == 0) PetscCall(PetscMkdtemp(tmpdir));
    384:       PetscCallMPI(MPI_Bcast(tmpdir, 11, MPI_CHAR, 0, PETSC_COMM_WORLD));
    385:       PetscCall(MatMumpsSetIcntl(F, 22, 1)); // enable out-of-core
    386:       PetscCall(MatMumpsSetOocTmpDir(F, tmpdir));
    387:     }
    388:   }
    389: #endif

    391:   /* Get info from matrix factors */
    392:   PetscCall(KSPSetUp(ksp));

    394: #if defined(PETSC_HAVE_MUMPS)
    395:   if (flg_mumps || flg_mumps_ch) {
    396:     PetscInt  icntl, infog34, num_null_pivots, *null_pivots;
    397:     PetscReal cntl, rinfo12, rinfo13;
    398:     icntl = 3;
    399:     PetscCall(MatMumpsGetCntl(F, icntl, &cntl));

    401:     /* compute determinant and check for any null pivots*/
    402:     if (rank == 0) {
    403:       PetscCall(MatMumpsGetInfog(F, 34, &infog34));
    404:       PetscCall(MatMumpsGetRinfog(F, 12, &rinfo12));
    405:       PetscCall(MatMumpsGetRinfog(F, 13, &rinfo13));
    406:       PetscCall(MatMumpsGetNullPivots(F, &num_null_pivots, &null_pivots));
    407:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Mumps row pivot threshold = %g\n", (double)cntl));
    408:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Mumps determinant = (%g, %g) * 2^%" PetscInt_FMT " \n", (double)rinfo12, (double)rinfo13, infog34));
    409:       if (num_null_pivots > 0) {
    410:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Mumps num of null pivots detected = %" PetscInt_FMT "\n", num_null_pivots));
    411:         PetscCall(PetscSortInt(num_null_pivots, null_pivots)); /* just make the printf deterministic */
    412:         for (j = 0; j < num_null_pivots; j++) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Mumps row with null pivots is = %" PetscInt_FMT "\n", null_pivots[j]));
    413:       }
    414:       PetscCall(PetscFree(null_pivots));
    415:     }
    416:   }
    417: #endif

    419:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    420:                       Solve the linear system
    421:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    422:   PetscCall(KSPSolve(ksp, b, x));

    424:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    425:                       Check solution and clean up
    426:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    427:   PetscCall(VecAXPY(x, -1.0, u));
    428:   PetscCall(VecNorm(x, NORM_2, &norm));
    429:   PetscCall(KSPGetIterationNumber(ksp, &its));

    431:   /*
    432:      Print convergence information.  PetscPrintf() produces a single
    433:      print statement from all processes that share a communicator.
    434:      An alternative is PetscFPrintf(), which prints to a file.
    435:   */
    436:   if (norm < 1.e-12) {
    437:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error < 1.e-12 iterations %" PetscInt_FMT "\n", its));
    438:   } else {
    439:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g iterations %" PetscInt_FMT "\n", (double)norm, its));
    440:   }

    442: #if defined(PETSC_HAVE_MUMPS)
    443:   // Get the OCC tmpdir via F before it is destroyed in KSPDestroy().
    444:   if (test_mumps_ooc_api) {
    445:     const char *dir;

    447:     PetscCall(MatMumpsGetOocTmpDir(F, &dir));
    448:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, " OOC_TMPDIR = %s\n", dir));
    449:   }
    450: #endif

    452:   /*
    453:      Free work space.  All PETSc objects should be destroyed when they
    454:      are no longer needed.
    455:   */
    456:   PetscCall(KSPDestroy(&ksp));

    458: #if defined(PETSC_HAVE_MUMPS)
    459:   // Files created by MUMPS under the OOC tmpdir were automatically deleted by MUMPS with JOB_END (-2)
    460:   // during KSP/PCDestroy(), but we need to remove the tmpdir created by us.
    461:   if (test_mumps_ooc_api && rank == 0) PetscCall(PetscRMTree(tmpdir));
    462: #endif

    464:   PetscCall(VecDestroy(&u));
    465:   PetscCall(VecDestroy(&x));
    466:   PetscCall(VecDestroy(&b));
    467:   PetscCall(MatDestroy(&A));

    469:   /*
    470:      Always call PetscFinalize() before exiting a program.  This routine
    471:        - finalizes the PETSc libraries as well as MPI
    472:        - provides summary and diagnostic information if certain runtime
    473:          options are chosen (e.g., -log_view).
    474:   */
    475:   PetscCall(PetscFinalize());
    476:   return 0;
    477: }

    479: /*TEST

    481:    test:
    482:       args: -use_petsc_lu
    483:       output_file: output/ex52_2.out

    485:    testset:
    486:       suffix: mumps
    487:       nsize: 3
    488:       requires: mumps
    489:       output_file: output/ex52_1.out
    490:       args: -use_mumps_lu

    492:       test:
    493:         suffix: mumps

    495:       test:
    496:         suffix: mumps_ooc
    497:         args: -mat_mumps_icntl_22 1 -mat_mumps_ooc_tmpdir /tmp

    499:       test:
    500:         suffix: mumps_ooc_api
    501:         args: -mat_mumps_icntl_22 1 -test_mumps_ooc_api
    502:         filter: grep -v "OOC_TMPDIR"

    504:    test:
    505:       suffix: mumps_2
    506:       nsize: 3
    507:       requires: mumps
    508:       args: -use_mumps_ch
    509:       output_file: output/ex52_1.out

    511:    test:
    512:       suffix: mumps_3
    513:       nsize: 3
    514:       requires: mumps
    515:       args: -use_mumps_ch -mat_type sbaij
    516:       output_file: output/ex52_1.out

    518:    test:
    519:       suffix: mumps_4
    520:       nsize: 3
    521:       requires: mumps !complex !single
    522:       args: -use_mumps_lu -m 50 -n 50 -print_mumps_memory
    523:       output_file: output/ex52_4.out

    525:    test:
    526:       suffix: mumps_5
    527:       nsize: 3
    528:       requires: mumps !complex !single
    529:       args: -use_mumps_lu -m 50 -n 50 -zero_first_and_last_rows
    530:       output_file: output/ex52_5.out

    532:    test:
    533:       suffix: mumps_omp_2
    534:       nsize: 4
    535:       requires: mumps hwloc openmp pthread defined(PETSC_HAVE_MPI_PROCESS_SHARED_MEMORY)
    536:       args: -use_mumps_lu -mat_mumps_use_omp_threads 2
    537:       output_file: output/ex52_1.out

    539:    test:
    540:       suffix: mumps_omp_3
    541:       nsize: 4
    542:       requires: mumps hwloc openmp pthread defined(PETSC_HAVE_MPI_PROCESS_SHARED_MEMORY)
    543:       args: -use_mumps_ch -mat_mumps_use_omp_threads 3 -mat_mumps_icntl_48 0
    544:       # Ignore the warning since we are intentionally testing the imbalanced case
    545:       filter: grep -v "Warning: number of OpenMP threads"
    546:       output_file: output/ex52_1.out

    548:    test:
    549:       suffix: mumps_omp_4
    550:       nsize: 4
    551:       requires: mumps hwloc openmp pthread defined(PETSC_HAVE_MPI_PROCESS_SHARED_MEMORY)
    552:       # let PETSc guess a proper number for threads
    553:       args: -use_mumps_ch -mat_type sbaij -mat_mumps_use_omp_threads
    554:       output_file: output/ex52_1.out

    556:    testset:
    557:       suffix: strumpack_2
    558:       nsize: {{1 2}}
    559:       requires: strumpack
    560:       args: -use_strumpack_lu
    561:       output_file: output/ex52_3.out

    563:       test:
    564:         suffix: aij
    565:         args: -mat_type aij

    567:       test:
    568:         requires: kokkos_kernels
    569:         suffix: kok
    570:         args: -mat_type aijkokkos

    572:       test:
    573:         requires: cuda
    574:         suffix: cuda
    575:         args: -mat_type aijcusparse

    577:       test:
    578:         requires: hip
    579:         suffix: hip
    580:         args: -mat_type aijhipsparse

    582:    test:
    583:       suffix: strumpack_ilu
    584:       nsize: {{1 2}}
    585:       requires: strumpack
    586:       args: -use_strumpack_ilu
    587:       output_file: output/ex52_3.out

    589:    testset:
    590:       suffix: superlu_dist
    591:       nsize: {{1 2}}
    592:       requires: superlu superlu_dist
    593:       args: -use_superlu_lu
    594:       output_file: output/ex52_2.out

    596:       test:
    597:         suffix: aij
    598:         args: -mat_type aij

    600:       test:
    601:         requires: kokkos_kernels
    602:         suffix: kok
    603:         args: -mat_type aijkokkos

    605:       test:
    606:         requires: cuda
    607:         suffix: cuda
    608:         args: -mat_type aijcusparse

    610:       test:
    611:         requires: hip
    612:         suffix: hip
    613:         args: -mat_type aijhipsparse

    615:    test:
    616:       suffix: superlu_ilu
    617:       requires: superlu superlu_dist
    618:       args: -use_superlu_ilu
    619:       output_file: output/ex52_2.out

    621: TEST*/