Actual source code: istrumpack.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: PetscCallExternalVoid("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, "MatSTRUMPACKSetCompLeafSize_C", NULL));
39: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKGetCompLeafSize_C", NULL));
40: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetCompMinSepSize_C", NULL));
41: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKGetCompMinSepSize_C", NULL));
42: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetCompLossyPrecision_C", NULL));
43: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKGetCompLossyPrecision_C", NULL));
44: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKSetCompButterflyLevels_C", NULL));
45: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSTRUMPACKGetCompButterflyLevels_C", NULL));
46: PetscFunctionReturn(PETSC_SUCCESS);
47: }
49: static PetscErrorCode MatSTRUMPACKSetReordering_STRUMPACK(Mat F, MatSTRUMPACKReordering reordering)
50: {
51: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
53: PetscFunctionBegin;
54: PetscCallExternalVoid("STRUMPACK_set_reordering_method", STRUMPACK_set_reordering_method(*S, (STRUMPACK_REORDERING_STRATEGY)reordering));
55: PetscFunctionReturn(PETSC_SUCCESS);
56: }
57: static PetscErrorCode MatSTRUMPACKGetReordering_STRUMPACK(Mat F, MatSTRUMPACKReordering *reordering)
58: {
59: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
61: PetscFunctionBegin;
62: PetscCallExternalVoid("STRUMPACK_reordering_method", *reordering = (MatSTRUMPACKReordering)STRUMPACK_reordering_method(*S));
63: PetscFunctionReturn(PETSC_SUCCESS);
64: }
66: static PetscErrorCode MatSTRUMPACKSetColPerm_STRUMPACK(Mat F, PetscBool cperm)
67: {
68: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
70: PetscFunctionBegin;
71: PetscCallExternalVoid("STRUMPACK_set_matching", STRUMPACK_set_matching(*S, cperm ? STRUMPACK_MATCHING_MAX_DIAGONAL_PRODUCT_SCALING : STRUMPACK_MATCHING_NONE));
72: PetscFunctionReturn(PETSC_SUCCESS);
73: }
74: static PetscErrorCode MatSTRUMPACKGetColPerm_STRUMPACK(Mat F, PetscBool *cperm)
75: {
76: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
78: PetscFunctionBegin;
79: PetscCallExternalVoid("STRUMPACK_matching", *cperm = (PetscBool)(STRUMPACK_matching(*S) != STRUMPACK_MATCHING_NONE));
80: PetscFunctionReturn(PETSC_SUCCESS);
81: }
83: static PetscErrorCode MatSTRUMPACKSetGPU_STRUMPACK(Mat F, PetscBool gpu)
84: {
85: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
87: PetscFunctionBegin;
88: if (gpu) {
89: PetscCheck(PetscDefined(STRUMPACK_USE_CUDA) || PetscDefined(STRUMPACK_USE_HIP), PetscObjectComm((PetscObject)F), PETSC_ERR_SUP, "Strumpack was not configured with GPU support");
90: PetscCallExternalVoid("STRUMPACK_enable_gpu", STRUMPACK_enable_gpu(*S));
91: } else PetscCallExternalVoid("STRUMPACK_disable_gpu", STRUMPACK_disable_gpu(*S));
92: PetscFunctionReturn(PETSC_SUCCESS);
93: }
94: static PetscErrorCode MatSTRUMPACKGetGPU_STRUMPACK(Mat F, PetscBool *gpu)
95: {
96: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
98: PetscFunctionBegin;
99: PetscCallExternalVoid("STRUMPACK_use_gpu", *gpu = (PetscBool)STRUMPACK_use_gpu(*S));
100: PetscFunctionReturn(PETSC_SUCCESS);
101: }
103: static PetscErrorCode MatSTRUMPACKSetCompression_STRUMPACK(Mat F, MatSTRUMPACKCompressionType comp)
104: {
105: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
107: PetscFunctionBegin;
108: #if !defined(STRUMPACK_USE_BPACK)
109: 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");
110: #endif
111: #if !defined(STRUMPACK_USE_ZFP)
112: 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");
113: #endif
114: PetscCallExternalVoid("STRUMPACK_set_compression", STRUMPACK_set_compression(*S, (STRUMPACK_COMPRESSION_TYPE)comp));
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
117: static PetscErrorCode MatSTRUMPACKGetCompression_STRUMPACK(Mat F, MatSTRUMPACKCompressionType *comp)
118: {
119: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
121: PetscFunctionBegin;
122: PetscCallExternalVoid("STRUMPACK_compression", *comp = (MatSTRUMPACKCompressionType)STRUMPACK_compression(*S));
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: static PetscErrorCode MatSTRUMPACKSetCompRelTol_STRUMPACK(Mat F, PetscReal rtol)
127: {
128: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
130: PetscFunctionBegin;
131: PetscCallExternalVoid("STRUMPACK_set_compression_rel_tol", STRUMPACK_set_compression_rel_tol(*S, rtol));
132: PetscFunctionReturn(PETSC_SUCCESS);
133: }
134: static PetscErrorCode MatSTRUMPACKGetCompRelTol_STRUMPACK(Mat F, PetscReal *rtol)
135: {
136: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
138: PetscFunctionBegin;
139: PetscCallExternalVoid("STRUMPACK_compression_rel_tol", *rtol = (PetscReal)STRUMPACK_compression_rel_tol(*S));
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: static PetscErrorCode MatSTRUMPACKSetCompAbsTol_STRUMPACK(Mat F, PetscReal atol)
144: {
145: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
147: PetscFunctionBegin;
148: PetscCallExternalVoid("STRUMPACK_set_compression_abs_tol", STRUMPACK_set_compression_abs_tol(*S, atol));
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
151: static PetscErrorCode MatSTRUMPACKGetCompAbsTol_STRUMPACK(Mat F, PetscReal *atol)
152: {
153: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
155: PetscFunctionBegin;
156: PetscCallExternalVoid("STRUMPACK_compression_abs_tol", *atol = (PetscReal)STRUMPACK_compression_abs_tol(*S));
157: PetscFunctionReturn(PETSC_SUCCESS);
158: }
160: static PetscErrorCode MatSTRUMPACKSetCompLeafSize_STRUMPACK(Mat F, PetscInt leaf_size)
161: {
162: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
164: PetscFunctionBegin;
165: PetscCallExternalVoid("STRUMPACK_set_compression_leaf_size", STRUMPACK_set_compression_leaf_size(*S, leaf_size));
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
168: static PetscErrorCode MatSTRUMPACKGetCompLeafSize_STRUMPACK(Mat F, PetscInt *leaf_size)
169: {
170: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
172: PetscFunctionBegin;
173: PetscCallExternalVoid("STRUMPACK_compression_leaf_size", *leaf_size = (PetscInt)STRUMPACK_compression_leaf_size(*S));
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: static PetscErrorCode MatSTRUMPACKSetGeometricNxyz_STRUMPACK(Mat F, PetscInt nx, PetscInt ny, PetscInt nz)
178: {
179: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
181: PetscFunctionBegin;
182: if (nx < 1) {
183: PetscCheck(nx == PETSC_DECIDE || nx == PETSC_DEFAULT, PetscObjectComm((PetscObject)F), PETSC_ERR_ARG_OUTOFRANGE, "nx < 1");
184: nx = 1;
185: }
186: PetscCallExternalVoid("STRUMPACK_set_nx", STRUMPACK_set_nx(*S, nx));
187: if (ny < 1) {
188: PetscCheck(ny == PETSC_DECIDE || ny == PETSC_DEFAULT, PetscObjectComm((PetscObject)F), PETSC_ERR_ARG_OUTOFRANGE, "ny < 1");
189: ny = 1;
190: }
191: PetscCallExternalVoid("STRUMPACK_set_ny", STRUMPACK_set_ny(*S, ny));
192: if (nz < 1) {
193: PetscCheck(nz == PETSC_DECIDE || nz == PETSC_DEFAULT, PetscObjectComm((PetscObject)F), PETSC_ERR_ARG_OUTOFRANGE, "nz < 1");
194: nz = 1;
195: }
196: PetscCallExternalVoid("STRUMPACK_set_nz", STRUMPACK_set_nz(*S, nz));
197: PetscFunctionReturn(PETSC_SUCCESS);
198: }
199: static PetscErrorCode MatSTRUMPACKSetGeometricComponents_STRUMPACK(Mat F, PetscInt nc)
200: {
201: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
203: PetscFunctionBegin;
204: PetscCallExternalVoid("STRUMPACK_set_components", STRUMPACK_set_components(*S, nc));
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
207: static PetscErrorCode MatSTRUMPACKSetGeometricWidth_STRUMPACK(Mat F, PetscInt w)
208: {
209: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
211: PetscFunctionBegin;
212: PetscCallExternalVoid("STRUMPACK_set_separator_width", STRUMPACK_set_separator_width(*S, w));
213: PetscFunctionReturn(PETSC_SUCCESS);
214: }
216: static PetscErrorCode MatSTRUMPACKSetCompMinSepSize_STRUMPACK(Mat F, PetscInt min_sep_size)
217: {
218: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
220: PetscFunctionBegin;
221: PetscCallExternalVoid("STRUMPACK_set_compression_min_sep_size", STRUMPACK_set_compression_min_sep_size(*S, min_sep_size));
222: PetscFunctionReturn(PETSC_SUCCESS);
223: }
224: static PetscErrorCode MatSTRUMPACKGetCompMinSepSize_STRUMPACK(Mat F, PetscInt *min_sep_size)
225: {
226: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
228: PetscFunctionBegin;
229: PetscCallExternalVoid("STRUMPACK_compression_min_sep_size", *min_sep_size = (PetscInt)STRUMPACK_compression_min_sep_size(*S));
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
233: static PetscErrorCode MatSTRUMPACKSetCompLossyPrecision_STRUMPACK(Mat F, PetscInt lossy_prec)
234: {
235: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
237: PetscFunctionBegin;
238: PetscCallExternalVoid("STRUMPACK_set_compression_lossy_precision", STRUMPACK_set_compression_lossy_precision(*S, lossy_prec));
239: PetscFunctionReturn(PETSC_SUCCESS);
240: }
241: static PetscErrorCode MatSTRUMPACKGetCompLossyPrecision_STRUMPACK(Mat F, PetscInt *lossy_prec)
242: {
243: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
245: PetscFunctionBegin;
246: PetscCallExternalVoid("STRUMPACK_compression_lossy_precision", *lossy_prec = (PetscInt)STRUMPACK_compression_lossy_precision(*S));
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
250: static PetscErrorCode MatSTRUMPACKSetCompButterflyLevels_STRUMPACK(Mat F, PetscInt bfly_lvls)
251: {
252: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
254: PetscFunctionBegin;
255: PetscCallExternalVoid("STRUMPACK_set_compression_butterfly_levels", STRUMPACK_set_compression_butterfly_levels(*S, bfly_lvls));
256: PetscFunctionReturn(PETSC_SUCCESS);
257: }
258: static PetscErrorCode MatSTRUMPACKGetCompButterflyLevels_STRUMPACK(Mat F, PetscInt *bfly_lvls)
259: {
260: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
262: PetscFunctionBegin;
263: PetscCallExternalVoid("STRUMPACK_compression_butterfly_levels", *bfly_lvls = (PetscInt)STRUMPACK_compression_butterfly_levels(*S));
264: PetscFunctionReturn(PETSC_SUCCESS);
265: }
267: static PetscErrorCode MatSolve_STRUMPACK(Mat A, Vec b_mpi, Vec x)
268: {
269: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)A->data;
270: STRUMPACK_RETURN_CODE sp_err;
271: const PetscScalar *bptr;
272: PetscScalar *xptr;
274: PetscFunctionBegin;
275: PetscCall(VecGetArray(x, &xptr));
276: PetscCall(VecGetArrayRead(b_mpi, &bptr));
278: PetscCallExternalVoid("STRUMPACK_solve", sp_err = STRUMPACK_solve(*S, (PetscScalar *)bptr, xptr, 0));
279: switch (sp_err) {
280: case STRUMPACK_SUCCESS:
281: break;
282: case STRUMPACK_MATRIX_NOT_SET: {
283: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: matrix was not set");
284: break;
285: }
286: case STRUMPACK_REORDERING_ERROR: {
287: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: matrix reordering failed");
288: break;
289: }
290: default:
291: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: solve failed");
292: }
293: PetscCall(VecRestoreArray(x, &xptr));
294: PetscCall(VecRestoreArrayRead(b_mpi, &bptr));
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }
298: static PetscErrorCode MatMatSolve_STRUMPACK(Mat A, Mat B_mpi, Mat X)
299: {
300: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)A->data;
301: STRUMPACK_RETURN_CODE sp_err;
302: PetscBool flg;
303: PetscInt m = A->rmap->n, nrhs;
304: const PetscScalar *bptr;
305: PetscScalar *xptr;
307: PetscFunctionBegin;
308: PetscCall(PetscObjectTypeCompareAny((PetscObject)B_mpi, &flg, MATSEQDENSE, MATMPIDENSE, NULL));
309: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix B must be MATDENSE matrix");
310: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &flg, MATSEQDENSE, MATMPIDENSE, NULL));
311: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix X must be MATDENSE matrix");
313: PetscCall(MatGetSize(B_mpi, NULL, &nrhs));
314: PetscCall(MatDenseGetArray(X, &xptr));
315: PetscCall(MatDenseGetArrayRead(B_mpi, &bptr));
317: PetscCallExternalVoid("STRUMPACK_matsolve", sp_err = STRUMPACK_matsolve(*S, nrhs, bptr, m, xptr, m, 0));
318: switch (sp_err) {
319: case STRUMPACK_SUCCESS:
320: break;
321: case STRUMPACK_MATRIX_NOT_SET: {
322: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: matrix was not set");
323: break;
324: }
325: case STRUMPACK_REORDERING_ERROR: {
326: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: matrix reordering failed");
327: break;
328: }
329: default:
330: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: solve failed");
331: }
332: PetscCall(MatDenseRestoreArrayRead(B_mpi, &bptr));
333: PetscCall(MatDenseRestoreArray(X, &xptr));
334: PetscFunctionReturn(PETSC_SUCCESS);
335: }
337: static PetscErrorCode MatView_STRUMPACK(Mat A, PetscViewer viewer)
338: {
339: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)A->data;
340: MatSTRUMPACKCompressionType compressionType;
341: MatSTRUMPACKReordering reorderingType;
342: PetscBool cperm, gpu;
343: PetscBool isascii;
344: PetscViewerFormat format;
346: PetscFunctionBegin;
347: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
348: if (isascii) {
349: PetscCall(PetscViewerGetFormat(viewer, &format));
350: if (format == PETSC_VIEWER_ASCII_INFO) {
351: PetscCall(PetscViewerASCIIPrintf(viewer, "STRUMPACK sparse solver\n"));
352: PetscCallExternalVoid("STRUMPACK_compression", compressionType = (MatSTRUMPACKCompressionType)STRUMPACK_compression(*S));
353: PetscCall(PetscViewerASCIIPrintf(viewer, " Compression type %s\n", MatSTRUMPACKCompressionTypes[compressionType]));
354: PetscCallExternalVoid("STRUMPACK_reordering_method", reorderingType = (MatSTRUMPACKReordering)STRUMPACK_reordering_method(*S));
355: PetscCall(PetscViewerASCIIPrintf(viewer, " Reordering %s\n", MatSTRUMPACKReorderingTypes[reorderingType]));
356: PetscCallExternalVoid("STRUMPACK_matching", cperm = (PetscBool)(STRUMPACK_matching(*S) != STRUMPACK_MATCHING_NONE));
357: PetscCall(PetscViewerASCIIPrintf(viewer, " Column permutation: %s\n", cperm ? "True" : "False"));
358: PetscCallExternalVoid("STRUMPACK_use_gpu", gpu = (PetscBool)STRUMPACK_use_gpu(*S));
359: PetscCall(PetscViewerASCIIPrintf(viewer, " Use GPU: %s\n", gpu ? "True" : "False"));
360: }
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 Aloc;
370: const PetscScalar *av;
371: const PetscInt *ai = NULL, *aj = NULL;
372: PetscInt M = A->rmap->N, m = A->rmap->n, dummy;
373: PetscBool ismpiaij, isseqaij, flg;
375: PetscFunctionBegin;
376: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &isseqaij));
377: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
378: if (ismpiaij) PetscCall(MatMPIAIJGetLocalMat(A, MAT_INITIAL_MATRIX, &Aloc));
379: else if (isseqaij) {
380: PetscCall(PetscObjectReference((PetscObject)A));
381: Aloc = A;
382: } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for type %s", ((PetscObject)A)->type_name);
384: PetscCall(MatGetRowIJ(Aloc, 0, PETSC_FALSE, PETSC_FALSE, &dummy, &ai, &aj, &flg));
385: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "GetRowIJ failed");
386: PetscCall(MatSeqAIJGetArrayRead(Aloc, &av));
388: if (ismpiaij) {
389: const PetscInt *dist = NULL;
390: PetscCall(MatGetOwnershipRanges(A, &dist));
391: PetscCallExternalVoid("STRUMPACK_set_distributed_csr_matrix", STRUMPACK_set_distributed_csr_matrix(*S, &m, ai, aj, av, dist, 0));
392: } else if (isseqaij) {
393: PetscCallExternalVoid("STRUMPACK_set_csr_matrix", STRUMPACK_set_csr_matrix(*S, &M, ai, aj, av, 0));
394: }
396: PetscCall(MatRestoreRowIJ(Aloc, 0, PETSC_FALSE, PETSC_FALSE, &dummy, &ai, &aj, &flg));
397: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "RestoreRowIJ failed");
398: PetscCall(MatSeqAIJRestoreArrayRead(Aloc, &av));
399: PetscCall(MatDestroy(&Aloc));
401: /* Reorder and Factor the matrix. */
402: /* TODO figure out how to avoid reorder if the matrix values changed, but the pattern remains the same. */
403: PetscCallExternalVoid("STRUMPACK_reorder", sp_err = STRUMPACK_reorder(*S));
404: PetscCheck(sp_err != STRUMPACK_MATRIX_NOT_SET, PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: matrix was not set");
405: PetscCheck(sp_err != STRUMPACK_REORDERING_ERROR, PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: matrix reordering failed");
406: PetscCheck(sp_err == STRUMPACK_SUCCESS, PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: unknown error while reordering");
407: PetscCallExternalVoid("STRUMPACK_factor", sp_err = STRUMPACK_factor(*S));
408: PetscCheck(sp_err == STRUMPACK_SUCCESS, PETSC_COMM_SELF, PETSC_ERR_LIB, "STRUMPACK error: factorization failed");
409: F->assembled = PETSC_TRUE;
410: F->preallocated = PETSC_TRUE;
411: PetscFunctionReturn(PETSC_SUCCESS);
412: }
414: static PetscErrorCode MatLUFactorSymbolic_STRUMPACK(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
415: {
416: PetscBool flg, set, verb;
417: PetscReal ctol;
418: PetscInt min_sep_size, leaf_size, nxyz[3], nrdims, nc, w;
419: #if defined(STRUMPACK_USE_ZFP)
420: PetscInt lossy_prec;
421: #endif
422: #if defined(STRUMPACK_USE_BPACK)
423: PetscInt bfly_lvls;
424: #endif
425: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver *)F->data;
426: STRUMPACK_REORDERING_STRATEGY ndcurrent, ndvalue;
427: STRUMPACK_COMPRESSION_TYPE compcurrent, compvalue;
429: PetscFunctionBegin;
430: PetscOptionsBegin(PetscObjectComm((PetscObject)F), ((PetscObject)F)->prefix, "STRUMPACK Options", "Mat");
432: PetscCall(PetscOptionsBool("-mat_strumpack_verbose", "Print Strumpack's verbose output", "", PETSC_FALSE, &verb, &flg));
433: if (flg) PetscCallExternalVoid("STRUMPACK_set_verbose", STRUMPACK_set_verbose(*S, verb));
435: /* By default, no compression is done. Compression is enabled when the user enables it with */
436: /* -mat_strumpack_compression with anything else than NONE, or when selecting ilu */
437: /* preconditioning, in which case we default to STRUMPACK_BLR compression. */
438: /* When compression is enabled, the STRUMPACK solver becomes an incomplete */
439: /* (or approximate) LU factorization. */
440: PetscCallExternalVoid("STRUMPACK_compression", compcurrent = STRUMPACK_compression(*S));
441: PetscCall(PetscOptionsEnum("-mat_strumpack_compression", "Rank-structured compression type", "None", MatSTRUMPACKCompressionTypes, (PetscEnum)compcurrent, (PetscEnum *)&compvalue, &set));
442: if (set) PetscCall(MatSTRUMPACKSetCompression(F, (MatSTRUMPACKCompressionType)compvalue));
443: else if (F->factortype == MAT_FACTOR_ILU) PetscCallExternalVoid("STRUMPACK_set_compression", STRUMPACK_set_compression(*S, STRUMPACK_BLR));
445: PetscCallExternalVoid("STRUMPACK_compression_rel_tol", ctol = (PetscReal)STRUMPACK_compression_rel_tol(*S));
446: PetscCall(PetscOptionsReal("-mat_strumpack_compression_rel_tol", "Relative compression tolerance", "None", ctol, &ctol, &set));
447: if (set) PetscCallExternalVoid("STRUMPACK_set_compression_rel_tol", STRUMPACK_set_compression_rel_tol(*S, (double)ctol));
449: PetscCallExternalVoid("STRUMPACK_compression_abs_tol", ctol = (PetscReal)STRUMPACK_compression_abs_tol(*S));
450: PetscCall(PetscOptionsReal("-mat_strumpack_compression_abs_tol", "Absolute compression tolerance", "None", ctol, &ctol, &set));
451: if (set) PetscCallExternalVoid("STRUMPACK_set_compression_abs_tol", STRUMPACK_set_compression_abs_tol(*S, (double)ctol));
453: PetscCallExternalVoid("STRUMPACK_compression_min_sep_size", min_sep_size = (PetscInt)STRUMPACK_compression_min_sep_size(*S));
454: PetscCall(PetscOptionsInt("-mat_strumpack_compression_min_sep_size", "Minimum size of separator for compression", "None", min_sep_size, &min_sep_size, &set));
455: if (set) PetscCallExternalVoid("STRUMPACK_set_compression_min_sep_size", STRUMPACK_set_compression_min_sep_size(*S, (int)min_sep_size));
457: PetscCallExternalVoid("STRUMPACK_compression_leaf_size", leaf_size = (PetscInt)STRUMPACK_compression_leaf_size(*S));
458: PetscCall(PetscOptionsInt("-mat_strumpack_compression_leaf_size", "Size of diagonal blocks in rank-structured approximation", "None", leaf_size, &leaf_size, &set));
459: if (set) PetscCallExternalVoid("STRUMPACK_set_compression_leaf_size", STRUMPACK_set_compression_leaf_size(*S, (int)leaf_size));
461: #if defined(STRUMPACK_USE_ZFP)
462: PetscCallExternalVoid("STRUMPACK_compression_lossy_precision", lossy_prec = (PetscInt)STRUMPACK_compression_lossy_precision(*S));
463: PetscCall(PetscOptionsInt("-mat_strumpack_compression_lossy_precision", "Number of bitplanes to use in lossy compression", "None", lossy_prec, &lossy_prec, &set));
464: if (set) PetscCallExternalVoid("STRUMPACK_set_compression_lossy_precision", STRUMPACK_set_compression_lossy_precision(*S, (int)lossy_prec));
465: #endif
467: #if defined(STRUMPACK_USE_BPACK)
468: PetscCallExternalVoid("STRUMPACK_compression_butterfly_levels", bfly_lvls = (PetscInt)STRUMPACK_compression_butterfly_levels(*S));
469: 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));
470: if (set) PetscCallExternalVoid("STRUMPACK_set_compression_butterfly_levels", STRUMPACK_set_compression_butterfly_levels(*S, (int)bfly_lvls));
471: #endif
473: #if defined(STRUMPACK_USE_CUDA) || defined(STRUMPACK_USE_HIP)
474: PetscCallExternalVoid("STRUMPACK_use_gpu", flg = (STRUMPACK_use_gpu(*S) == 0) ? PETSC_FALSE : PETSC_TRUE);
475: PetscCall(PetscOptionsBool("-mat_strumpack_gpu", "Enable GPU acceleration (not supported for all compression types)", "None", flg, &flg, &set));
476: if (set) PetscCall(MatSTRUMPACKSetGPU(F, flg));
477: #endif
479: PetscCallExternalVoid("STRUMPACK_matching", flg = (STRUMPACK_matching(*S) == STRUMPACK_MATCHING_NONE) ? PETSC_FALSE : PETSC_TRUE);
480: PetscCall(PetscOptionsBool("-mat_strumpack_colperm", "Find a col perm to get nonzero diagonal", "None", flg, &flg, &set));
481: if (set) PetscCallExternalVoid("STRUMPACK_set_matching", STRUMPACK_set_matching(*S, flg ? STRUMPACK_MATCHING_MAX_DIAGONAL_PRODUCT_SCALING : STRUMPACK_MATCHING_NONE));
483: PetscCallExternalVoid("STRUMPACK_reordering_method", ndcurrent = STRUMPACK_reordering_method(*S));
484: PetscCall(PetscOptionsEnum("-mat_strumpack_reordering", "Sparsity reducing matrix reordering", "None", MatSTRUMPACKReorderingTypes, (PetscEnum)ndcurrent, (PetscEnum *)&ndvalue, &set));
485: if (set) PetscCallExternalVoid("STRUMPACK_set_reordering_method", STRUMPACK_set_reordering_method(*S, ndvalue));
487: /* geometric ordering, for a regular 1D/2D/3D mesh in the natural ordering, */
488: /* with nc DOF's per gridpoint, and possibly a wider stencil */
489: nrdims = 3;
490: nxyz[0] = nxyz[1] = nxyz[2] = 1;
491: PetscCall(PetscOptionsIntArray("-mat_strumpack_geometric_xyz", "Mesh sizes nx,ny,nz (Use 1 for default)", "", nxyz, &nrdims, &set));
492: if (set) {
493: PetscCheck(nrdims == 1 || nrdims == 2 || nrdims == 3, PetscObjectComm((PetscObject)F), PETSC_ERR_ARG_OUTOFRANGE, "-mat_strumpack_geometric_xyz requires 1, 2, or 3 values");
494: PetscCall(MatSTRUMPACKSetGeometricNxyz(F, nxyz[0], nxyz[1], nxyz[2]));
495: }
496: PetscCall(PetscOptionsInt("-mat_strumpack_geometric_components", "Number of components per mesh point, for geometric nested dissection ordering", "None", 1, &nc, &set));
497: if (set) PetscCallExternalVoid("STRUMPACK_set_components", STRUMPACK_set_components(*S, (int)nc));
498: 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));
499: if (set) PetscCallExternalVoid("STRUMPACK_set_separator_width", STRUMPACK_set_separator_width(*S, (int)w));
501: PetscCallExternalVoid("STRUMPACK_use_METIS_NodeNDP", flg = (STRUMPACK_use_METIS_NodeNDP(*S) == 0) ? PETSC_FALSE : PETSC_TRUE);
502: PetscCall(PetscOptionsBool("-mat_strumpack_metis_nodeNDP", "Use METIS_NodeNDP instead of METIS_NodeND, for a more balanced tree", "None", flg, &flg, &set));
503: if (set) {
504: if (flg) PetscCallExternalVoid("STRUMPACK_enable_METIS_NodeNDP", STRUMPACK_enable_METIS_NodeNDP(*S));
505: else PetscCallExternalVoid("STRUMPACK_disable_METIS_NodeNDP", STRUMPACK_disable_METIS_NodeNDP(*S));
506: }
508: /* Disable the outer iterative solver from STRUMPACK. */
509: /* When STRUMPACK is used as a direct solver, it will by default do iterative refinement. */
510: /* When STRUMPACK is used as an approximate factorization preconditioner (by enabling */
511: /* low-rank compression), it will use its own preconditioned GMRES. Here we can disable */
512: /* the outer iterative solver, as PETSc uses STRUMPACK from within a KSP. */
513: PetscCallExternalVoid("STRUMPACK_set_Krylov_solver", STRUMPACK_set_Krylov_solver(*S, STRUMPACK_DIRECT));
514: PetscOptionsEnd();
515: PetscFunctionReturn(PETSC_SUCCESS);
516: }
518: static PetscErrorCode MatFactorGetSolverType_aij_strumpack(Mat A, MatSolverType *type)
519: {
520: PetscFunctionBegin;
521: *type = MATSOLVERSTRUMPACK;
522: PetscFunctionReturn(PETSC_SUCCESS);
523: }
525: /*MC
526: MATSOLVERSTRUMPACK = "strumpack" - A solver package providing a direct sparse solver (`PCLU`)
527: and a preconditioner (`PCILU`) using low-rank compression via the external package STRUMPACK <https://portal.nersc.gov/project/sparse/strumpack/master>.
529: Use `./configure --download-strumpack --download-metis` to have PETSc installed with STRUMPACK.
531: For full functionality, add `--download-slate --download-magma --download-parmetis --download-ptscotch --download-zfp --download-butterflypack`.
533: SLATE provides GPU support in the multi-GPU setting, providing ScaLAPACK functionality but with GPU acceleration.
535: MAGMA can optionally be used for on node GPU support instead cuBLAS/cuSOLVER, and performs slightly better.
537: ParMETIS and PTScotch can be used for parallel fill-reducing ordering.
539: ZFP is used for floating point compression of the sparse factors (`lossy` or `lossless` compression).
541: ButterflyPACK is used for `MAT_STRUMPACK_COMPRESSION_TYPE_HODLR` (Hierarchically Off-Diagonal Low Rank) and Hierarchically Off-Diagonal Butterfly
542: compression of the sparse factors.
544: Options Database Keys:
545: + -mat_strumpack_verbose (true|false) - Enable verbose output
546: . -mat_strumpack_compression (none|hss|blr|hodlr|blr_hodlr|zfp_blr_hodlr|lossless|lossy) - Type of rank-structured compression in sparse LU factors
547: . -mat_strumpack_compression_rel_tol rel_tol - Relative compression tolerance, when using `-pctype ilu`
548: . -mat_strumpack_compression_abs_tol abs_tol - Absolute compression tolerance, when using `-pctype ilu`
549: . -mat_strumpack_compression_min_sep_size min_sep - Minimum size of separator for rank-structured compression, when using `-pctype ilu`
550: . -mat_strumpack_compression_leaf_size size - Size of diagonal blocks in rank-structured approximation, when using `-pctype ilu`
551: . -mat_strumpack_compression_lossy_precision pre - Precision when using lossy compression [1-64], when using `-pctype ilu`, `MAT_STRUMPACK_COMPRESSION_TYPE_LOSSY` compression (requires ZFP support)
552: . -mat_strumpack_compression_butterfly_levels levels - Number of levels in the hierarchically off-diagonal matrix for which to use butterfly, when using `-pctype ilu`, `MAT_STRUMPACK_COMPRESSION_TYPE_BLR_HODLR` compression (requires ButterflyPACK support)
553: . -mat_strumpack_gpu (true|false) - Enable GPU acceleration in numerical factorization (not supported for all compression types)
554: . -mat_strumpack_colperm (true|false) - Permute matrix to make diagonal nonzeros
555: . -mat_strumpack_reordering (natural|metis|parmetis|scotch|ptscotch|rcm|geometric|amd|mmd|and|mlf|spectral) - Sparsity reducing matrix reordering
556: . -mat_strumpack_geometric_xyz m,n,p - Mesh x,y,z dimensions, for use with `MAT_STRUMPACK_GEOMETRIC` ordering
557: . -mat_strumpack_geometric_components dof - Number of components per mesh point, for `MAT_STRUMPACK_GEOMETRIC` ordering
558: . -mat_strumpack_geometric_width width - Width of the separator of the mesh, for geometric nested dissection ordering
559: - -mat_strumpack_metis_nodeNDP - Use METIS_NodeNDP instead of METIS_NodeND, for a more balanced tree
561: Level: beginner
563: Notes:
564: Recommended use is 1 MPI process per GPU.
566: Use `-pc_type lu` `-pc_factor_mat_solver_type strumpack` to use this as an exact (direct) solver.
568: 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).
570: Works with `MATAIJ` matrices
572: `MAT_STRUMPACK_COMPRESSION_TYPE_HODLR`, `MAT_STRUMPACK_COMPRESSION_TYPE_BLR_HODLR`, and `MAT_STRUMPACK_COMPRESSION_TYPE_ZFP_BLR_HODLR` compression
573: require STRUMPACK to be configured with ButterflyPACK support (`--download-butterflypack`).
575: `MAT_STRUMPACK_COMPRESSION_TYPE_LOSSLESS`, `MAT_STRUMPACK_COMPRESSION_TYPE_LOSSY`, and `MAT_STRUMPACK_COMPRESSION_TYPE_ZFP_BLR_HODLR` compression require
576: STRUMPACK to be configured with ZFP support (`--download-zfp`).
578: Developer Note:
579: `MatView_Info_STRUMPACK()` should display all the STRUMPACK options used
581: .seealso: `MATSOLVERSTRUMPACK`, [](ch_matrices), `Mat`, `PCLU`, `PCILU`, `MATSOLVERSUPERLU_DIST`, `MATSOLVERMUMPS`, `PCFactorSetMatSolverType()`, `MatSolverType`,
582: `MatGetFactor()`, `MatSTRUMPACKSetReordering()`, `MatSTRUMPACKReordering`, `MatSTRUMPACKCompressionType`, `MatSTRUMPACKSetColPerm()`
583: M*/
584: static PetscErrorCode MatGetFactor_aij_strumpack(Mat A, MatFactorType ftype, Mat *F)
585: {
586: Mat B;
587: PetscBool flg;
588: #if defined(STRUMPACK_USE_SLATE_SCALAPACK)
589: PetscMPIInt mpithreads;
590: #endif
591: STRUMPACK_SparseSolver *S;
592: STRUMPACK_INTERFACE iface;
593: const STRUMPACK_PRECISION table[2][2][2] = {
594: {{STRUMPACK_FLOATCOMPLEX_64, STRUMPACK_DOUBLECOMPLEX_64}, {STRUMPACK_FLOAT_64, STRUMPACK_DOUBLE_64}},
595: {{STRUMPACK_FLOATCOMPLEX, STRUMPACK_DOUBLECOMPLEX}, {STRUMPACK_FLOAT, STRUMPACK_DOUBLE} }
596: };
597: const STRUMPACK_PRECISION prec = table[(sizeof(PetscInt) == 8) ? 0 : 1][(PETSC_SCALAR == PETSC_COMPLEX) ? 0 : 1][(PETSC_REAL == PETSC_FLOAT) ? 0 : 1];
599: PetscFunctionBegin;
600: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
601: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
602: PetscCall(PetscStrallocpy("strumpack", &((PetscObject)B)->type_name));
603: PetscCall(MatSetUp(B));
604: PetscCall(MatSeqAIJSetPreallocation(B, 0, NULL));
605: PetscCall(MatMPIAIJSetPreallocation(B, 0, NULL, 0, NULL));
606: PetscCheck(ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU, PETSC_COMM_SELF, PETSC_ERR_SUP, "Factor type not supported");
607: B->ops->lufactorsymbolic = MatLUFactorSymbolic_STRUMPACK;
608: B->ops->ilufactorsymbolic = MatLUFactorSymbolic_STRUMPACK;
609: B->ops->getinfo = MatGetInfo_External;
610: B->ops->view = MatView_STRUMPACK;
611: B->ops->destroy = MatDestroy_STRUMPACK;
612: B->ops->getdiagonal = MatGetDiagonal_STRUMPACK;
613: B->ops->lufactornumeric = MatLUFactorNumeric_STRUMPACK;
614: B->ops->solve = MatSolve_STRUMPACK;
615: B->ops->matsolve = MatMatSolve_STRUMPACK;
616: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_aij_strumpack));
617: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetReordering_C", MatSTRUMPACKSetReordering_STRUMPACK));
618: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKGetReordering_C", MatSTRUMPACKGetReordering_STRUMPACK));
619: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetColPerm_C", MatSTRUMPACKSetColPerm_STRUMPACK));
620: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKGetColPerm_C", MatSTRUMPACKGetColPerm_STRUMPACK));
621: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetGeometricNxyz_C", MatSTRUMPACKSetGeometricNxyz_STRUMPACK));
622: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetGeometricComponents_C", MatSTRUMPACKSetGeometricComponents_STRUMPACK));
623: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetGeometricWidth_C", MatSTRUMPACKSetGeometricWidth_STRUMPACK));
624: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetGPU_C", MatSTRUMPACKSetGPU_STRUMPACK));
625: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKGetGPU_C", MatSTRUMPACKGetGPU_STRUMPACK));
626: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetCompression_C", MatSTRUMPACKSetCompression_STRUMPACK));
627: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKGetCompression_C", MatSTRUMPACKGetCompression_STRUMPACK));
628: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetCompRelTol_C", MatSTRUMPACKSetCompRelTol_STRUMPACK));
629: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKGetCompRelTol_C", MatSTRUMPACKGetCompRelTol_STRUMPACK));
630: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetCompAbsTol_C", MatSTRUMPACKSetCompAbsTol_STRUMPACK));
631: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKGetCompAbsTol_C", MatSTRUMPACKGetCompAbsTol_STRUMPACK));
632: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetCompLeafSize_C", MatSTRUMPACKSetCompLeafSize_STRUMPACK));
633: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKGetCompLeafSize_C", MatSTRUMPACKGetCompLeafSize_STRUMPACK));
634: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetCompMinSepSize_C", MatSTRUMPACKSetCompMinSepSize_STRUMPACK));
635: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKGetCompMinSepSize_C", MatSTRUMPACKGetCompMinSepSize_STRUMPACK));
636: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetCompLossyPrecision_C", MatSTRUMPACKSetCompLossyPrecision_STRUMPACK));
637: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKGetCompLossyPrecision_C", MatSTRUMPACKGetCompLossyPrecision_STRUMPACK));
638: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKSetCompButterflyLevels_C", MatSTRUMPACKSetCompButterflyLevels_STRUMPACK));
639: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSTRUMPACKGetCompButterflyLevels_C", MatSTRUMPACKGetCompButterflyLevels_STRUMPACK));
640: B->factortype = ftype;
642: PetscCall(PetscFree(B->solvertype));
643: PetscCall(PetscStrallocpy(MATSOLVERSTRUMPACK, &B->solvertype));
645: #if defined(STRUMPACK_USE_SLATE_SCALAPACK)
646: PetscCallMPI(MPI_Query_thread(&mpithreads));
647: PetscCheck(mpithreads == MPI_THREAD_MULTIPLE, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP_SYS, "SLATE requires MPI_THREAD_MULTIPLE");
648: #endif
649: PetscCall(PetscNew(&S));
650: B->data = S;
651: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &flg)); /* A might be MATSEQAIJCUSPARSE */
652: iface = flg ? STRUMPACK_MT : STRUMPACK_MPI_DIST;
653: PetscCallExternalVoid("STRUMPACK_init", STRUMPACK_init(S, PetscObjectComm((PetscObject)A), prec, iface, 0, NULL, PETSC_FALSE));
654: *F = B;
655: PetscFunctionReturn(PETSC_SUCCESS);
656: }
658: PETSC_INTERN PetscErrorCode MatSolverTypeRegister_STRUMPACK(void)
659: {
660: PetscFunctionBegin;
661: PetscCall(MatSolverTypeRegister(MATSOLVERSTRUMPACK, MATMPIAIJ, MAT_FACTOR_LU, MatGetFactor_aij_strumpack));
662: PetscCall(MatSolverTypeRegister(MATSOLVERSTRUMPACK, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_aij_strumpack));
663: PetscCall(MatSolverTypeRegister(MATSOLVERSTRUMPACK, MATMPIAIJ, MAT_FACTOR_ILU, MatGetFactor_aij_strumpack));
664: PetscCall(MatSolverTypeRegister(MATSOLVERSTRUMPACK, MATSEQAIJ, MAT_FACTOR_ILU, MatGetFactor_aij_strumpack));
665: PetscFunctionReturn(PETSC_SUCCESS);
666: }