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: }