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";
8
9 #include <petscksp.h>
10
11 #if defined(PETSC_HAVE_MUMPS)
12 /* Subroutine contributed by Varun Hiremath */
printMumpsMemoryInfo(Mat F)13 PetscErrorCode printMumpsMemoryInfo(Mat F)
14 {
15 PetscInt maxMem, sumMem;
16
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
25
main(int argc,char ** args)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;
49
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));
69
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));
76
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.
83
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.
88
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 }
115
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());
125
126 /* A is symmetric. Set symmetric flag to enable ICC/Cholesky preconditioner */
127 PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
128
129 /* Create parallel vectors */
130 PetscCall(MatCreateVecs(A, &u, &b));
131 PetscCall(VecDuplicate(u, &x));
132
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));
149
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));
156
157 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158 Create the linear solver and set various options
159 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
160
161 /*
162 Create linear solver context
163 */
164 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
165 PetscCall(KSPSetOperators(ksp, A, A));
166
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) {
184 PetscCall(PCSetType(pc, PCLU));
185 } else if (flg_mumps_ch) {
186 PetscCall(MatSetOption(A, MAT_SPD, PETSC_TRUE)); /* set MUMPS id%SYM=1 */
187 PetscCall(PCSetType(pc, PCCHOLESKY));
188 }
189 PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS));
190 PetscCall(PCFactorSetUpMatSolverType(pc)); /* call MatGetFactor() to create F */
191 PetscCall(PCFactorGetMatrix(pc, &F));
192 PetscCall(MatMumpsSetIcntl(F, 24, 1));
193 PetscCall(MatMumpsGetIcntl(F, 24, &ival));
194 PetscCheck(ival == 1, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "ICNTL(24) = %" PetscInt_FMT " (!= 1)", ival);
195 PetscCall(MatMumpsSetCntl(F, 3, 1e-6));
196 PetscCall(MatMumpsGetCntl(F, 3, &val));
197 PetscCheck(PetscEqualReal(val, 1e-6), PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "CNTL(3) = %g (!= %g)", (double)val, 1e-6);
198 if (flg_mumps) {
199 /* Zero the first and last rows in the rank, they should then show up in corresponding null pivot rows output via
200 MatMumpsGetNullPivots */
201 flg = PETSC_FALSE;
202 PetscCall(PetscOptionsGetBool(NULL, NULL, "-zero_first_and_last_rows", &flg, NULL));
203 if (flg) {
204 PetscInt rows[2];
205 rows[0] = Istart; /* first row of the rank */
206 rows[1] = Iend - 1; /* last row of the rank */
207 PetscCall(MatZeroRows(A, 2, rows, 0.0, NULL, NULL));
208 }
209 /* Get memory estimates from MUMPS' MatLUFactorSymbolic(), e.g. INFOG(16), INFOG(17).
210 KSPSetUp() below will do nothing inside MatLUFactorSymbolic() */
211 MatFactorInfo 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 }
217
218 /* sequential ordering */
219 icntl = 7;
220 ival = 2;
221 PetscCall(MatMumpsSetIcntl(F, icntl, ival));
222
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);
229
230 /* compute determinant of A */
231 PetscCall(MatMumpsSetIcntl(F, 33, 1));
232 }
233 #endif
234
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
271
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
284
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));
305
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));
319
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));
323
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));
328
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));*/
336
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
344
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;
357
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));
366
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 }
374
375 PetscCall(KSPSetFromOptions(ksp));
376
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
390
391 /* Get info from matrix factors */
392 PetscCall(KSPSetUp(ksp));
393
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));
400
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
418
419 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
420 Solve the linear system
421 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
422 PetscCall(KSPSolve(ksp, b, x));
423
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));
430
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 }
441
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;
446
447 PetscCall(MatMumpsGetOocTmpDir(F, &dir));
448 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " OOC_TMPDIR = %s\n", dir));
449 }
450 #endif
451
452 /*
453 Free work space. All PETSc objects should be destroyed when they
454 are no longer needed.
455 */
456 PetscCall(KSPDestroy(&ksp));
457
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
463
464 PetscCall(VecDestroy(&u));
465 PetscCall(VecDestroy(&x));
466 PetscCall(VecDestroy(&b));
467 PetscCall(MatDestroy(&A));
468
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 }
478
479 /*TEST
480
481 test:
482 args: -use_petsc_lu
483 output_file: output/ex52_2.out
484
485 testset:
486 suffix: mumps
487 nsize: 3
488 requires: mumps
489 output_file: output/ex52_1.out
490 args: -use_mumps_lu
491
492 test:
493 suffix: mumps
494
495 test:
496 suffix: mumps_ooc
497 args: -mat_mumps_icntl_22 1 -mat_mumps_ooc_tmpdir /tmp
498
499 test:
500 suffix: mumps_ooc_api
501 args: -mat_mumps_icntl_22 1 -test_mumps_ooc_api
502 filter: grep -v "OOC_TMPDIR"
503
504 test:
505 suffix: mumps_2
506 nsize: 3
507 requires: mumps
508 args: -use_mumps_ch
509 output_file: output/ex52_1.out
510
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
517
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
524
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
531
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
538
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
547
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
555
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
562
563 test:
564 suffix: aij
565 args: -mat_type aij
566
567 test:
568 requires: kokkos_kernels
569 suffix: kok
570 args: -mat_type aijkokkos
571
572 test:
573 requires: cuda
574 suffix: cuda
575 args: -mat_type aijcusparse
576
577 test:
578 requires: hip
579 suffix: hip
580 args: -mat_type aijhipsparse
581
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
588
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
595
596 test:
597 suffix: aij
598 args: -mat_type aij
599
600 test:
601 requires: kokkos_kernels
602 suffix: kok
603 args: -mat_type aijkokkos
604
605 test:
606 requires: cuda
607 suffix: cuda
608 args: -mat_type aijcusparse
609
610 test:
611 requires: hip
612 suffix: hip
613 args: -mat_type aijhipsparse
614
615 test:
616 suffix: superlu_ilu
617 requires: superlu superlu_dist
618 args: -use_superlu_ilu
619 output_file: output/ex52_2.out
620
621 TEST*/
622