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, "MatSTRUMPACKSetColPerm_C", NULL));
25: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetHSSRelTol_C", NULL));
26: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetHSSAbsTol_C", NULL));
27: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetHSSMaxRank_C", NULL));
28: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetHSSLeafSize_C", NULL));
29: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetHSSMinSepSize_C", NULL));
31: PetscFunctionReturn(PETSC_SUCCESS);
32: }
34: static PetscErrorCode MatSTRUMPACKSetReordering_STRUMPACK(Mat F, MatSTRUMPACKReordering reordering)
35: {
36: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
38: PetscFunctionBegin;
39: PetscStackCallExternalVoid("STRUMPACK_reordering_method", STRUMPACK_set_reordering_method(*S, (STRUMPACK_REORDERING_STRATEGY)reordering));
40: PetscFunctionReturn(PETSC_SUCCESS);
41: }
43: /*@
44: MatSTRUMPACKSetReordering - Set STRUMPACK fill-reducing reordering
46: Input Parameters:
47: + F - the factored matrix obtained by calling `MatGetFactor()`
48: - reordering - the code to be used to find the fill-reducing reordering, see `MatSTRUMPACKReordering`
50: Options Database Key:
51: . -mat_strumpack_reordering <METIS> - Sparsity reducing matrix reordering, see `MatSTRUMPACKReordering`
53: Level: beginner
55: References:
56: . * - STRUMPACK manual
58: .seealso: [](ch_matrices), `Mat`, `MatSTRUMPACKReordering`, `MatGetFactor()`, `MatSTRUMPACKSetColPerm()`, `MatSTRUMPACKSetHSSRelTol()`,
59: `MatSTRUMPACKSetHSSAbsTol()`, `MatSTRUMPACKSetHSSMaxRank()`, `MatSTRUMPACKSetHSSLeafSize()`, `MatSTRUMPACKSetHSSMinSepSize()`
60: @*/
61: PetscErrorCode MatSTRUMPACKSetReordering(Mat F, MatSTRUMPACKReordering reordering)
62: {
63: PetscFunctionBegin;
66: PetscTryMethod(F, "MatSTRUMPACKSetReordering_C", (Mat, MatSTRUMPACKReordering), (F, reordering));
67: PetscFunctionReturn(PETSC_SUCCESS);
68: }
70: static PetscErrorCode MatSTRUMPACKSetColPerm_STRUMPACK(Mat F, PetscBool cperm)
71: {
72: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
74: PetscFunctionBegin;
75: PetscStackCallExternalVoid("STRUMPACK_set_mc64job", STRUMPACK_set_mc64job(*S, cperm ? 5 : 0));
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: /*@
80: MatSTRUMPACKSetColPerm - Set whether STRUMPACK should try to permute the columns of the matrix in order to get a nonzero diagonal
82: Logically Collective
84: Input Parameters:
85: + F - the factored matrix obtained by calling `MatGetFactor()`
86: - cperm - `PETSC_TRUE` to permute (internally) the columns of the matrix
88: Options Database Key:
89: . -mat_strumpack_colperm <cperm> - true to use the permutation
91: Level: beginner
93: References:
94: . * - STRUMPACK manual
96: .seealso: [](ch_matrices), `MatSTRUMPACKSetReordering()`, `Mat`, `MatGetFactor()`, `MatSTRUMPACKSetHSSRelTol()`, `MatSTRUMPACKSetHSSAbsTol()`,
97: `MatSTRUMPACKSetHSSMaxRank()`, `MatSTRUMPACKSetHSSLeafSize()`, `MatSTRUMPACKSetHSSMinSepSize()`
98: @*/
99: PetscErrorCode MatSTRUMPACKSetColPerm(Mat F, PetscBool cperm)
100: {
101: PetscFunctionBegin;
104: PetscTryMethod(F, "MatSTRUMPACKSetColPerm_C", (Mat, PetscBool), (F, cperm));
105: PetscFunctionReturn(PETSC_SUCCESS);
106: }
108: static PetscErrorCode MatSTRUMPACKSetHSSRelTol_STRUMPACK(Mat F, PetscReal rtol)
109: {
110: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
112: PetscFunctionBegin;
113: PetscStackCallExternalVoid("STRUMPACK_set_HSS_rel_tol", STRUMPACK_set_HSS_rel_tol(*S, rtol));
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: /*@
118: MatSTRUMPACKSetHSSRelTol - Set STRUMPACK relative tolerance for HSS compression
120: Logically Collective
122: Input Parameters:
123: + F - the factored matrix obtained by calling `MatGetFactor()`
124: - rtol - relative compression tolerance
126: Options Database Key:
127: . -mat_strumpack_hss_rel_tol <1e-2> - Relative compression tolerance (None)
129: Level: beginner
131: References:
132: . * - STRUMPACK manual
134: .seealso: [](ch_matrices), `Mat`, `MatGetFactor()`, `MatSTRUMPACKSetReordering()`, `MatSTRUMPACKSetColPerm()`, `MatSTRUMPACKSetHSSAbsTol()`,
135: `MatSTRUMPACKSetHSSMaxRank()`, `MatSTRUMPACKSetHSSLeafSize()`, `MatSTRUMPACKSetHSSMinSepSize()`
136: @*/
137: PetscErrorCode MatSTRUMPACKSetHSSRelTol(Mat F, PetscReal rtol)
138: {
139: PetscFunctionBegin;
142: PetscTryMethod(F, "MatSTRUMPACKSetHSSRelTol_C", (Mat, PetscReal), (F, rtol));
143: PetscFunctionReturn(PETSC_SUCCESS);
144: }
146: static PetscErrorCode MatSTRUMPACKSetHSSAbsTol_STRUMPACK(Mat F, PetscReal atol)
147: {
148: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
150: PetscFunctionBegin;
151: PetscStackCallExternalVoid("STRUMPACK_set_HSS_abs_tol", STRUMPACK_set_HSS_abs_tol(*S, atol));
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
155: /*@
156: MatSTRUMPACKSetHSSAbsTol - Set STRUMPACK absolute tolerance for HSS compression
158: Logically Collective
160: Input Parameters:
161: + F - the factored matrix obtained by calling `MatGetFactor()`
162: - atol - absolute compression tolerance
164: Options Database Key:
165: . -mat_strumpack_hss_abs_tol <1e-8> - Absolute compression tolerance (None)
167: Level: beginner
169: References:
170: . * - STRUMPACK manual
172: .seealso: [](ch_matrices), `Mat`, `MatGetFactor()`, `MatSTRUMPACKSetReordering()`, `MatSTRUMPACKSetColPerm()`, `MatSTRUMPACKSetHSSRelTol()`,
173: `MatSTRUMPACKSetHSSMaxRank()`, `MatSTRUMPACKSetHSSLeafSize()`, `MatSTRUMPACKSetHSSMinSepSize()`
174: @*/
175: PetscErrorCode MatSTRUMPACKSetHSSAbsTol(Mat F, PetscReal atol)
176: {
177: PetscFunctionBegin;
180: PetscTryMethod(F, "MatSTRUMPACKSetHSSAbsTol_C", (Mat, PetscReal), (F, atol));
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: static PetscErrorCode MatSTRUMPACKSetHSSMaxRank_STRUMPACK(Mat F, PetscInt hssmaxrank)
185: {
186: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
188: PetscFunctionBegin;
189: PetscStackCallExternalVoid("STRUMPACK_set_HSS_max_rank", STRUMPACK_set_HSS_max_rank(*S, hssmaxrank));
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }
193: /*@
194: MatSTRUMPACKSetHSSMaxRank - Set STRUMPACK maximum HSS rank
196: Logically Collective
198: Input Parameters:
199: + F - the factored matrix obtained by calling `MatGetFactor()`
200: - hssmaxrank - maximum rank used in low-rank approximation
202: Options Database Key:
203: . -mat_strumpack_max_rank - Maximum rank in HSS compression, when using pctype ilu (None)
205: Level: beginner
207: References:
208: . * - STRUMPACK manual
210: .seealso: [](ch_matrices), `Mat`, `MatGetFactor()`, `MatSTRUMPACKSetReordering()`, `MatSTRUMPACKSetColPerm()`, `MatSTRUMPACKSetHSSRelTol()`,
211: `MatSTRUMPACKSetHSSAbsTol()`, `MatSTRUMPACKSetHSSLeafSize()`, `MatSTRUMPACKSetHSSMinSepSize()`
212: @*/
213: PetscErrorCode MatSTRUMPACKSetHSSMaxRank(Mat F, PetscInt hssmaxrank)
214: {
215: PetscFunctionBegin;
218: PetscTryMethod(F, "MatSTRUMPACKSetHSSMaxRank_C", (Mat, PetscInt), (F, hssmaxrank));
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: static PetscErrorCode MatSTRUMPACKSetHSSLeafSize_STRUMPACK(Mat F, PetscInt leaf_size)
223: {
224: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
226: PetscFunctionBegin;
227: PetscStackCallExternalVoid("STRUMPACK_set_HSS_leaf_size", STRUMPACK_set_HSS_leaf_size(*S, leaf_size));
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: /*@
232: MatSTRUMPACKSetHSSLeafSize - Set STRUMPACK HSS leaf size
234: Logically Collective
236: Input Parameters:
237: + F - the factored matrix obtained by calling `MatGetFactor()`
238: - leaf_size - Size of diagonal blocks in HSS approximation
240: Options Database Key:
241: . -mat_strumpack_leaf_size - Size of diagonal blocks in HSS approximation, when using pctype ilu (None)
243: Level: beginner
245: References:
246: . * - STRUMPACK manual
248: .seealso: [](ch_matrices), `Mat`, `MatGetFactor()`, `MatSTRUMPACKSetReordering()`, `MatSTRUMPACKSetColPerm()`, `MatSTRUMPACKSetHSSRelTol()`,
249: `MatSTRUMPACKSetHSSAbsTol()`, `MatSTRUMPACKSetHSSMaxRank()`, `MatSTRUMPACKSetHSSMinSepSize()`
250: @*/
251: PetscErrorCode MatSTRUMPACKSetHSSLeafSize(Mat F, PetscInt leaf_size)
252: {
253: PetscFunctionBegin;
256: PetscTryMethod(F, "MatSTRUMPACKSetHSSLeafSize_C", (Mat, PetscInt), (F, leaf_size));
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }
260: static PetscErrorCode MatSTRUMPACKSetHSSMinSepSize_STRUMPACK(Mat F, PetscInt hssminsize)
261: {
262: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
264: PetscFunctionBegin;
265: PetscStackCallExternalVoid("STRUMPACK_set_HSS_min_sep_size", STRUMPACK_set_HSS_min_sep_size(*S, hssminsize));
266: PetscFunctionReturn(PETSC_SUCCESS);
267: }
269: /*@
270: MatSTRUMPACKSetHSSMinSepSize - Set STRUMPACK minimum separator size for low-rank approximation
272: Logically Collective
274: Input Parameters:
275: + F - the factored matrix obtained by calling `MatGetFactor()`
276: - hssminsize - minimum dense matrix size for low-rank approximation
278: Options Database Key:
279: . -mat_strumpack_hss_min_sep_size <hssminsize> - set the minimum separator size
281: Level: beginner
283: References:
284: . * - STRUMPACK manual
286: .seealso: [](ch_matrices), `Mat`, `MatGetFactor()`, `MatSTRUMPACKSetReordering()`, `MatSTRUMPACKSetColPerm()`, `MatSTRUMPACKSetHSSRelTol()`,
287: `MatSTRUMPACKSetHSSAbsTol()`, `MatSTRUMPACKSetHSSMaxRank()`, `MatSTRUMPACKSetHSSLeafSize()`
288: @*/
289: PetscErrorCode MatSTRUMPACKSetHSSMinSepSize(Mat F, PetscInt hssminsize)
290: {
291: PetscFunctionBegin;
294: PetscTryMethod(F, "MatSTRUMPACKSetHSSMinSepSize_C", (Mat, PetscInt), (F, hssminsize));
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }
298: static PetscErrorCode MatSolve_STRUMPACK(Mat A, Vec b_mpi, Vec x)
299: {
300: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)A->data;
301: STRUMPACK_RETURN_CODE sp_err;
302: const PetscScalar *bptr;
303: PetscScalar *xptr;
305: PetscFunctionBegin;
306: PetscCall(VecGetArray(x, &xptr));
307: PetscCall(VecGetArrayRead(b_mpi, &bptr));
309: PetscStackCallExternalVoid("STRUMPACK_solve", sp_err = STRUMPACK_solve(*S, (PetscScalar *)bptr, xptr, 0));
310: switch (sp_err) {
311: case STRUMPACK_SUCCESS:
312: break;
313: case STRUMPACK_MATRIX_NOT_SET: {
314: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: matrix was not set");
315: break;
316: }
317: case STRUMPACK_REORDERING_ERROR: {
318: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: matrix reordering failed");
319: break;
320: }
321: default:
322: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: solve failed");
323: }
324: PetscCall(VecRestoreArray(x, &xptr));
325: PetscCall(VecRestoreArrayRead(b_mpi, &bptr));
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }
329: static PetscErrorCode MatMatSolve_STRUMPACK(Mat A, Mat B_mpi, Mat X)
330: {
331: PetscBool flg;
333: PetscFunctionBegin;
334: PetscCall(PetscObjectTypeCompareAny((PetscObject)B_mpi, &flg, MATSEQDENSE, MATMPIDENSE, NULL));
335: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix B must be MATDENSE matrix");
336: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &flg, MATSEQDENSE, MATMPIDENSE, NULL));
337: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix X must be MATDENSE matrix");
338: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatMatSolve_STRUMPACK() is not implemented yet");
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }
342: static PetscErrorCode MatView_Info_STRUMPACK(Mat A, PetscViewer viewer)
343: {
344: PetscFunctionBegin;
345: /* check if matrix is strumpack type */
346: if (A->ops->solve != MatSolve_STRUMPACK) PetscFunctionReturn(PETSC_SUCCESS);
347: PetscCall(PetscViewerASCIIPrintf(viewer, "STRUMPACK sparse solver!\n"));
348: PetscFunctionReturn(PETSC_SUCCESS);
349: }
351: static PetscErrorCode MatView_STRUMPACK(Mat A, PetscViewer viewer)
352: {
353: PetscBool iascii;
354: PetscViewerFormat format;
356: PetscFunctionBegin;
357: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
358: if (iascii) {
359: PetscCall(PetscViewerGetFormat(viewer, &format));
360: if (format == PETSC_VIEWER_ASCII_INFO) PetscCall(MatView_Info_STRUMPACK(A, viewer));
361: }
362: PetscFunctionReturn(PETSC_SUCCESS);
363: }
365: static PetscErrorCode MatLUFactorNumeric_STRUMPACK(Mat F, Mat A, const MatFactorInfo *info)
366: {
367: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
368: STRUMPACK_RETURN_CODE sp_err;
369: Mat_SeqAIJ *A_d, *A_o;
370: Mat_MPIAIJ *mat;
371: PetscInt M = A->rmap->N, m = A->rmap->n;
372: PetscBool flg;
373: const PetscScalar *A_d_a, *A_o_a;
375: PetscFunctionBegin;
376: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &flg));
377: if (flg) { /* A might be MATMPIAIJCUSPARSE etc */
378: mat = (Mat_MPIAIJ *)A->data;
379: PetscCall(MatSeqAIJGetArrayRead(mat->A, &A_d_a)); /* Make sure mat is sync'ed to host */
380: PetscCall(MatSeqAIJGetArrayRead(mat->B, &A_o_a));
381: A_d = (Mat_SeqAIJ *)(mat->A)->data;
382: A_o = (Mat_SeqAIJ *)(mat->B)->data;
383: PetscStackCallExternalVoid("STRUMPACK_set_MPIAIJ_matrix", STRUMPACK_set_MPIAIJ_matrix(*S, &m, A_d->i, A_d->j, A_d_a, A_o->i, A_o->j, A_o_a, mat->garray));
384: PetscCall(MatSeqAIJRestoreArrayRead(mat->A, &A_d_a));
385: PetscCall(MatSeqAIJRestoreArrayRead(mat->B, &A_o_a));
386: } else { /* A might be MATSEQAIJCUSPARSE etc */
387: PetscCall(MatSeqAIJGetArrayRead(A, &A_d_a));
388: A_d = (Mat_SeqAIJ *)A->data;
389: PetscStackCallExternalVoid("STRUMPACK_set_csr_matrix", STRUMPACK_set_csr_matrix(*S, &M, A_d->i, A_d->j, A_d_a, 0));
390: PetscCall(MatSeqAIJRestoreArrayRead(A, &A_d_a));
391: }
393: /* Reorder and Factor the matrix. */
394: /* TODO figure out how to avoid reorder if the matrix values changed, but the pattern remains the same. */
395: PetscStackCallExternalVoid("STRUMPACK_reorder", sp_err = STRUMPACK_reorder(*S));
396: PetscStackCallExternalVoid("STRUMPACK_factor", sp_err = STRUMPACK_factor(*S));
397: switch (sp_err) {
398: case STRUMPACK_SUCCESS:
399: break;
400: case STRUMPACK_MATRIX_NOT_SET: {
401: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: matrix was not set");
402: break;
403: }
404: case STRUMPACK_REORDERING_ERROR: {
405: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: matrix reordering failed");
406: break;
407: }
408: default:
409: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: factorization failed");
410: }
411: F->assembled = PETSC_TRUE;
412: F->preallocated = PETSC_TRUE;
413: PetscFunctionReturn(PETSC_SUCCESS);
414: }
416: static PetscErrorCode MatLUFactorSymbolic_STRUMPACK(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
417: {
418: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
419: PetscBool flg, set;
420: PetscReal ctol;
421: PetscInt hssminsize, max_rank, leaf_size;
422: STRUMPACK_REORDERING_STRATEGY ndcurrent, ndvalue;
423: STRUMPACK_KRYLOV_SOLVER itcurrent, itsolver;
424: const char *const STRUMPACKNDTypes[] = {"NATURAL", "METIS", "PARMETIS", "SCOTCH", "PTSCOTCH", "RCM", "STRUMPACKNDTypes", "", 0};
425: const char *const SolverTypes[] = {"AUTO", "NONE", "REFINE", "PREC_GMRES", "GMRES", "PREC_BICGSTAB", "BICGSTAB", "SolverTypes", "", 0};
427: PetscFunctionBegin;
428: /* Set options to F */
429: PetscOptionsBegin(PetscObjectComm((PetscObject)F), ((PetscObject)F)->prefix, "STRUMPACK Options", "Mat");
430: PetscStackCallExternalVoid("STRUMPACK_HSS_rel_tol", ctol = (PetscReal)STRUMPACK_HSS_rel_tol(*S));
431: PetscCall(PetscOptionsReal("-mat_strumpack_hss_rel_tol", "Relative compression tolerance", "None", ctol, &ctol, &set));
432: if (set) PetscStackCallExternalVoid("STRUMPACK_set_HSS_rel_tol", STRUMPACK_set_HSS_rel_tol(*S, (double)ctol));
434: PetscStackCallExternalVoid("STRUMPACK_HSS_abs_tol", ctol = (PetscReal)STRUMPACK_HSS_abs_tol(*S));
435: PetscCall(PetscOptionsReal("-mat_strumpack_hss_abs_tol", "Absolute compression tolerance", "None", ctol, &ctol, &set));
436: if (set) PetscStackCallExternalVoid("STRUMPACK_set_HSS_abs_tol", STRUMPACK_set_HSS_abs_tol(*S, (double)ctol));
438: PetscStackCallExternalVoid("STRUMPACK_mc64job", flg = (STRUMPACK_mc64job(*S) == 0) ? PETSC_FALSE : PETSC_TRUE);
439: PetscCall(PetscOptionsBool("-mat_strumpack_colperm", "Find a col perm to get nonzero diagonal", "None", flg, &flg, &set));
440: if (set) PetscStackCallExternalVoid("STRUMPACK_set_mc64job", STRUMPACK_set_mc64job(*S, flg ? 5 : 0));
442: PetscStackCallExternalVoid("STRUMPACK_HSS_min_sep_size", hssminsize = (PetscInt)STRUMPACK_HSS_min_sep_size(*S));
443: PetscCall(PetscOptionsInt("-mat_strumpack_hss_min_sep_size", "Minimum size of separator for HSS compression", "None", hssminsize, &hssminsize, &set));
444: if (set) PetscStackCallExternalVoid("STRUMPACK_set_HSS_min_sep_size", STRUMPACK_set_HSS_min_sep_size(*S, (int)hssminsize));
446: PetscStackCallExternalVoid("STRUMPACK_HSS_max_rank", max_rank = (PetscInt)STRUMPACK_HSS_max_rank(*S));
447: PetscCall(PetscOptionsInt("-mat_strumpack_max_rank", "Maximum rank in HSS approximation", "None", max_rank, &max_rank, &set));
448: if (set) PetscStackCallExternalVoid("STRUMPACK_set_HSS_max_rank", STRUMPACK_set_HSS_max_rank(*S, (int)max_rank));
450: PetscStackCallExternalVoid("STRUMPACK_HSS_leaf_size", leaf_size = (PetscInt)STRUMPACK_HSS_leaf_size(*S));
451: PetscCall(PetscOptionsInt("-mat_strumpack_leaf_size", "Size of diagonal blocks in HSS approximation", "None", leaf_size, &leaf_size, &set));
452: if (set) PetscStackCallExternalVoid("STRUMPACK_set_HSS_leaf_size", STRUMPACK_set_HSS_leaf_size(*S, (int)leaf_size));
454: PetscStackCallExternalVoid("STRUMPACK_reordering_method", ndcurrent = STRUMPACK_reordering_method(*S));
455: PetscCall(PetscOptionsEnum("-mat_strumpack_reordering", "Sparsity reducing matrix reordering", "None", STRUMPACKNDTypes, (PetscEnum)ndcurrent, (PetscEnum *)&ndvalue, &set));
456: if (set) PetscStackCallExternalVoid("STRUMPACK_set_reordering_method", STRUMPACK_set_reordering_method(*S, ndvalue));
458: /* Disable the outer iterative solver from STRUMPACK. */
459: /* When STRUMPACK is used as a direct solver, it will by default do iterative refinement. */
460: /* When STRUMPACK is used as an approximate factorization preconditioner (by enabling */
461: /* low-rank compression), it will use it's own preconditioned GMRES. Here we can disable */
462: /* the outer iterative solver, as PETSc uses STRUMPACK from within a KSP. */
463: PetscStackCallExternalVoid("STRUMPACK_set_Krylov_solver", STRUMPACK_set_Krylov_solver(*S, STRUMPACK_DIRECT));
465: PetscStackCallExternalVoid("STRUMPACK_Krylov_solver", itcurrent = STRUMPACK_Krylov_solver(*S));
466: PetscCall(PetscOptionsEnum("-mat_strumpack_iterative_solver", "Select iterative solver from STRUMPACK", "None", SolverTypes, (PetscEnum)itcurrent, (PetscEnum *)&itsolver, &set));
467: if (set) PetscStackCallExternalVoid("STRUMPACK_set_Krylov_solver", STRUMPACK_set_Krylov_solver(*S, itsolver));
468: PetscOptionsEnd();
470: F->ops->lufactornumeric = MatLUFactorNumeric_STRUMPACK;
471: F->ops->solve = MatSolve_STRUMPACK;
472: F->ops->matsolve = MatMatSolve_STRUMPACK;
473: PetscFunctionReturn(PETSC_SUCCESS);
474: }
476: static PetscErrorCode MatFactorGetSolverType_aij_strumpack(Mat A, MatSolverType *type)
477: {
478: PetscFunctionBegin;
479: *type = MATSOLVERSTRUMPACK;
480: PetscFunctionReturn(PETSC_SUCCESS);
481: }
483: /*MC
484: MATSOLVERSSTRUMPACK = "strumpack" - A solver package providing a direct sparse solver (PCLU)
485: and a preconditioner (PCILU) using low-rank compression via the external package STRUMPACK.
487: Consult the STRUMPACK-sparse manual for more info.
489: Use ` ./configure --download-strumpack` to have PETSc installed with STRUMPACK
491: Use `-pc_type lu` `-pc_factor_mat_solver_type strumpack` to use this as an exact (direct) solver.
493: Use `-pc_type ilu` `-pc_factor_mat_solver_type strumpack` to enable low-rank compression (i.e, use as a preconditioner).
495: Works with `MATAIJ` matrices
497: Options Database Keys:
498: + -mat_strumpack_verbose - verbose info
499: . -mat_strumpack_hss_rel_tol <1e-2> - Relative compression tolerance (None)
500: . -mat_strumpack_hss_abs_tol <1e-8> - Absolute compression tolerance (None)
501: . -mat_strumpack_colperm <TRUE> - Permute matrix to make diagonal nonzeros (None)
502: . -mat_strumpack_hss_min_sep_size <256> - Minimum size of separator for HSS compression (None)
503: . -mat_strumpack_max_rank - Maximum rank in HSS compression, when using pctype ilu (None)
504: . -mat_strumpack_leaf_size - Size of diagonal blocks in HSS approximation, when using pctype ilu (None)
505: . -mat_strumpack_reordering <METIS> - Sparsity reducing matrix reordering see `MatSTRUMPACKReordering`
506: - -mat_strumpack_iterative_solver <DIRECT> - Select iterative solver from STRUMPACK (choose one of) `AUTO`, `DIRECT`, `REFINE`, `PREC_GMRES`,
507: `GMRES`, `PREC_BICGSTAB`, `BICGSTAB`
509: Level: beginner
511: .seealso: [](ch_matrices), `Mat`, `PCLU`, `PCILU`, `MATSOLVERSUPERLU_DIST`, `MATSOLVERMUMPS`, `PCFactorSetMatSolverType()`, `MatSolverType`,
512: `MatGetFactor()`, `MatSTRUMPACKSetReordering()`, `MatSTRUMPACKSetColPerm()`, `MatSTRUMPACKSetHSSRelTol()`,
513: `MatSTRUMPACKSetHSSAbsTol()`, `MatSTRUMPACKSetHSSMaxRank()`, `MatSTRUMPACKSetHSSLeafSize()`, `MatSTRUMPACKSetHSSMinSepSize()`
514: M*/
515: static PetscErrorCode MatGetFactor_aij_strumpack(Mat A, MatFactorType ftype, Mat *F)
516: {
517: Mat B;
518: PetscInt M = A->rmap->N, N = A->cmap->N;
519: PetscBool verb, flg;
520: STRUMPACK_SparseSolver *S;
521: STRUMPACK_INTERFACE iface;
522: const STRUMPACK_PRECISION table[2][2][2] = {
523: {{STRUMPACK_FLOATCOMPLEX_64, STRUMPACK_DOUBLECOMPLEX_64}, {STRUMPACK_FLOAT_64, STRUMPACK_DOUBLE_64}},
524: {{STRUMPACK_FLOATCOMPLEX, STRUMPACK_DOUBLECOMPLEX}, {STRUMPACK_FLOAT, STRUMPACK_DOUBLE} }
525: };
526: const STRUMPACK_PRECISION prec = table[(sizeof(PetscInt) == 8) ? 0 : 1][(PETSC_SCALAR == PETSC_COMPLEX) ? 0 : 1][(PETSC_REAL == PETSC_FLOAT) ? 0 : 1];
528: PetscFunctionBegin;
529: /* Create the factorization matrix */
530: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
531: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, M, N));
532: PetscCall(PetscStrallocpy("strumpack", &((PetscObject)B)->type_name));
533: PetscCall(MatSetUp(B));
534: B->trivialsymbolic = PETSC_TRUE;
535: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
536: B->ops->lufactorsymbolic = MatLUFactorSymbolic_STRUMPACK;
537: B->ops->ilufactorsymbolic = MatLUFactorSymbolic_STRUMPACK;
538: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Factor type not supported");
539: B->ops->getinfo = MatGetInfo_External;
540: B->ops->view = MatView_STRUMPACK;
541: B->ops->destroy = MatDestroy_STRUMPACK;
542: B->ops->getdiagonal = MatGetDiagonal_STRUMPACK;
543: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_aij_strumpack));
544: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetReordering_C", MatSTRUMPACKSetReordering_STRUMPACK));
545: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetColPerm_C", MatSTRUMPACKSetColPerm_STRUMPACK));
546: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetHSSRelTol_C", MatSTRUMPACKSetHSSRelTol_STRUMPACK));
547: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetHSSAbsTol_C", MatSTRUMPACKSetHSSAbsTol_STRUMPACK));
548: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetHSSMaxRank_C", MatSTRUMPACKSetHSSMaxRank_STRUMPACK));
549: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetHSSLeafSize_C", MatSTRUMPACKSetHSSLeafSize_STRUMPACK));
550: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetHSSMinSepSize_C", MatSTRUMPACKSetHSSMinSepSize_STRUMPACK));
551: B->factortype = ftype;
552: PetscCall(PetscNew(&S));
553: B->data = S;
555: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &flg)); /* A might be MATSEQAIJCUSPARSE */
556: iface = flg ? STRUMPACK_MT : STRUMPACK_MPI_DIST;
558: PetscOptionsBegin(PetscObjectComm((PetscObject)B), ((PetscObject)B)->prefix, "STRUMPACK Options", "Mat");
559: verb = PetscLogPrintInfo ? PETSC_TRUE : PETSC_FALSE;
560: PetscCall(PetscOptionsBool("-mat_strumpack_verbose", "Print STRUMPACK information", "None", verb, &verb, NULL));
562: PetscStackCallExternalVoid("STRUMPACK_init", STRUMPACK_init(S, PetscObjectComm((PetscObject)A), prec, iface, 0, NULL, verb));
564: if (ftype == MAT_FACTOR_ILU) {
565: /* When enabling HSS compression, the STRUMPACK solver becomes an incomplete */
566: /* (or approximate) LU factorization. */
567: PetscStackCallExternalVoid("STRUMPACK_enable_HSS", STRUMPACK_enable_HSS(*S));
568: }
569: PetscOptionsEnd();
571: /* set solvertype */
572: PetscCall(PetscFree(B->solvertype));
573: PetscCall(PetscStrallocpy(MATSOLVERSTRUMPACK, &B->solvertype));
575: *F = B;
576: PetscFunctionReturn(PETSC_SUCCESS);
577: }
579: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_STRUMPACK(void)
580: {
581: PetscFunctionBegin;
582: PetscCall(MatSolverTypeRegister(MATSOLVERSTRUMPACK, MATMPIAIJ, MAT_FACTOR_LU, MatGetFactor_aij_strumpack));
583: PetscCall(MatSolverTypeRegister(MATSOLVERSTRUMPACK, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_aij_strumpack));
584: PetscCall(MatSolverTypeRegister(MATSOLVERSTRUMPACK, MATMPIAIJ, MAT_FACTOR_ILU, MatGetFactor_aij_strumpack));
585: PetscCall(MatSolverTypeRegister(MATSOLVERSTRUMPACK, MATSEQAIJ, MAT_FACTOR_ILU, MatGetFactor_aij_strumpack));
586: PetscFunctionReturn(PETSC_SUCCESS);
587: }