Actual source code: strumpack.c
1: #include <../src/mat/impls/aij/seq/aij.h>
2: #include <../src/mat/impls/aij/mpi/mpiaij.h>
3: #include <StrumpackSparseSolver.h>
5: static PetscErrorCode MatGetDiagonal_STRUMPACK(Mat A, Vec v)
6: {
7: PetscFunctionBegin;
8: SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Mat type: STRUMPACK factor");
9: PetscFunctionReturn(PETSC_SUCCESS);
10: }
12: static PetscErrorCode MatDestroy_STRUMPACK(Mat A)
13: {
14: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)A->data;
16: PetscFunctionBegin;
17: /* Deallocate STRUMPACK storage */
18: PetscStackCallExternalVoid("STRUMPACK_destroy", STRUMPACK_destroy(S));
19: PetscCall(PetscFree(A->data));
21: /* clear composed functions */
22: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL));
23: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetReordering_C", NULL));
24: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKGetReordering_C", NULL));
25: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetColPerm_C", NULL));
26: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKGetColPerm_C", NULL));
27: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetGeometricNxyz_C", NULL));
28: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetGeometricComponents_C", NULL));
29: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetGeometricWidth_C", NULL));
30: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetGPU_C", NULL));
31: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKGetGPU_C", NULL));
32: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetCompression_C", NULL));
33: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKGetCompression_C", NULL));
34: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetCompRelTol_C", NULL));
35: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKGetCompRelTol_C", NULL));
36: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetCompAbsTol_C", NULL));
37: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKGetCompAbsTol_C", NULL));
38: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetCompMaxRank_C", NULL));
39: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKGetCompMaxRank_C", NULL));
40: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetCompLeafSize_C", NULL));
41: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKGetCompLeafSize_C", NULL));
42: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetCompMinSepSize_C", NULL));
43: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKGetCompMinSepSize_C", NULL));
44: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetCompLossyPrecision_C", NULL));
45: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKGetCompLossyPrecision_C", NULL));
46: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetCompButterflyLevels_C", NULL));
47: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKGetCompButterflyLevels_C", NULL));
49: PetscFunctionReturn(PETSC_SUCCESS);
50: }
52: static PetscErrorCode MatSTRUMPACKSetReordering_STRUMPACK(Mat F, MatSTRUMPACKReordering reordering)
53: {
54: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
56: PetscFunctionBegin;
57: PetscStackCallExternalVoid("STRUMPACK_set_reordering_method", STRUMPACK_set_reordering_method(*S, (STRUMPACK_REORDERING_STRATEGY)reordering));
58: PetscFunctionReturn(PETSC_SUCCESS);
59: }
60: static PetscErrorCode MatSTRUMPACKGetReordering_STRUMPACK(Mat F, MatSTRUMPACKReordering *reordering)
61: {
62: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
64: PetscFunctionBegin;
65: PetscStackCallExternalVoid("STRUMPACK_reordering_method", *reordering = (MatSTRUMPACKReordering)STRUMPACK_reordering_method(*S));
66: PetscFunctionReturn(PETSC_SUCCESS);
67: }
69: /*@
70: MatSTRUMPACKSetReordering - Set STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master> fill-reducing reordering
72: Logically Collective
74: Input Parameters:
75: + F - the factored matrix obtained by calling `MatGetFactor()` from PETSc-STRUMPACK interface
76: - reordering - the code to be used to find the fill-reducing reordering
78: Options Database Key:
79: . -mat_strumpack_reordering <METIS> - Sparsity reducing matrix reordering, see `MatSTRUMPACKReordering`
81: Level: intermediate
83: .seealso: `MATSOLVERSTRUMPACK`, [](ch_matrices), `Mat`, `MatSTRUMPACKReordering`, `MatGetFactor()`, `MatSTRUMPACKSetColPerm()`, `MatSTRUMPACKGetReordering()`
84: @*/
85: PetscErrorCode MatSTRUMPACKSetReordering(Mat F, MatSTRUMPACKReordering reordering)
86: {
87: PetscFunctionBegin;
90: PetscTryMethod(F, "MatSTRUMPACKSetReordering_C", (Mat, MatSTRUMPACKReordering), (F, reordering));
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
93: /*@
94: MatSTRUMPACKGetReordering - Get STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master> fill-reducing reordering
96: Logically Collective
98: Input Parameters:
99: . F - the factored matrix obtained by calling `MatGetFactor()` from PETSc-STRUMPACK interface
101: Output Parameter:
102: . reordering - the code to be used to find the fill-reducing reordering
104: Level: intermediate
106: .seealso: `MATSOLVERSTRUMPACK`, [](ch_matrices), `Mat`, `MatSTRUMPACKReordering`, `MatGetFactor()`, `MatSTRUMPACKSetReordering()`, `MatSTRUMPACKSetColPerm()`
107: @*/
108: PetscErrorCode MatSTRUMPACKGetReordering(Mat F, MatSTRUMPACKReordering *reordering)
109: {
110: PetscFunctionBegin;
112: PetscTryMethod(F, "MatSTRUMPACKGetReordering_C", (Mat, MatSTRUMPACKReordering *), (F, reordering));
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: static PetscErrorCode MatSTRUMPACKSetColPerm_STRUMPACK(Mat F, PetscBool cperm)
118: {
119: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
121: PetscFunctionBegin;
122: PetscStackCallExternalVoid("STRUMPACK_set_matching", STRUMPACK_set_matching(*S, cperm ? STRUMPACK_MATCHING_MAX_DIAGONAL_PRODUCT_SCALING : STRUMPACK_MATCHING_NONE));
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
125: static PetscErrorCode MatSTRUMPACKGetColPerm_STRUMPACK(Mat F, PetscBool *cperm)
126: {
127: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
129: PetscFunctionBegin;
130: PetscStackCallExternalVoid("STRUMPACK_matching", *cperm = (PetscBool)(STRUMPACK_matching(*S) != STRUMPACK_MATCHING_NONE));
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
134: /*@
135: MatSTRUMPACKSetColPerm - Set whether STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master>
136: should try to permute the columns of the matrix in order to get a nonzero diagonal
138: Logically Collective
140: Input Parameters:
141: + F - the factored matrix obtained by calling `MatGetFactor()`
142: - cperm - `PETSC_TRUE` to permute (internally) the columns of the matrix
144: Options Database Key:
145: . -mat_strumpack_colperm <cperm> - true to use the permutation
147: Level: intermediate
149: .seealso: `MATSOLVERSTRUMPACK`, [](ch_matrices), `MatSTRUMPACKSetReordering()`, `Mat`, `MatGetFactor()`, `MatSTRUMPACKGetColPerm()`
150: @*/
151: PetscErrorCode MatSTRUMPACKSetColPerm(Mat F, PetscBool cperm)
152: {
153: PetscFunctionBegin;
156: PetscTryMethod(F, "MatSTRUMPACKSetColPerm_C", (Mat, PetscBool), (F, cperm));
157: PetscFunctionReturn(PETSC_SUCCESS);
158: }
159: /*@
160: MatSTRUMPACKGetColPerm - Get whether STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master>
161: will try to permute the columns of the matrix in order to get a nonzero diagonal
163: Logically Collective
165: Input Parameters:
166: . F - the factored matrix obtained by calling `MatGetFactor()`
168: Output Parameter:
169: . cperm - Indicates whether STRUMPACK will permute columns
171: Level: intermediate
173: .seealso: `MATSOLVERSTRUMPACK`, [](ch_matrices), `MatSTRUMPACKSetReordering()`, `Mat`, `MatGetFactor()`, `MatSTRUMPACKSetColPerm()`
174: @*/
175: PetscErrorCode MatSTRUMPACKGetColPerm(Mat F, PetscBool *cperm)
176: {
177: PetscFunctionBegin;
179: PetscTryMethod(F, "MatSTRUMPACKGetColPerm_C", (Mat, PetscBool *), (F, cperm));
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: static PetscErrorCode MatSTRUMPACKSetGPU_STRUMPACK(Mat F, PetscBool gpu)
185: {
186: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
188: PetscFunctionBegin;
189: if (gpu) {
190: #if !(defined(STRUMPACK_USE_CUDA) || defined(STRUMPACK_USE_HIP) || defined(STRUMPACK_USE_SYCL))
191: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Warning: strumpack was not configured with GPU support\n"));
192: #endif
193: PetscStackCallExternalVoid("STRUMPACK_enable_gpu", STRUMPACK_enable_gpu(*S));
194: } else PetscStackCallExternalVoid("STRUMPACK_disable_gpu", STRUMPACK_disable_gpu(*S));
195: PetscFunctionReturn(PETSC_SUCCESS);
196: }
197: static PetscErrorCode MatSTRUMPACKGetGPU_STRUMPACK(Mat F, PetscBool *gpu)
198: {
199: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
201: PetscFunctionBegin;
202: PetscStackCallExternalVoid("STRUMPACK_use_gpu", *gpu = (PetscBool)STRUMPACK_use_gpu(*S));
203: PetscFunctionReturn(PETSC_SUCCESS);
204: }
206: /*@
207: MatSTRUMPACKSetGPU - Set whether STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master>
208: should enable GPU acceleration (not supported for all compression types)
210: Logically Collective
212: Input Parameters:
213: + F - the factored matrix obtained by calling `MatGetFactor()` from PETSc-STRUMPACK interface
214: - gpu - whether or not to use GPU acceleration
216: Options Database Key:
217: . -mat_strumpack_gpu <gpu> - true to use gpu offload
219: Level: intermediate
221: .seealso: `MATSOLVERSTRUMPACK`, `MatGetFactor()`, `MatSTRUMPACKGetGPU()`
222: @*/
223: PetscErrorCode MatSTRUMPACKSetGPU(Mat F, PetscBool gpu)
224: {
225: PetscFunctionBegin;
228: PetscTryMethod(F, "MatSTRUMPACKSetGPU_C", (Mat, PetscBool), (F, gpu));
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
231: /*@
232: MatSTRUMPACKGetGPU - Get whether STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master>
233: will try to use GPU acceleration (not supported for all compression types)
235: Logically Collective
237: Input Parameters:
238: . F - the factored matrix obtained by calling `MatGetFactor()` from PETSc-STRUMPACK interface
240: Output Parameter:
241: . gpu - whether or not STRUMPACK will try to use GPU acceleration
243: Level: intermediate
245: .seealso: `MATSOLVERSTRUMPACK`, `MatGetFactor()`, `MatSTRUMPACKSetGPU()`
246: @*/
247: PetscErrorCode MatSTRUMPACKGetGPU(Mat F, PetscBool *gpu)
248: {
249: PetscFunctionBegin;
251: PetscTryMethod(F, "MatSTRUMPACKGetGPU_C", (Mat, PetscBool *), (F, gpu));
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: static PetscErrorCode MatSTRUMPACKSetCompression_STRUMPACK(Mat F, MatSTRUMPACKCompressionType comp)
257: {
258: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
260: PetscFunctionBegin;
261: #if !defined(STRUMPACK_USE_BPACK)
262: PetscCheck(comp != MAT_STRUMPACK_COMPRESSION_TYPE_HODLR && comp != MAT_STRUMPACK_COMPRESSION_TYPE_BLR_HODLR && comp != MAT_STRUMPACK_COMPRESSION_TYPE_ZFP_BLR_HODLR, PetscObjectComm((PetscObject)F), PETSC_ERR_SUP, "Compression scheme requires ButterflyPACK, please reconfigure with --download-butterflypack");
263: #endif
264: #if !defined(STRUMPACK_USE_ZFP)
265: PetscCheck(comp != MAT_STRUMPACK_COMPRESSION_TYPE_ZFP_BLR_HODLR && comp != MAT_STRUMPACK_COMPRESSION_TYPE_LOSSLESS && comp != MAT_STRUMPACK_COMPRESSION_TYPE_LOSSY, PetscObjectComm((PetscObject)F), PETSC_ERR_SUP, "Compression scheme requires ZFP, please reconfigure with --download-zfp");
266: #endif
267: PetscStackCallExternalVoid("STRUMPACK_set_compression", STRUMPACK_set_compression(*S, (STRUMPACK_COMPRESSION_TYPE)comp));
268: PetscFunctionReturn(PETSC_SUCCESS);
269: }
270: static PetscErrorCode MatSTRUMPACKGetCompression_STRUMPACK(Mat F, MatSTRUMPACKCompressionType *comp)
271: {
272: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
274: PetscFunctionBegin;
275: PetscStackCallExternalVoid("STRUMPACK_compression", *comp = (MatSTRUMPACKCompressionType)STRUMPACK_compression(*S));
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: /*@
280: MatSTRUMPACKSetCompression - Set STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master> compression type
282: Input Parameters:
283: + F - the factored matrix obtained by calling `MatGetFactor()` from PETSc-STRUMPACK interface
284: - comp - Type of compression to be used in the approximate sparse factorization
286: Options Database Key:
287: . -mat_strumpack_compression <NONE> - Type of rank-structured compression in sparse LU factors (choose one of) NONE HSS BLR HODLR BLR_HODLR ZFP_BLR_HODLR LOSSLESS LOSSY
289: Level: intermediate
291: Note:
292: Default for `comp` is `MAT_STRUMPACK_COMPRESSION_TYPE_NONE` for `-pc_type lu` and `MAT_STRUMPACK_COMPRESSION_TYPE_BLR`
293: for `-pc_type ilu`
295: .seealso: `MATSOLVERSTRUMPACK`, `MatGetFactor()`, `MatSTRUMPACKCompressionType`, `MatSTRUMPACKGetCompression()`
296: @*/
297: PetscErrorCode MatSTRUMPACKSetCompression(Mat F, MatSTRUMPACKCompressionType comp)
298: {
299: PetscFunctionBegin;
302: PetscTryMethod(F, "MatSTRUMPACKSetCompression_C", (Mat, MatSTRUMPACKCompressionType), (F, comp));
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
305: /*@
306: MatSTRUMPACKGetCompression - Get STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master> compression type
308: Input Parameters:
309: . F - the factored matrix obtained by calling `MatGetFactor()` from PETSc-STRUMPACK interface
311: Output Parameter:
312: . comp - Type of compression to be used in the approximate sparse factorization
314: Level: intermediate
316: Note:
317: Default is `MAT_STRUMPACK_COMPRESSION_TYPE_NONE` for `-pc_type lu` and `MAT_STRUMPACK_COMPRESSION_TYPE_BLR` for `-pc_type ilu`
319: .seealso: `MATSOLVERSTRUMPACK`, `MatGetFactor()`, `MatSTRUMPACKCompressionType`, `MatSTRUMPACKSetCompression()`
320: @*/
321: PetscErrorCode MatSTRUMPACKGetCompression(Mat F, MatSTRUMPACKCompressionType *comp)
322: {
323: PetscFunctionBegin;
325: PetscTryMethod(F, "MatSTRUMPACKGetCompression_C", (Mat, MatSTRUMPACKCompressionType *), (F, comp));
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
330: static PetscErrorCode MatSTRUMPACKSetCompRelTol_STRUMPACK(Mat F, PetscReal rtol)
331: {
332: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
334: PetscFunctionBegin;
335: PetscStackCallExternalVoid("STRUMPACK_set_compression_rel_tol", STRUMPACK_set_compression_rel_tol(*S, rtol));
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
338: static PetscErrorCode MatSTRUMPACKGetCompRelTol_STRUMPACK(Mat F, PetscReal *rtol)
339: {
340: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
342: PetscFunctionBegin;
343: PetscStackCallExternalVoid("STRUMPACK_compression_rel_tol", *rtol = (PetscReal)STRUMPACK_compression_rel_tol(*S));
344: PetscFunctionReturn(PETSC_SUCCESS);
345: }
347: /*@
348: MatSTRUMPACKSetCompRelTol - Set STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master> relative tolerance for compression
350: Logically Collective
352: Input Parameters:
353: + F - the factored matrix obtained by calling `MatGetFactor()`
354: - rtol - relative compression tolerance
356: Options Database Key:
357: . -mat_strumpack_compression_rel_tol <1e-4> - Relative compression tolerance, when using `-pctype ilu`
359: Level: intermediate
361: .seealso: `MATSOLVERSTRUMPACK`, [](ch_matrices), `Mat`, `MatGetFactor()`, `MatSTRUMPACKGetCompRelTol()`, `MatSTRUMPACKSetReordering()`, `MatSTRUMPACKSetColPerm()`
362: @*/
363: PetscErrorCode MatSTRUMPACKSetCompRelTol(Mat F, PetscReal rtol)
364: {
365: PetscFunctionBegin;
368: PetscTryMethod(F, "MatSTRUMPACKSetCompRelTol_C", (Mat, PetscReal), (F, rtol));
369: PetscFunctionReturn(PETSC_SUCCESS);
370: }
371: /*@
372: MatSTRUMPACKGetCompRelTol - Get STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master> relative tolerance for compression
374: Logically Collective
376: Input Parameters:
377: . F - the factored matrix obtained by calling `MatGetFactor()`
379: Output Parameter:
380: . rtol - relative compression tolerance
382: Level: intermediate
384: .seealso: `MATSOLVERSTRUMPACK`, [](ch_matrices), `Mat`, `MatGetFactor()`, `MatSTRUMPACKSetCompRelTol()`, `MatSTRUMPACKSetReordering()`, `MatSTRUMPACKSetColPerm()`
385: @*/
386: PetscErrorCode MatSTRUMPACKGetCompRelTol(Mat F, PetscReal *rtol)
387: {
388: PetscFunctionBegin;
390: PetscTryMethod(F, "MatSTRUMPACKGetCompRelTol_C", (Mat, PetscReal *), (F, rtol));
392: PetscFunctionReturn(PETSC_SUCCESS);
393: }
395: static PetscErrorCode MatSTRUMPACKSetCompAbsTol_STRUMPACK(Mat F, PetscReal atol)
396: {
397: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
399: PetscFunctionBegin;
400: PetscStackCallExternalVoid("STRUMPACK_set_compression_abs_tol", STRUMPACK_set_compression_abs_tol(*S, atol));
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
403: static PetscErrorCode MatSTRUMPACKGetCompAbsTol_STRUMPACK(Mat F, PetscReal *atol)
404: {
405: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
407: PetscFunctionBegin;
408: PetscStackCallExternalVoid("STRUMPACK_compression_abs_tol", *atol = (PetscReal)STRUMPACK_compression_abs_tol(*S));
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }
412: /*@
413: MatSTRUMPACKSetCompAbsTol - Set STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master> absolute tolerance for compression
415: Logically Collective
417: Input Parameters:
418: + F - the factored matrix obtained by calling `MatGetFactor()`
419: - atol - absolute compression tolerance
421: Options Database Key:
422: . -mat_strumpack_compression_abs_tol <1e-10> - Absolute compression tolerance, when using `-pctype ilu`
424: Level: intermediate
426: .seealso: `MATSOLVERSTRUMPACK`, [](ch_matrices), `Mat`, `MatGetFactor()`, `MatSTRUMPACKGetCompAbsTol()`, `MatSTRUMPACKSetReordering()`, `MatSTRUMPACKSetColPerm()`
427: @*/
428: PetscErrorCode MatSTRUMPACKSetCompAbsTol(Mat F, PetscReal atol)
429: {
430: PetscFunctionBegin;
433: PetscTryMethod(F, "MatSTRUMPACKSetCompAbsTol_C", (Mat, PetscReal), (F, atol));
434: PetscFunctionReturn(PETSC_SUCCESS);
435: }
436: /*@
437: MatSTRUMPACKGetCompAbsTol - Get STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master> absolute tolerance for compression
439: Logically Collective
441: Input Parameters:
442: . F - the factored matrix obtained by calling `MatGetFactor()`
444: Output Parameter:
445: . atol - absolute compression tolerance
447: Level: intermediate
449: .seealso: `MATSOLVERSTRUMPACK`, [](ch_matrices), `Mat`, `MatGetFactor()`, `MatSTRUMPACKSetCompAbsTol()`, `MatSTRUMPACKSetReordering()`, `MatSTRUMPACKSetColPerm()`
450: @*/
451: PetscErrorCode MatSTRUMPACKGetCompAbsTol(Mat F, PetscReal *atol)
452: {
453: PetscFunctionBegin;
455: PetscTryMethod(F, "MatSTRUMPACKGetCompAbsTol_C", (Mat, PetscReal *), (F, atol));
457: PetscFunctionReturn(PETSC_SUCCESS);
458: }
460: static PetscErrorCode MatSTRUMPACKSetCompLeafSize_STRUMPACK(Mat F, PetscInt leaf_size)
461: {
462: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
464: PetscFunctionBegin;
465: PetscStackCallExternalVoid("STRUMPACK_set_compression_leaf_size", STRUMPACK_set_compression_leaf_size(*S, leaf_size));
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
468: static PetscErrorCode MatSTRUMPACKGetCompLeafSize_STRUMPACK(Mat F, PetscInt *leaf_size)
469: {
470: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
472: PetscFunctionBegin;
473: PetscStackCallExternalVoid("STRUMPACK_compression_leaf_size", *leaf_size = (PetscInt)STRUMPACK_compression_leaf_size(*S));
474: PetscFunctionReturn(PETSC_SUCCESS);
475: }
477: /*@
478: MatSTRUMPACKSetCompLeafSize - Set STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master> leaf size for HSS, BLR, HODLR...
480: Logically Collective
482: Input Parameters:
483: + F - the factored matrix obtained by calling `MatGetFactor()` from PETSc-STRUMPACK interface
484: - leaf_size - Size of diagonal blocks in rank-structured approximation
486: Options Database Key:
487: . -mat_strumpack_compression_leaf_size - Size of diagonal blocks in rank-structured approximation, when using `-pctype ilu`
489: Level: intermediate
491: .seealso: `MATSOLVERSTRUMPACK`, [](ch_matrices), `Mat`, `MatGetFactor()`, `MatSTRUMPACKGetCompLeafSize()`, `MatSTRUMPACKSetReordering()`, `MatSTRUMPACKSetColPerm()`
492: @*/
493: PetscErrorCode MatSTRUMPACKSetCompLeafSize(Mat F, PetscInt leaf_size)
494: {
495: PetscFunctionBegin;
498: PetscTryMethod(F, "MatSTRUMPACKSetCompLeafSize_C", (Mat, PetscInt), (F, leaf_size));
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
501: /*@
502: MatSTRUMPACKGetCompLeafSize - Get STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master> leaf size for HSS, BLR, HODLR...
504: Logically Collective
506: Input Parameters:
507: . F - the factored matrix obtained by calling `MatGetFactor()` from PETSc-STRUMPACK interface
509: Output Parameter:
510: . leaf_size - Size of diagonal blocks in rank-structured approximation
512: Level: intermediate
514: .seealso: `MATSOLVERSTRUMPACK`, [](ch_matrices), `Mat`, `MatGetFactor()`, `MatSTRUMPACKSetCompLeafSize()`, `MatSTRUMPACKSetReordering()`, `MatSTRUMPACKSetColPerm()`
515: @*/
516: PetscErrorCode MatSTRUMPACKGetCompLeafSize(Mat F, PetscInt *leaf_size)
517: {
518: PetscFunctionBegin;
520: PetscTryMethod(F, "MatSTRUMPACKGetCompLeafSize_C", (Mat, PetscInt *), (F, leaf_size));
522: PetscFunctionReturn(PETSC_SUCCESS);
523: }
525: static PetscErrorCode MatSTRUMPACKSetGeometricNxyz_STRUMPACK(Mat F, PetscInt nx, PetscInt ny, PetscInt nz)
526: {
527: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
529: PetscFunctionBegin;
530: if (nx < 1) {
531: PetscCheck(nx == PETSC_DECIDE || nx == PETSC_DEFAULT, PetscObjectComm((PetscObject)F), PETSC_ERR_ARG_OUTOFRANGE, "nx < 1");
532: nx = 1;
533: }
534: PetscStackCallExternalVoid("STRUMPACK_set_nx", STRUMPACK_set_nx(*S, nx));
535: if (ny < 1) {
536: PetscCheck(ny == PETSC_DECIDE || ny == PETSC_DEFAULT, PetscObjectComm((PetscObject)F), PETSC_ERR_ARG_OUTOFRANGE, "ny < 1");
537: ny = 1;
538: }
539: PetscStackCallExternalVoid("STRUMPACK_set_ny", STRUMPACK_set_ny(*S, ny));
540: if (nz < 1) {
541: PetscCheck(nz == PETSC_DECIDE || nz == PETSC_DEFAULT, PetscObjectComm((PetscObject)F), PETSC_ERR_ARG_OUTOFRANGE, "nz < 1");
542: nz = 1;
543: }
544: PetscStackCallExternalVoid("STRUMPACK_set_nz", STRUMPACK_set_nz(*S, nz));
545: PetscFunctionReturn(PETSC_SUCCESS);
546: }
547: static PetscErrorCode MatSTRUMPACKSetGeometricComponents_STRUMPACK(Mat F, PetscInt nc)
548: {
549: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
551: PetscFunctionBegin;
552: PetscStackCallExternalVoid("STRUMPACK_set_components", STRUMPACK_set_components(*S, nc));
553: PetscFunctionReturn(PETSC_SUCCESS);
554: }
555: static PetscErrorCode MatSTRUMPACKSetGeometricWidth_STRUMPACK(Mat F, PetscInt w)
556: {
557: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
559: PetscFunctionBegin;
560: PetscStackCallExternalVoid("STRUMPACK_set_separator_width", STRUMPACK_set_separator_width(*S, w));
561: PetscFunctionReturn(PETSC_SUCCESS);
562: }
564: /*@
565: MatSTRUMPACKSetGeometricNxyz - Set STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master> mesh x, y and z dimensions, for use with GEOMETRIC ordering.
567: Logically Collective
569: Input Parameters:
570: + F - the factored matrix obtained by calling `MatGetFactor()` from PETSc-STRUMPACK interface
571: . nx - x dimension of the mesh
572: . ny - y dimension of the mesh
573: - nz - z dimension of the mesh
575: Level: intermediate
577: Note:
578: If the mesh is two (or one) dimensional one can use 1, `PETSC_DECIDE` or `PETSC_DEFAULT`
579: for the missing z (and y) dimensions.
581: .seealso: `MATSOLVERSTRUMPACK`, `MatGetFactor()`
582: @*/
583: PetscErrorCode MatSTRUMPACKSetGeometricNxyz(Mat F, PetscInt nx, PetscInt ny, PetscInt nz)
584: {
585: PetscFunctionBegin;
590: PetscTryMethod(F, "MatSTRUMPACKSetGeometricNxyz_C", (Mat, PetscInt, PetscInt, PetscInt), (F, nx, ny, nz));
591: PetscFunctionReturn(PETSC_SUCCESS);
592: }
593: /*@
594: MatSTRUMPACKSetGeometricComponents - Set STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master>
595: number of degrees of freedom per mesh point, for use with GEOMETRIC ordering.
597: Logically Collective
599: Input Parameters:
600: + F - the factored matrix obtained by calling `MatGetFactor()` from PETSc-STRUMPACK interface
601: - nc - Number of components/dof's per grid point
603: Options Database Key:
604: . -mat_strumpack_geometric_components <1> - Number of components per mesh point, for geometric nested dissection ordering
606: Level: intermediate
608: .seealso: `MATSOLVERSTRUMPACK`, `MatGetFactor()`
609: @*/
610: PetscErrorCode MatSTRUMPACKSetGeometricComponents(Mat F, PetscInt nc)
611: {
612: PetscFunctionBegin;
615: PetscTryMethod(F, "MatSTRUMPACKSetGeometricComponents_C", (Mat, PetscInt), (F, nc));
616: PetscFunctionReturn(PETSC_SUCCESS);
617: }
618: /*@
619: MatSTRUMPACKSetGeometricWidth - Set STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master> width of the separator, for use with GEOMETRIC ordering.
621: Logically Collective
623: Input Parameters:
624: + F - the factored matrix obtained by calling `MatGetFactor()` from PETSc-STRUMPACK interface
625: - w - width of the separator
627: Options Database Key:
628: . -mat_strumpack_geometric_width <1> - Width of the separator of the mesh, for geometric nested dissection ordering
630: Level: intermediate
632: .seealso: `MATSOLVERSTRUMPACK`, `MatGetFactor()`
633: @*/
634: PetscErrorCode MatSTRUMPACKSetGeometricWidth(Mat F, PetscInt w)
635: {
636: PetscFunctionBegin;
639: PetscTryMethod(F, "MatSTRUMPACKSetGeometricWidth_C", (Mat, PetscInt), (F, w));
640: PetscFunctionReturn(PETSC_SUCCESS);
641: }
643: static PetscErrorCode MatSTRUMPACKSetCompMinSepSize_STRUMPACK(Mat F, PetscInt min_sep_size)
644: {
645: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
647: PetscFunctionBegin;
648: PetscStackCallExternalVoid("STRUMPACK_set_compression_min_sep_size", STRUMPACK_set_compression_min_sep_size(*S, min_sep_size));
649: PetscFunctionReturn(PETSC_SUCCESS);
650: }
651: static PetscErrorCode MatSTRUMPACKGetCompMinSepSize_STRUMPACK(Mat F, PetscInt *min_sep_size)
652: {
653: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
655: PetscFunctionBegin;
656: PetscStackCallExternalVoid("STRUMPACK_compression_min_sep_size", *min_sep_size = (PetscInt)STRUMPACK_compression_min_sep_size(*S));
657: PetscFunctionReturn(PETSC_SUCCESS);
658: }
660: /*@
661: MatSTRUMPACKSetCompMinSepSize - Set STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master> minimum separator size for low-rank approximation
663: Logically Collective
665: Input Parameters:
666: + F - the factored matrix obtained by calling `MatGetFactor()` from PETSc-STRUMPACK interface
667: - min_sep_size - minimum dense matrix size for low-rank approximation
669: Options Database Key:
670: . -mat_strumpack_compression_min_sep_size <min_sep_size> - Minimum size of dense sub-block for low-rank compression
672: Level: intermediate
674: .seealso: `MATSOLVERSTRUMPACK`, `MatGetFactor()`, `MatSTRUMPACKGetCompMinSepSize()`
675: @*/
676: PetscErrorCode MatSTRUMPACKSetCompMinSepSize(Mat F, PetscInt min_sep_size)
677: {
678: PetscFunctionBegin;
681: PetscTryMethod(F, "MatSTRUMPACKSetCompMinSepSize_C", (Mat, PetscInt), (F, min_sep_size));
682: PetscFunctionReturn(PETSC_SUCCESS);
683: }
684: /*@
685: MatSTRUMPACKGetCompMinSepSize - Get STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master> minimum separator size for low-rank approximation
687: Logically Collective
689: Input Parameters:
690: . F - the factored matrix obtained by calling `MatGetFactor()` from PETSc-STRUMPACK interface
692: Output Parameter:
693: . min_sep_size - minimum dense matrix size for low-rank approximation
695: Level: intermediate
697: .seealso: `MATSOLVERSTRUMPACK`, `MatGetFactor()`, `MatSTRUMPACKSetCompMinSepSize()`
698: @*/
699: PetscErrorCode MatSTRUMPACKGetCompMinSepSize(Mat F, PetscInt *min_sep_size)
700: {
701: PetscFunctionBegin;
703: PetscTryMethod(F, "MatSTRUMPACKGetCompMinSepSize_C", (Mat, PetscInt *), (F, min_sep_size));
705: PetscFunctionReturn(PETSC_SUCCESS);
706: }
708: static PetscErrorCode MatSTRUMPACKSetCompLossyPrecision_STRUMPACK(Mat F, PetscInt lossy_prec)
709: {
710: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
712: PetscFunctionBegin;
713: PetscStackCallExternalVoid("STRUMPACK_set_compression_lossy_precision", STRUMPACK_set_compression_lossy_precision(*S, lossy_prec));
714: PetscFunctionReturn(PETSC_SUCCESS);
715: }
716: static PetscErrorCode MatSTRUMPACKGetCompLossyPrecision_STRUMPACK(Mat F, PetscInt *lossy_prec)
717: {
718: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
720: PetscFunctionBegin;
721: PetscStackCallExternalVoid("STRUMPACK_compression_lossy_precision", *lossy_prec = (PetscInt)STRUMPACK_compression_lossy_precision(*S));
722: PetscFunctionReturn(PETSC_SUCCESS);
723: }
725: /*@
726: MatSTRUMPACKSetCompLossyPrecision - Set STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master> precision for lossy compression (requires ZFP support)
728: Logically Collective
730: Input Parameters:
731: + F - the factored matrix obtained by calling `MatGetFactor()` from PETSc-STRUMPACK interface
732: - lossy_prec - Number of bitplanes to use in lossy compression
734: Options Database Key:
735: . -mat_strumpack_compression_lossy_precision <lossy_prec> - Precision when using lossy compression [1-64], when using `-pctype ilu -mat_strumpack_compression MAT_STRUMPACK_COMPRESSION_TYPE_LOSSY`
737: Level: intermediate
739: .seealso: `MATSOLVERSTRUMPACK`, `MatGetFactor()`, `MatSTRUMPACKGetCompLossyPrecision()`
740: @*/
741: PetscErrorCode MatSTRUMPACKSetCompLossyPrecision(Mat F, PetscInt lossy_prec)
742: {
743: PetscFunctionBegin;
746: PetscTryMethod(F, "MatSTRUMPACKSetCompLossyPrecision_C", (Mat, PetscInt), (F, lossy_prec));
747: PetscFunctionReturn(PETSC_SUCCESS);
748: }
749: /*@
750: MatSTRUMPACKGetCompLossyPrecision - Get STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master> precision for lossy compression (requires ZFP support)
752: Logically Collective
754: Input Parameters:
755: . F - the factored matrix obtained by calling `MatGetFactor()` from PETSc-STRUMPACK interface
757: Output Parameter:
758: . lossy_prec - Number of bitplanes to use in lossy compression
760: Level: intermediate
762: .seealso: `MATSOLVERSTRUMPACK`, `MatGetFactor()`, `MatSTRUMPACKSetCompLossyPrecision()`
763: @*/
764: PetscErrorCode MatSTRUMPACKGetCompLossyPrecision(Mat F, PetscInt *lossy_prec)
765: {
766: PetscFunctionBegin;
768: PetscTryMethod(F, "MatSTRUMPACKGetCompLossyPrecision_C", (Mat, PetscInt *), (F, lossy_prec));
770: PetscFunctionReturn(PETSC_SUCCESS);
771: }
773: static PetscErrorCode MatSTRUMPACKSetCompButterflyLevels_STRUMPACK(Mat F, PetscInt bfly_lvls)
774: {
775: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
777: PetscFunctionBegin;
778: PetscStackCallExternalVoid("STRUMPACK_set_compression_butterfly_levels", STRUMPACK_set_compression_butterfly_levels(*S, bfly_lvls));
779: PetscFunctionReturn(PETSC_SUCCESS);
780: }
781: static PetscErrorCode MatSTRUMPACKGetCompButterflyLevels_STRUMPACK(Mat F, PetscInt *bfly_lvls)
782: {
783: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
785: PetscFunctionBegin;
786: PetscStackCallExternalVoid("STRUMPACK_compression_butterfly_levels", *bfly_lvls = (PetscInt)STRUMPACK_compression_butterfly_levels(*S));
787: PetscFunctionReturn(PETSC_SUCCESS);
788: }
790: /*@
791: MatSTRUMPACKSetCompButterflyLevels - Set STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master>
792: number of butterfly levels in HODLR compression (requires ButterflyPACK support)
794: Logically Collective
796: Input Parameters:
797: + F - the factored matrix obtained by calling `MatGetFactor()` from PETSc-STRUMPACK interface
798: - bfly_lvls - Number of levels of butterfly compression in HODLR compression
800: Options Database Key:
801: . -mat_strumpack_compression_butterfly_levels <bfly_lvls> - Number of levels in the hierarchically off-diagonal matrix for which to use butterfly,
802: when using `-pctype ilu`, (BLR_)HODLR compression
804: Level: intermediate
806: .seealso: `MATSOLVERSTRUMPACK`, `MatGetFactor()`, `MatSTRUMPACKGetCompButterflyLevels()`
807: @*/
808: PetscErrorCode MatSTRUMPACKSetCompButterflyLevels(Mat F, PetscInt bfly_lvls)
809: {
810: PetscFunctionBegin;
813: PetscTryMethod(F, "MatSTRUMPACKSetButterflyLevels_C", (Mat, PetscInt), (F, bfly_lvls));
814: PetscFunctionReturn(PETSC_SUCCESS);
815: }
816: /*@
817: MatSTRUMPACKGetCompButterflyLevels - Get STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master>
818: number of butterfly levels in HODLR compression (requires ButterflyPACK support)
820: Logically Collective
822: Input Parameters:
823: . F - the factored matrix obtained by calling `MatGetFactor()` from PETSc-STRUMPACK interface
825: Output Parameter:
826: . bfly_lvls - Number of levels of butterfly compression in HODLR compression
828: Level: intermediate
830: .seealso: `MATSOLVERSTRUMPACK`, `MatGetFactor()`, `MatSTRUMPACKSetCompButterflyLevels()`
831: @*/
832: PetscErrorCode MatSTRUMPACKGetCompButterflyLevels(Mat F, PetscInt *bfly_lvls)
833: {
834: PetscFunctionBegin;
836: PetscTryMethod(F, "MatSTRUMPACKGetButterflyLevels_C", (Mat, PetscInt *), (F, bfly_lvls));
838: PetscFunctionReturn(PETSC_SUCCESS);
839: }
841: static PetscErrorCode MatSolve_STRUMPACK(Mat A, Vec b_mpi, Vec x)
842: {
843: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)A->data;
844: STRUMPACK_RETURN_CODE sp_err;
845: const PetscScalar *bptr;
846: PetscScalar *xptr;
848: PetscFunctionBegin;
849: PetscCall(VecGetArray(x, &xptr));
850: PetscCall(VecGetArrayRead(b_mpi, &bptr));
852: PetscStackCallExternalVoid("STRUMPACK_solve", sp_err = STRUMPACK_solve(*S, (PetscScalar *)bptr, xptr, 0));
853: switch (sp_err) {
854: case STRUMPACK_SUCCESS:
855: break;
856: case STRUMPACK_MATRIX_NOT_SET: {
857: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: matrix was not set");
858: break;
859: }
860: case STRUMPACK_REORDERING_ERROR: {
861: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: matrix reordering failed");
862: break;
863: }
864: default:
865: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: solve failed");
866: }
867: PetscCall(VecRestoreArray(x, &xptr));
868: PetscCall(VecRestoreArrayRead(b_mpi, &bptr));
869: PetscFunctionReturn(PETSC_SUCCESS);
870: }
872: static PetscErrorCode MatMatSolve_STRUMPACK(Mat A, Mat B_mpi, Mat X)
873: {
874: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)A->data;
875: STRUMPACK_RETURN_CODE sp_err;
876: PetscBool flg;
877: PetscInt m = A->rmap->n, nrhs;
878: const PetscScalar *bptr;
879: PetscScalar *xptr;
881: PetscFunctionBegin;
882: PetscCall(PetscObjectTypeCompareAny((PetscObject)B_mpi, &flg, MATSEQDENSE, MATMPIDENSE, NULL));
883: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix B must be MATDENSE matrix");
884: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &flg, MATSEQDENSE, MATMPIDENSE, NULL));
885: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix X must be MATDENSE matrix");
887: PetscCall(MatGetSize(B_mpi, NULL, &nrhs));
888: PetscCall(MatDenseGetArray(X, &xptr));
889: PetscCall(MatDenseGetArrayRead(B_mpi, &bptr));
891: PetscStackCallExternalVoid("STRUMPACK_matsolve", sp_err = STRUMPACK_matsolve(*S, nrhs, bptr, m, xptr, m, 0));
892: switch (sp_err) {
893: case STRUMPACK_SUCCESS:
894: break;
895: case STRUMPACK_MATRIX_NOT_SET: {
896: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: matrix was not set");
897: break;
898: }
899: case STRUMPACK_REORDERING_ERROR: {
900: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: matrix reordering failed");
901: break;
902: }
903: default:
904: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: solve failed");
905: }
906: PetscCall(MatDenseRestoreArrayRead(B_mpi, &bptr));
907: PetscCall(MatDenseRestoreArray(X, &xptr));
909: PetscFunctionReturn(PETSC_SUCCESS);
910: }
912: static PetscErrorCode MatView_Info_STRUMPACK(Mat A, PetscViewer viewer)
913: {
914: PetscFunctionBegin;
915: /* check if matrix is strumpack type */
916: if (A->ops->solve != MatSolve_STRUMPACK) PetscFunctionReturn(PETSC_SUCCESS);
917: PetscCall(PetscViewerASCIIPrintf(viewer, "STRUMPACK sparse solver!\n"));
918: PetscFunctionReturn(PETSC_SUCCESS);
919: }
921: static PetscErrorCode MatView_STRUMPACK(Mat A, PetscViewer viewer)
922: {
923: PetscBool iascii;
924: PetscViewerFormat format;
926: PetscFunctionBegin;
927: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
928: if (iascii) {
929: PetscCall(PetscViewerGetFormat(viewer, &format));
930: if (format == PETSC_VIEWER_ASCII_INFO) PetscCall(MatView_Info_STRUMPACK(A, viewer));
931: }
932: PetscFunctionReturn(PETSC_SUCCESS);
933: }
935: static PetscErrorCode MatLUFactorNumeric_STRUMPACK(Mat F, Mat A, const MatFactorInfo *info)
936: {
937: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
938: STRUMPACK_RETURN_CODE sp_err;
939: Mat Aloc;
940: const PetscScalar *av;
941: const PetscInt *ai = NULL, *aj = NULL;
942: PetscInt M = A->rmap->N, m = A->rmap->n, dummy;
943: PetscBool ismpiaij, isseqaij, flg;
945: PetscFunctionBegin;
946: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &isseqaij));
947: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
948: if (ismpiaij) {
949: PetscCall(MatMPIAIJGetLocalMat(A, MAT_INITIAL_MATRIX, &Aloc));
950: } else if (isseqaij) {
951: PetscCall(PetscObjectReference((PetscObject)A));
952: Aloc = A;
953: } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for type %s", ((PetscObject)A)->type_name);
955: PetscCall(MatGetRowIJ(Aloc, 0, PETSC_FALSE, PETSC_FALSE, &dummy, &ai, &aj, &flg));
956: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "GetRowIJ failed");
957: PetscCall(MatSeqAIJGetArrayRead(Aloc, &av));
959: if (ismpiaij) {
960: const PetscInt *dist = NULL;
961: PetscCall(MatGetOwnershipRanges(A, &dist));
962: PetscStackCallExternalVoid("STRUMPACK_set_distributed_csr_matrix", STRUMPACK_set_distributed_csr_matrix(*S, &m, ai, aj, av, dist, 0));
963: } else if (isseqaij) {
964: PetscStackCallExternalVoid("STRUMPACK_set_csr_matrix", STRUMPACK_set_csr_matrix(*S, &M, ai, aj, av, 0));
965: }
967: PetscCall(MatRestoreRowIJ(Aloc, 0, PETSC_FALSE, PETSC_FALSE, &dummy, &ai, &aj, &flg));
968: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "RestoreRowIJ failed");
969: PetscCall(MatSeqAIJRestoreArrayRead(Aloc, &av));
970: PetscCall(MatDestroy(&Aloc));
972: /* Reorder and Factor the matrix. */
973: /* TODO figure out how to avoid reorder if the matrix values changed, but the pattern remains the same. */
974: PetscStackCallExternalVoid("STRUMPACK_reorder", sp_err = STRUMPACK_reorder(*S));
975: PetscStackCallExternalVoid("STRUMPACK_factor", sp_err = STRUMPACK_factor(*S));
976: switch (sp_err) {
977: case STRUMPACK_SUCCESS:
978: break;
979: case STRUMPACK_MATRIX_NOT_SET: {
980: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: matrix was not set");
981: break;
982: }
983: case STRUMPACK_REORDERING_ERROR: {
984: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: matrix reordering failed");
985: break;
986: }
987: default:
988: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: factorization failed");
989: }
990: F->assembled = PETSC_TRUE;
991: F->preallocated = PETSC_TRUE;
992: PetscFunctionReturn(PETSC_SUCCESS);
993: }
995: static PetscErrorCode MatLUFactorSymbolic_STRUMPACK(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
996: {
997: PetscFunctionBegin;
998: F->ops->lufactornumeric = MatLUFactorNumeric_STRUMPACK;
999: F->ops->solve = MatSolve_STRUMPACK;
1000: F->ops->matsolve = MatMatSolve_STRUMPACK;
1001: PetscFunctionReturn(PETSC_SUCCESS);
1002: }
1004: static PetscErrorCode MatFactorGetSolverType_aij_strumpack(Mat A, MatSolverType *type)
1005: {
1006: PetscFunctionBegin;
1007: *type = MATSOLVERSTRUMPACK;
1008: PetscFunctionReturn(PETSC_SUCCESS);
1009: }
1011: /*MC
1012: MATSOLVERSTRUMPACK = "strumpack" - A solver package providing a direct sparse solver (PCLU)
1013: and a preconditioner (PCILU) using low-rank compression via the external package STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master>.
1015: Use `./configure --download-strumpack --download-metis` to have PETSc installed with STRUMPACK.
1017: For full functionality, add `--download-slate --download-magma --download-parmetis --download-ptscotch --download-zfp --download-butterflypack`.
1018: SLATE provides GPU support in the multi-GPU setting, providing ScaLAPACK functionality but with GPU acceleration.
1019: MAGMA can optionally be used for on node GPU support instead cuBLAS/cuSOLVER, and performs slightly better.
1020: ParMETIS and PTScotch can be used for parallel fill-reducing ordering.
1021: ZFP is used for floating point compression of the sparse factors (LOSSY or LOSSLESS compression).
1022: ButterflyPACK is used for HODLR (Hierarchically Off-Diagonal Low Rank) and HODBF (Hierarchically Off-Diagonal Butterfly) compression of the sparse factors.
1024: Options Database Keys:
1025: + -mat_strumpack_verbose - Enable verbose output
1026: . -mat_strumpack_compression - Type of rank-structured compression in sparse LU factors (choose one of) NONE HSS BLR HODLR BLR_HODLR ZFP_BLR_HODLR LOSSLESS LOSSY
1027: . -mat_strumpack_compression_rel_tol - Relative compression tolerance, when using `-pctype ilu`
1028: . -mat_strumpack_compression_abs_tol - Absolute compression tolerance, when using `-pctype ilu`
1029: . -mat_strumpack_compression_min_sep_size - Minimum size of separator for rank-structured compression, when using `-pctype ilu`
1030: . -mat_strumpack_compression_leaf_size - Size of diagonal blocks in rank-structured approximation, when using `-pctype ilu`
1031: . -mat_strumpack_compression_lossy_precision - Precision when using lossy compression [1-64], when using `-pctype ilu`, compression LOSSY (requires ZFP support)
1032: . -mat_strumpack_compression_butterfly_levels - Number of levels in the hierarchically off-diagonal matrix for which to use butterfly, when using `-pctype ilu`, (BLR_)HODLR compression (requires ButterflyPACK support)
1033: . -mat_strumpack_gpu - Enable GPU acceleration in numerical factorization (not supported for all compression types)
1034: . -mat_strumpack_colperm <TRUE> - Permute matrix to make diagonal nonzeros
1035: . -mat_strumpack_reordering <METIS> - Sparsity reducing matrix reordering (choose one of) NATURAL METIS PARMETIS SCOTCH PTSCOTCH RCM GEOMETRIC AMD MMD AND MLF SPECTRAL
1036: . -mat_strumpack_geometric_xyz <1,1,1> - Mesh x,y,z dimensions, for use with GEOMETRIC ordering
1037: . -mat_strumpack_geometric_components <1> - Number of components per mesh point, for geometric nested dissection ordering
1038: . -mat_strumpack_geometric_width <1> - Width of the separator of the mesh, for geometric nested dissection ordering
1039: - -mat_strumpack_metis_nodeNDP - Use METIS_NodeNDP instead of METIS_NodeND, for a more balanced tree
1041: Level: beginner
1043: Notes:
1044: Recommended use is 1 MPI process per GPU.
1046: Use `-pc_type lu` `-pc_factor_mat_solver_type strumpack` to use this as an exact (direct) solver.
1048: Use `-pc_type ilu` `-pc_factor_mat_solver_type strumpack` to enable low-rank compression (i.e, use as a preconditioner), by default using block low rank (BLR).
1050: Works with `MATAIJ` matrices
1052: HODLR, BLR_HODBF and ZFP_BLR_HODLR compression require STRUMPACK to be configured with ButterflyPACK support (`--download-butterflypack`).
1054: LOSSY, LOSSLESS and ZFP_BLR_HODLR compression require STRUMPACK to be configured with ZFP support (`--download-zfp`).
1056: .seealso: `MATSOLVERSTRUMPACK`, [](ch_matrices), `Mat`, `PCLU`, `PCILU`, `MATSOLVERSUPERLU_DIST`, `MATSOLVERMUMPS`, `PCFactorSetMatSolverType()`, `MatSolverType`,
1057: `MatGetFactor()`, `MatSTRUMPACKSetReordering()`, `MatSTRUMPACKReordering`, `MatSTRUMPACKCompressionType`, `MatSTRUMPACKSetColPerm()`.
1058: M*/
1059: static PetscErrorCode MatGetFactor_aij_strumpack(Mat A, MatFactorType ftype, Mat *F)
1060: {
1061: Mat B;
1062: PetscInt M = A->rmap->N, N = A->cmap->N;
1063: PetscBool verb, flg, set;
1064: PetscReal ctol;
1065: PetscInt min_sep_size, leaf_size, nxyz[3], nrdims, nc, w;
1066: #if defined(STRUMPACK_USE_ZFP)
1067: PetscInt lossy_prec;
1068: #endif
1069: #if defined(STRUMPACK_USE_BPACK)
1070: PetscInt bfly_lvls;
1071: #endif
1072: #if defined(STRUMPACK_USE_SLATE_SCALAPACK)
1073: PetscMPIInt mpithreads;
1074: #endif
1075: STRUMPACK_SparseSolver *S;
1076: STRUMPACK_INTERFACE iface;
1077: STRUMPACK_REORDERING_STRATEGY ndcurrent, ndvalue;
1078: STRUMPACK_COMPRESSION_TYPE compcurrent, compvalue;
1079: const STRUMPACK_PRECISION table[2][2][2] = {
1080: {{STRUMPACK_FLOATCOMPLEX_64, STRUMPACK_DOUBLECOMPLEX_64}, {STRUMPACK_FLOAT_64, STRUMPACK_DOUBLE_64}},
1081: {{STRUMPACK_FLOATCOMPLEX, STRUMPACK_DOUBLECOMPLEX}, {STRUMPACK_FLOAT, STRUMPACK_DOUBLE} }
1082: };
1083: const STRUMPACK_PRECISION prec = table[(sizeof(PetscInt) == 8) ? 0 : 1][(PETSC_SCALAR == PETSC_COMPLEX) ? 0 : 1][(PETSC_REAL == PETSC_FLOAT) ? 0 : 1];
1084: const char *const STRUMPACKNDTypes[] = {"NATURAL", "METIS", "PARMETIS", "SCOTCH", "PTSCOTCH", "RCM", "GEOMETRIC", "AMD", "MMD", "AND", "MLF", "SPECTRAL", "STRUMPACKNDTypes", "", 0};
1085: const char *const CompTypes[] = {"NONE", "HSS", "BLR", "HODLR", "BLR_HODLR", "ZFP_BLR_HODLR", "LOSSLESS", "LOSSY", "CompTypes", "", 0};
1087: PetscFunctionBegin;
1088: #if defined(STRUMPACK_USE_SLATE_SCALAPACK)
1089: PetscCallMPI(MPI_Query_thread(&mpithreads));
1090: PetscCheck(mpithreads == MPI_THREAD_MULTIPLE, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP_SYS, "SLATE requires MPI_THREAD_MULTIPLE");
1091: #endif
1092: /* Create the factorization matrix */
1093: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
1094: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, M, N));
1095: PetscCall(PetscStrallocpy("strumpack", &((PetscObject)B)->type_name));
1096: PetscCall(MatSetUp(B));
1097: PetscCall(MatSeqAIJSetPreallocation(B, 0, NULL));
1098: PetscCall(MatMPIAIJSetPreallocation(B, 0, NULL, 0, NULL));
1099: B->trivialsymbolic = PETSC_TRUE;
1100: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
1101: B->ops->lufactorsymbolic = MatLUFactorSymbolic_STRUMPACK;
1102: B->ops->ilufactorsymbolic = MatLUFactorSymbolic_STRUMPACK;
1103: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Factor type not supported");
1104: B->ops->getinfo = MatGetInfo_External;
1105: B->ops->view = MatView_STRUMPACK;
1106: B->ops->destroy = MatDestroy_STRUMPACK;
1107: B->ops->getdiagonal = MatGetDiagonal_STRUMPACK;
1108: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_aij_strumpack));
1109: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetReordering_C", MatSTRUMPACKSetReordering_STRUMPACK));
1110: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKGetReordering_C", MatSTRUMPACKGetReordering_STRUMPACK));
1111: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetColPerm_C", MatSTRUMPACKSetColPerm_STRUMPACK));
1112: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKGetColPerm_C", MatSTRUMPACKGetColPerm_STRUMPACK));
1113: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetGeometricNxyz_C", MatSTRUMPACKSetGeometricNxyz_STRUMPACK));
1114: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetGeometricComponents_C", MatSTRUMPACKSetGeometricComponents_STRUMPACK));
1115: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetGeometricWidth_C", MatSTRUMPACKSetGeometricWidth_STRUMPACK));
1116: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetGPU_C", MatSTRUMPACKSetGPU_STRUMPACK));
1117: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKGetGPU_C", MatSTRUMPACKGetGPU_STRUMPACK));
1118: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetCompression_C", MatSTRUMPACKSetCompression_STRUMPACK));
1119: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKGetCompression_C", MatSTRUMPACKGetCompression_STRUMPACK));
1120: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetCompRelTol_C", MatSTRUMPACKSetCompRelTol_STRUMPACK));
1121: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKGetCompRelTol_C", MatSTRUMPACKGetCompRelTol_STRUMPACK));
1122: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetCompAbsTol_C", MatSTRUMPACKSetCompAbsTol_STRUMPACK));
1123: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKGetCompAbsTol_C", MatSTRUMPACKGetCompAbsTol_STRUMPACK));
1124: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetCompLeafSize_C", MatSTRUMPACKSetCompLeafSize_STRUMPACK));
1125: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKGetCompLeafSize_C", MatSTRUMPACKGetCompLeafSize_STRUMPACK));
1126: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetCompMinSepSize_C", MatSTRUMPACKSetCompMinSepSize_STRUMPACK));
1127: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKGetCompMinSepSize_C", MatSTRUMPACKGetCompMinSepSize_STRUMPACK));
1128: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetCompLossyPrecision_C", MatSTRUMPACKSetCompLossyPrecision_STRUMPACK));
1129: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKGetCompLossyPrecision_C", MatSTRUMPACKGetCompLossyPrecision_STRUMPACK));
1130: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetCompButterflyLevels_C", MatSTRUMPACKSetCompButterflyLevels_STRUMPACK));
1131: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKGetCompButterflyLevels_C", MatSTRUMPACKGetCompButterflyLevels_STRUMPACK));
1132: B->factortype = ftype;
1134: /* set solvertype */
1135: PetscCall(PetscFree(B->solvertype));
1136: PetscCall(PetscStrallocpy(MATSOLVERSTRUMPACK, &B->solvertype));
1138: PetscCall(PetscNew(&S));
1139: B->data = S;
1141: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &flg)); /* A might be MATSEQAIJCUSPARSE */
1142: iface = flg ? STRUMPACK_MT : STRUMPACK_MPI_DIST;
1144: PetscOptionsBegin(PetscObjectComm((PetscObject)B), ((PetscObject)B)->prefix, "STRUMPACK Options", "Mat");
1145: verb = PetscLogPrintInfo ? PETSC_TRUE : PETSC_FALSE;
1146: PetscCall(PetscOptionsBool("-mat_strumpack_verbose", "Print STRUMPACK information", "None", verb, &verb, NULL));
1148: PetscStackCallExternalVoid("STRUMPACK_init", STRUMPACK_init(S, PetscObjectComm((PetscObject)A), prec, iface, 0, NULL, verb));
1150: /* By default, no compression is done. Compression is enabled when the user enables it with */
1151: /* -mat_strumpack_compression with anything else than NONE, or when selecting ilu */
1152: /* preconditioning, in which case we default to STRUMPACK_BLR compression. */
1153: /* When compression is enabled, the STRUMPACK solver becomes an incomplete */
1154: /* (or approximate) LU factorization. */
1155: PetscStackCallExternalVoid("STRUMPACK_compression", compcurrent = STRUMPACK_compression(*S));
1156: PetscCall(PetscOptionsEnum("-mat_strumpack_compression", "Rank-structured compression type", "None", CompTypes, (PetscEnum)compcurrent, (PetscEnum *)&compvalue, &set));
1157: if (set) {
1158: PetscCall(MatSTRUMPACKSetCompression(B, (MatSTRUMPACKCompressionType)compvalue));
1159: } else {
1160: if (ftype == MAT_FACTOR_ILU) { PetscStackCallExternalVoid("STRUMPACK_set_compression", STRUMPACK_set_compression(*S, STRUMPACK_BLR)); }
1161: }
1163: PetscStackCallExternalVoid("STRUMPACK_compression_rel_tol", ctol = (PetscReal)STRUMPACK_compression_rel_tol(*S));
1164: PetscCall(PetscOptionsReal("-mat_strumpack_compression_rel_tol", "Relative compression tolerance", "None", ctol, &ctol, &set));
1165: if (set) PetscStackCallExternalVoid("STRUMPACK_set_compression_rel_tol", STRUMPACK_set_compression_rel_tol(*S, (double)ctol));
1167: PetscStackCallExternalVoid("STRUMPACK_compression_abs_tol", ctol = (PetscReal)STRUMPACK_compression_abs_tol(*S));
1168: PetscCall(PetscOptionsReal("-mat_strumpack_compression_abs_tol", "Absolute compression tolerance", "None", ctol, &ctol, &set));
1169: if (set) PetscStackCallExternalVoid("STRUMPACK_set_compression_abs_tol", STRUMPACK_set_compression_abs_tol(*S, (double)ctol));
1171: PetscStackCallExternalVoid("STRUMPACK_compression_min_sep_size", min_sep_size = (PetscInt)STRUMPACK_compression_min_sep_size(*S));
1172: PetscCall(PetscOptionsInt("-mat_strumpack_compression_min_sep_size", "Minimum size of separator for compression", "None", min_sep_size, &min_sep_size, &set));
1173: if (set) PetscStackCallExternalVoid("STRUMPACK_set_compression_min_sep_size", STRUMPACK_set_compression_min_sep_size(*S, (int)min_sep_size));
1175: PetscStackCallExternalVoid("STRUMPACK_compression_leaf_size", leaf_size = (PetscInt)STRUMPACK_compression_leaf_size(*S));
1176: PetscCall(PetscOptionsInt("-mat_strumpack_compression_leaf_size", "Size of diagonal blocks in rank-structured approximation", "None", leaf_size, &leaf_size, &set));
1177: if (set) PetscStackCallExternalVoid("STRUMPACK_set_compression_leaf_size", STRUMPACK_set_compression_leaf_size(*S, (int)leaf_size));
1179: #if defined(STRUMPACK_USE_ZFP)
1180: PetscStackCallExternalVoid("STRUMPACK_compression_lossy_precision", lossy_prec = (PetscInt)STRUMPACK_compression_lossy_precision(*S));
1181: PetscCall(PetscOptionsInt("-mat_strumpack_compression_lossy_precision", "Number of bitplanes to use in lossy compression", "None", lossy_prec, &lossy_prec, &set));
1182: if (set) PetscStackCallExternalVoid("STRUMPACK_set_compression_lossy_precision", STRUMPACK_set_compression_lossy_precision(*S, (int)lossy_prec));
1183: #endif
1185: #if defined(STRUMPACK_USE_BPACK)
1186: PetscStackCallExternalVoid("STRUMPACK_compression_butterfly_levels", bfly_lvls = (PetscInt)STRUMPACK_compression_butterfly_levels(*S));
1187: PetscCall(PetscOptionsInt("-mat_strumpack_compression_butterfly_levels", "Number of levels in the HODLR matrix for which to use butterfly compression", "None", bfly_lvls, &bfly_lvls, &set));
1188: if (set) PetscStackCallExternalVoid("STRUMPACK_set_compression_butterfly_levels", STRUMPACK_set_compression_butterfly_levels(*S, (int)bfly_lvls));
1189: #endif
1191: #if defined(STRUMPACK_USE_CUDA) || defined(STRUMPACK_USE_HIP) || defined(STRUMPACK_USE_SYCL)
1192: PetscStackCallExternalVoid("STRUMPACK_use_gpu", flg = (STRUMPACK_use_gpu(*S) == 0) ? PETSC_FALSE : PETSC_TRUE);
1193: PetscCall(PetscOptionsBool("-mat_strumpack_gpu", "Enable GPU acceleration (not supported for all compression types)", "None", flg, &flg, &set));
1194: if (set) MatSTRUMPACKSetGPU(B, flg);
1195: #endif
1197: PetscStackCallExternalVoid("STRUMPACK_matching", flg = (STRUMPACK_matching(*S) == STRUMPACK_MATCHING_NONE) ? PETSC_FALSE : PETSC_TRUE);
1198: PetscCall(PetscOptionsBool("-mat_strumpack_colperm", "Find a col perm to get nonzero diagonal", "None", flg, &flg, &set));
1199: if (set) PetscStackCallExternalVoid("STRUMPACK_set_matching", STRUMPACK_set_matching(*S, flg ? STRUMPACK_MATCHING_MAX_DIAGONAL_PRODUCT_SCALING : STRUMPACK_MATCHING_NONE));
1201: PetscStackCallExternalVoid("STRUMPACK_reordering_method", ndcurrent = STRUMPACK_reordering_method(*S));
1202: PetscCall(PetscOptionsEnum("-mat_strumpack_reordering", "Sparsity reducing matrix reordering", "None", STRUMPACKNDTypes, (PetscEnum)ndcurrent, (PetscEnum *)&ndvalue, &set));
1203: if (set) PetscStackCallExternalVoid("STRUMPACK_set_reordering_method", STRUMPACK_set_reordering_method(*S, ndvalue));
1205: /* geometric ordering, for a regular 1D/2D/3D mesh in the natural ordering, */
1206: /* with nc DOF's per gridpoint, and possibly a wider stencil */
1207: nrdims = 3;
1208: nxyz[0] = nxyz[1] = nxyz[2] = 1;
1209: PetscCall(PetscOptionsIntArray("-mat_strumpack_geometric_xyz", "Mesh sizes nx,ny,nz (Use 1 for default)", "", nxyz, &nrdims, &set));
1210: if (set) {
1211: PetscCheck(nrdims == 1 || nrdims == 2 || nrdims == 3, PetscObjectComm((PetscObject)F), PETSC_ERR_ARG_OUTOFRANGE, "'-mat_strumpack_geometric_xyz' requires 1, 2, or 3 values.");
1212: PetscCall(MatSTRUMPACKSetGeometricNxyz(B, (int)nxyz[0], (int)nxyz[1], (int)nxyz[2]));
1213: }
1214: PetscCall(PetscOptionsInt("-mat_strumpack_geometric_components", "Number of components per mesh point, for geometric nested dissection ordering", "None", 1, &nc, &set));
1215: if (set) PetscStackCallExternalVoid("STRUMPACK_set_components", STRUMPACK_set_components(*S, (int)nc));
1216: PetscCall(PetscOptionsInt("-mat_strumpack_geometric_width", "Width of the separator (for instance a 1D 3-point wide stencil needs a 1 point wide separator, a 1D 5-point stencil needs a 2 point wide separator), for geometric nested dissection ordering", "None", 1, &w, &set));
1217: if (set) PetscStackCallExternalVoid("STRUMPACK_set_separator_width", STRUMPACK_set_separator_width(*S, (int)w));
1219: PetscStackCallExternalVoid("STRUMPACK_use_METIS_NodeNDP", flg = (STRUMPACK_use_METIS_NodeNDP(*S) == 0) ? PETSC_FALSE : PETSC_TRUE);
1220: PetscCall(PetscOptionsBool("-mat_strumpack_metis_nodeNDP", "Use METIS_NodeNDP instead of METIS_NodeND, for a more balanced tree", "None", flg, &flg, &set));
1221: if (set) {
1222: if (flg) {
1223: PetscStackCallExternalVoid("STRUMPACK_enable_METIS_NodeNDP", STRUMPACK_enable_METIS_NodeNDP(*S));
1224: } else {
1225: PetscStackCallExternalVoid("STRUMPACK_disable_METIS_NodeNDP", STRUMPACK_disable_METIS_NodeNDP(*S));
1226: }
1227: }
1229: /* Disable the outer iterative solver from STRUMPACK. */
1230: /* When STRUMPACK is used as a direct solver, it will by default do iterative refinement. */
1231: /* When STRUMPACK is used as an approximate factorization preconditioner (by enabling */
1232: /* low-rank compression), it will use it's own preconditioned GMRES. Here we can disable */
1233: /* the outer iterative solver, as PETSc uses STRUMPACK from within a KSP. */
1234: PetscStackCallExternalVoid("STRUMPACK_set_Krylov_solver", STRUMPACK_set_Krylov_solver(*S, STRUMPACK_DIRECT));
1236: PetscOptionsEnd();
1238: *F = B;
1239: PetscFunctionReturn(PETSC_SUCCESS);
1240: }
1242: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_STRUMPACK(void)
1243: {
1244: PetscFunctionBegin;
1245: PetscCall(MatSolverTypeRegister(MATSOLVERSTRUMPACK, MATMPIAIJ, MAT_FACTOR_LU, MatGetFactor_aij_strumpack));
1246: PetscCall(MatSolverTypeRegister(MATSOLVERSTRUMPACK, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_aij_strumpack));
1247: PetscCall(MatSolverTypeRegister(MATSOLVERSTRUMPACK, MATMPIAIJ, MAT_FACTOR_ILU, MatGetFactor_aij_strumpack));
1248: PetscCall(MatSolverTypeRegister(MATSOLVERSTRUMPACK, MATSEQAIJ, MAT_FACTOR_ILU, MatGetFactor_aij_strumpack));
1249: PetscFunctionReturn(PETSC_SUCCESS);
1250: }