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: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Mat type: STRUMPACK factor");
8: return 0;
9: }
11: static PetscErrorCode MatDestroy_STRUMPACK(Mat A)
12: {
13: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)A->spptr;
14: PetscBool flg;
16: /* Deallocate STRUMPACK storage */
17: PetscStackCall("STRUMPACK_destroy",STRUMPACK_destroy(S));
18: PetscFree(A->spptr);
19: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
20: if (flg) {
21: MatDestroy_SeqAIJ(A);
22: } else {
23: MatDestroy_MPIAIJ(A);
24: }
26: /* clear composed functions */
27: PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL);
28: PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetReordering_C",NULL);
29: PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetColPerm_C",NULL);
30: PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSRelTol_C",NULL);
31: PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSAbsTol_C",NULL);
32: PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSMaxRank_C",NULL);
33: PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSLeafSize_C",NULL);
34: PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSMinSepSize_C",NULL);
36: return 0;
37: }
39: static PetscErrorCode MatSTRUMPACKSetReordering_STRUMPACK(Mat F,MatSTRUMPACKReordering reordering)
40: {
41: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
43: PetscStackCall("STRUMPACK_reordering_method",STRUMPACK_set_reordering_method(*S,(STRUMPACK_REORDERING_STRATEGY)reordering));
44: return 0;
45: }
47: /*@
48: MatSTRUMPACKSetReordering - Set STRUMPACK fill-reducing reordering
50: Input Parameters:
51: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
52: - reordering - the code to be used to find the fill-reducing reordering
53: Possible values: NATURAL=0 METIS=1 PARMETIS=2 SCOTCH=3 PTSCOTCH=4 RCM=5
55: Options Database:
56: . -mat_strumpack_reordering <METIS> - Sparsity reducing matrix reordering (choose one of) NATURAL METIS PARMETIS SCOTCH PTSCOTCH RCM (None)
58: Level: beginner
60: References:
61: . * - STRUMPACK manual
63: .seealso: MatGetFactor()
64: @*/
65: PetscErrorCode MatSTRUMPACKSetReordering(Mat F,MatSTRUMPACKReordering reordering)
66: {
69: PetscTryMethod(F,"MatSTRUMPACKSetReordering_C",(Mat,MatSTRUMPACKReordering),(F,reordering));
70: return 0;
71: }
73: static PetscErrorCode MatSTRUMPACKSetColPerm_STRUMPACK(Mat F,PetscBool cperm)
74: {
75: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
77: PetscStackCall("STRUMPACK_set_mc64job",STRUMPACK_set_mc64job(*S,cperm ? 5 : 0));
78: return 0;
79: }
81: /*@
82: MatSTRUMPACKSetColPerm - Set whether STRUMPACK should try to permute the columns of the matrix in order to get a nonzero diagonal
84: Logically Collective on Mat
86: Input Parameters:
87: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
88: - cperm - PETSC_TRUE to permute (internally) the columns of the matrix
90: Options Database:
91: . -mat_strumpack_colperm <cperm> - true to use the permutation
93: Level: beginner
95: References:
96: . * - STRUMPACK manual
98: .seealso: MatGetFactor()
99: @*/
100: PetscErrorCode MatSTRUMPACKSetColPerm(Mat F,PetscBool cperm)
101: {
104: PetscTryMethod(F,"MatSTRUMPACKSetColPerm_C",(Mat,PetscBool),(F,cperm));
105: return 0;
106: }
108: static PetscErrorCode MatSTRUMPACKSetHSSRelTol_STRUMPACK(Mat F,PetscReal rtol)
109: {
110: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
112: PetscStackCall("STRUMPACK_set_HSS_rel_tol", STRUMPACK_set_HSS_rel_tol(*S,rtol));
113: return 0;
114: }
116: /*@
117: MatSTRUMPACKSetHSSRelTol - Set STRUMPACK relative tolerance for HSS compression
119: Logically Collective on Mat
121: Input Parameters:
122: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
123: - rtol - relative compression tolerance
125: Options Database:
126: . -mat_strumpack_hss_rel_tol <1e-2> - Relative compression tolerance (None)
128: Level: beginner
130: References:
131: . * - STRUMPACK manual
133: .seealso: MatGetFactor()
134: @*/
135: PetscErrorCode MatSTRUMPACKSetHSSRelTol(Mat F,PetscReal rtol)
136: {
139: PetscTryMethod(F,"MatSTRUMPACKSetHSSRelTol_C",(Mat,PetscReal),(F,rtol));
140: return 0;
141: }
143: static PetscErrorCode MatSTRUMPACKSetHSSAbsTol_STRUMPACK(Mat F,PetscReal atol)
144: {
145: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
147: PetscStackCall("STRUMPACK_set_HSS_abs_tol", STRUMPACK_set_HSS_abs_tol(*S,atol));
148: return 0;
149: }
151: /*@
152: MatSTRUMPACKSetHSSAbsTol - Set STRUMPACK absolute tolerance for HSS compression
154: Logically Collective on Mat
156: Input Parameters:
157: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
158: - atol - absolute compression tolerance
160: Options Database:
161: . -mat_strumpack_hss_abs_tol <1e-8> - Absolute compression tolerance (None)
163: Level: beginner
165: References:
166: . * - STRUMPACK manual
168: .seealso: MatGetFactor()
169: @*/
170: PetscErrorCode MatSTRUMPACKSetHSSAbsTol(Mat F,PetscReal atol)
171: {
174: PetscTryMethod(F,"MatSTRUMPACKSetHSSAbsTol_C",(Mat,PetscReal),(F,atol));
175: return 0;
176: }
178: static PetscErrorCode MatSTRUMPACKSetHSSMaxRank_STRUMPACK(Mat F,PetscInt hssmaxrank)
179: {
180: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
182: PetscStackCall("STRUMPACK_set_HSS_max_rank", STRUMPACK_set_HSS_max_rank(*S,hssmaxrank));
183: return 0;
184: }
186: /*@
187: MatSTRUMPACKSetHSSMaxRank - Set STRUMPACK maximum HSS rank
189: Logically Collective on Mat
191: Input Parameters:
192: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
193: - hssmaxrank - maximum rank used in low-rank approximation
195: Options Database:
196: . -mat_strumpack_max_rank - Maximum rank in HSS compression, when using pctype ilu (None)
198: Level: beginner
200: References:
201: . * - STRUMPACK manual
203: .seealso: MatGetFactor()
204: @*/
205: PetscErrorCode MatSTRUMPACKSetHSSMaxRank(Mat F,PetscInt hssmaxrank)
206: {
209: PetscTryMethod(F,"MatSTRUMPACKSetHSSMaxRank_C",(Mat,PetscInt),(F,hssmaxrank));
210: return 0;
211: }
213: static PetscErrorCode MatSTRUMPACKSetHSSLeafSize_STRUMPACK(Mat F,PetscInt leaf_size)
214: {
215: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
217: PetscStackCall("STRUMPACK_set_HSS_leaf_size", STRUMPACK_set_HSS_leaf_size(*S,leaf_size));
218: return 0;
219: }
221: /*@
222: MatSTRUMPACKSetHSSLeafSize - Set STRUMPACK HSS leaf size
224: Logically Collective on Mat
226: Input Parameters:
227: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
228: - leaf_size - Size of diagonal blocks in HSS approximation
230: Options Database:
231: . -mat_strumpack_leaf_size - Size of diagonal blocks in HSS approximation, when using pctype ilu (None)
233: Level: beginner
235: References:
236: . * - STRUMPACK manual
238: .seealso: MatGetFactor()
239: @*/
240: PetscErrorCode MatSTRUMPACKSetHSSLeafSize(Mat F,PetscInt leaf_size)
241: {
244: PetscTryMethod(F,"MatSTRUMPACKSetHSSLeafSize_C",(Mat,PetscInt),(F,leaf_size));
245: return 0;
246: }
248: static PetscErrorCode MatSTRUMPACKSetHSSMinSepSize_STRUMPACK(Mat F,PetscInt hssminsize)
249: {
250: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
252: PetscStackCall("STRUMPACK_set_HSS_min_sep_size", STRUMPACK_set_HSS_min_sep_size(*S,hssminsize));
253: return 0;
254: }
256: /*@
257: MatSTRUMPACKSetHSSMinSepSize - Set STRUMPACK minimum separator size for low-rank approximation
259: Logically Collective on Mat
261: Input Parameters:
262: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
263: - hssminsize - minimum dense matrix size for low-rank approximation
265: Options Database:
266: . -mat_strumpack_hss_min_sep_size <hssminsize> - set the minimum separator size
268: Level: beginner
270: References:
271: . * - STRUMPACK manual
273: .seealso: MatGetFactor()
274: @*/
275: PetscErrorCode MatSTRUMPACKSetHSSMinSepSize(Mat F,PetscInt hssminsize)
276: {
279: PetscTryMethod(F,"MatSTRUMPACKSetHSSMinSepSize_C",(Mat,PetscInt),(F,hssminsize));
280: return 0;
281: }
283: static PetscErrorCode MatSolve_STRUMPACK(Mat A,Vec b_mpi,Vec x)
284: {
285: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)A->spptr;
286: STRUMPACK_RETURN_CODE sp_err;
287: const PetscScalar *bptr;
288: PetscScalar *xptr;
290: VecGetArray(x,&xptr);
291: VecGetArrayRead(b_mpi,&bptr);
293: PetscStackCall("STRUMPACK_solve",sp_err = STRUMPACK_solve(*S,(PetscScalar*)bptr,xptr,0));
294: switch (sp_err) {
295: case STRUMPACK_SUCCESS: break;
296: case STRUMPACK_MATRIX_NOT_SET: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix was not set"); break; }
297: case STRUMPACK_REORDERING_ERROR: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix reordering failed"); break; }
298: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: solve failed");
299: }
300: VecRestoreArray(x,&xptr);
301: VecRestoreArrayRead(b_mpi,&bptr);
302: return 0;
303: }
305: static PetscErrorCode MatMatSolve_STRUMPACK(Mat A,Mat B_mpi,Mat X)
306: {
307: PetscBool flg;
309: PetscObjectTypeCompareAny((PetscObject)B_mpi,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
311: PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
313: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatSolve_STRUMPACK() is not implemented yet");
314: return 0;
315: }
317: static PetscErrorCode MatView_Info_STRUMPACK(Mat A,PetscViewer viewer)
318: {
319: /* check if matrix is strumpack type */
320: if (A->ops->solve != MatSolve_STRUMPACK) return 0;
321: PetscViewerASCIIPrintf(viewer,"STRUMPACK sparse solver!\n");
322: return 0;
323: }
325: static PetscErrorCode MatView_STRUMPACK(Mat A,PetscViewer viewer)
326: {
327: PetscBool iascii;
328: PetscViewerFormat format;
330: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
331: if (iascii) {
332: PetscViewerGetFormat(viewer,&format);
333: if (format == PETSC_VIEWER_ASCII_INFO) MatView_Info_STRUMPACK(A,viewer);
334: }
335: return 0;
336: }
338: static PetscErrorCode MatLUFactorNumeric_STRUMPACK(Mat F,Mat A,const MatFactorInfo *info)
339: {
340: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
341: STRUMPACK_RETURN_CODE sp_err;
342: Mat_SeqAIJ *A_d,*A_o;
343: Mat_MPIAIJ *mat;
344: PetscInt M=A->rmap->N,m=A->rmap->n;
345: PetscBool flg;
347: PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);
348: if (flg) { /* A is MATMPIAIJ */
349: mat = (Mat_MPIAIJ*)A->data;
350: A_d = (Mat_SeqAIJ*)(mat->A)->data;
351: A_o = (Mat_SeqAIJ*)(mat->B)->data;
352: PetscStackCall("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));
353: } else { /* A is MATSEQAIJ */
354: A_d = (Mat_SeqAIJ*)A->data;
355: PetscStackCall("STRUMPACK_set_csr_matrix",STRUMPACK_set_csr_matrix(*S,&M,A_d->i,A_d->j,A_d->a,0));
356: }
358: /* Reorder and Factor the matrix. */
359: /* TODO figure out how to avoid reorder if the matrix values changed, but the pattern remains the same. */
360: PetscStackCall("STRUMPACK_reorder",sp_err = STRUMPACK_reorder(*S));
361: PetscStackCall("STRUMPACK_factor",sp_err = STRUMPACK_factor(*S));
362: switch (sp_err) {
363: case STRUMPACK_SUCCESS: break;
364: case STRUMPACK_MATRIX_NOT_SET: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix was not set"); break; }
365: case STRUMPACK_REORDERING_ERROR: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix reordering failed"); break; }
366: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: factorization failed");
367: }
368: F->assembled = PETSC_TRUE;
369: F->preallocated = PETSC_TRUE;
370: return 0;
371: }
373: static PetscErrorCode MatLUFactorSymbolic_STRUMPACK(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
374: {
375: F->ops->lufactornumeric = MatLUFactorNumeric_STRUMPACK;
376: F->ops->solve = MatSolve_STRUMPACK;
377: F->ops->matsolve = MatMatSolve_STRUMPACK;
378: return 0;
379: }
381: static PetscErrorCode MatFactorGetSolverType_aij_strumpack(Mat A,MatSolverType *type)
382: {
383: *type = MATSOLVERSTRUMPACK;
384: return 0;
385: }
387: /*MC
388: MATSOLVERSSTRUMPACK = "strumpack" - A solver package providing a direct sparse solver (PCLU)
389: and a preconditioner (PCILU) using low-rank compression via the external package STRUMPACK.
391: Consult the STRUMPACK-sparse manual for more info.
393: Use
394: ./configure --download-strumpack
395: to have PETSc installed with STRUMPACK
397: Use
398: -pc_type lu -pc_factor_mat_solver_type strumpack
399: to use this as an exact (direct) solver, use
400: -pc_type ilu -pc_factor_mat_solver_type strumpack
401: to enable low-rank compression (i.e, use as a preconditioner).
403: Works with AIJ matrices
405: Options Database Keys:
406: + -mat_strumpack_verbose - verbose info
407: . -mat_strumpack_hss_rel_tol <1e-2> - Relative compression tolerance (None)
408: . -mat_strumpack_hss_abs_tol <1e-8> - Absolute compression tolerance (None)
409: . -mat_strumpack_colperm <TRUE> - Permute matrix to make diagonal nonzeros (None)
410: . -mat_strumpack_hss_min_sep_size <256> - Minimum size of separator for HSS compression (None)
411: . -mat_strumpack_max_rank - Maximum rank in HSS compression, when using pctype ilu (None)
412: . -mat_strumpack_leaf_size - Size of diagonal blocks in HSS approximation, when using pctype ilu (None)
413: . -mat_strumpack_reordering <METIS> - Sparsity reducing matrix reordering (choose one of) NATURAL METIS PARMETIS SCOTCH PTSCOTCH RCM (None)
414: - -mat_strumpack_iterative_solver <DIRECT> - Select iterative solver from STRUMPACK (choose one of) AUTO DIRECT REFINE PREC_GMRES GMRES PREC_BICGSTAB BICGSTAB (None)
416: Level: beginner
418: .seealso: PCLU, PCILU, MATSOLVERSUPERLU_DIST, MATSOLVERMUMPS, PCFactorSetMatSolverType(), MatSolverType
419: M*/
420: static PetscErrorCode MatGetFactor_aij_strumpack(Mat A,MatFactorType ftype,Mat *F)
421: {
422: Mat B;
423: PetscInt M=A->rmap->N,N=A->cmap->N;
424: PetscBool verb,flg,set;
425: PetscReal ctol;
426: PetscInt hssminsize,max_rank,leaf_size;
427: STRUMPACK_SparseSolver *S;
428: STRUMPACK_INTERFACE iface;
429: STRUMPACK_REORDERING_STRATEGY ndcurrent,ndvalue;
430: STRUMPACK_KRYLOV_SOLVER itcurrent,itsolver;
431: const STRUMPACK_PRECISION table[2][2][2] =
432: {{{STRUMPACK_FLOATCOMPLEX_64, STRUMPACK_DOUBLECOMPLEX_64},
433: {STRUMPACK_FLOAT_64, STRUMPACK_DOUBLE_64}},
434: {{STRUMPACK_FLOATCOMPLEX, STRUMPACK_DOUBLECOMPLEX},
435: {STRUMPACK_FLOAT, STRUMPACK_DOUBLE}}};
436: const STRUMPACK_PRECISION prec = table[(sizeof(PetscInt)==8)?0:1][(PETSC_SCALAR==PETSC_COMPLEX)?0:1][(PETSC_REAL==PETSC_FLOAT)?0:1];
437: const char *const STRUMPACKNDTypes[] = {"NATURAL","METIS","PARMETIS","SCOTCH","PTSCOTCH","RCM","STRUMPACKNDTypes","",0};
438: const char *const SolverTypes[] = {"AUTO","NONE","REFINE","PREC_GMRES","GMRES","PREC_BICGSTAB","BICGSTAB","SolverTypes","",0};
439: PetscErrorCode ierr;
441: /* Create the factorization matrix */
442: MatCreate(PetscObjectComm((PetscObject)A),&B);
443: MatSetSizes(B,A->rmap->n,A->cmap->n,M,N);
444: MatSetType(B,((PetscObject)A)->type_name);
445: MatSeqAIJSetPreallocation(B,0,NULL);
446: MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
447: B->trivialsymbolic = PETSC_TRUE;
448: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
449: B->ops->lufactorsymbolic = MatLUFactorSymbolic_STRUMPACK;
450: B->ops->ilufactorsymbolic = MatLUFactorSymbolic_STRUMPACK;
451: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
452: B->ops->view = MatView_STRUMPACK;
453: B->ops->destroy = MatDestroy_STRUMPACK;
454: B->ops->getdiagonal = MatGetDiagonal_STRUMPACK;
455: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_aij_strumpack);
456: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetReordering_C",MatSTRUMPACKSetReordering_STRUMPACK);
457: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetColPerm_C",MatSTRUMPACKSetColPerm_STRUMPACK);
458: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSRelTol_C",MatSTRUMPACKSetHSSRelTol_STRUMPACK);
459: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSAbsTol_C",MatSTRUMPACKSetHSSAbsTol_STRUMPACK);
460: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSMaxRank_C",MatSTRUMPACKSetHSSMaxRank_STRUMPACK);
461: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSLeafSize_C",MatSTRUMPACKSetHSSLeafSize_STRUMPACK);
462: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSMinSepSize_C",MatSTRUMPACKSetHSSMinSepSize_STRUMPACK);
463: B->factortype = ftype;
464: PetscNewLog(B,&S);
465: B->spptr = S;
467: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
468: iface = flg ? STRUMPACK_MT : STRUMPACK_MPI_DIST;
470: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"STRUMPACK Options","Mat");
472: verb = PetscLogPrintInfo ? PETSC_TRUE : PETSC_FALSE;
473: PetscOptionsBool("-mat_strumpack_verbose","Print STRUMPACK information","None",verb,&verb,NULL);
475: PetscStackCall("STRUMPACK_init",STRUMPACK_init(S,PetscObjectComm((PetscObject)A),prec,iface,0,NULL,verb));
477: PetscStackCall("STRUMPACK_HSS_rel_tol",ctol = (PetscReal)STRUMPACK_HSS_rel_tol(*S));
478: PetscOptionsReal("-mat_strumpack_hss_rel_tol","Relative compression tolerance","None",ctol,&ctol,&set);
479: if (set) PetscStackCall("STRUMPACK_set_HSS_rel_tol",STRUMPACK_set_HSS_rel_tol(*S,(double)ctol));
481: PetscStackCall("STRUMPACK_HSS_abs_tol",ctol = (PetscReal)STRUMPACK_HSS_abs_tol(*S));
482: PetscOptionsReal("-mat_strumpack_hss_abs_tol","Absolute compression tolerance","None",ctol,&ctol,&set);
483: if (set) PetscStackCall("STRUMPACK_set_HSS_abs_tol",STRUMPACK_set_HSS_abs_tol(*S,(double)ctol));
485: PetscStackCall("STRUMPACK_mc64job",flg = (STRUMPACK_mc64job(*S) == 0) ? PETSC_FALSE : PETSC_TRUE);
486: PetscOptionsBool("-mat_strumpack_colperm","Find a col perm to get nonzero diagonal","None",flg,&flg,&set);
487: if (set) PetscStackCall("STRUMPACK_set_mc64job",STRUMPACK_set_mc64job(*S,flg ? 5 : 0));
489: PetscStackCall("STRUMPACK_HSS_min_sep_size",hssminsize = (PetscInt)STRUMPACK_HSS_min_sep_size(*S));
490: PetscOptionsInt("-mat_strumpack_hss_min_sep_size","Minimum size of separator for HSS compression","None",hssminsize,&hssminsize,&set);
491: if (set) PetscStackCall("STRUMPACK_set_HSS_min_sep_size",STRUMPACK_set_HSS_min_sep_size(*S,(int)hssminsize));
493: PetscStackCall("STRUMPACK_HSS_max_rank",max_rank = (PetscInt)STRUMPACK_HSS_max_rank(*S));
494: PetscOptionsInt("-mat_strumpack_max_rank","Maximum rank in HSS approximation","None",max_rank,&max_rank,&set);
495: if (set) PetscStackCall("STRUMPACK_set_HSS_max_rank",STRUMPACK_set_HSS_max_rank(*S,(int)max_rank));
497: PetscStackCall("STRUMPACK_HSS_leaf_size",leaf_size = (PetscInt)STRUMPACK_HSS_leaf_size(*S));
498: PetscOptionsInt("-mat_strumpack_leaf_size","Size of diagonal blocks in HSS approximation","None",leaf_size,&leaf_size,&set);
499: if (set) PetscStackCall("STRUMPACK_set_HSS_leaf_size",STRUMPACK_set_HSS_leaf_size(*S,(int)leaf_size));
501: PetscStackCall("STRUMPACK_reordering_method",ndcurrent = STRUMPACK_reordering_method(*S));
502: PetscOptionsEnum("-mat_strumpack_reordering","Sparsity reducing matrix reordering","None",STRUMPACKNDTypes,(PetscEnum)ndcurrent,(PetscEnum*)&ndvalue,&set);
503: if (set) PetscStackCall("STRUMPACK_set_reordering_method",STRUMPACK_set_reordering_method(*S,ndvalue));
505: if (ftype == MAT_FACTOR_ILU) {
506: /* When enabling HSS compression, the STRUMPACK solver becomes an incomplete */
507: /* (or approximate) LU factorization. */
508: PetscStackCall("STRUMPACK_enable_HSS",STRUMPACK_enable_HSS(*S));
509: }
511: /* Disable the outer iterative solver from STRUMPACK. */
512: /* When STRUMPACK is used as a direct solver, it will by default do iterative refinement. */
513: /* When STRUMPACK is used as an approximate factorization preconditioner (by enabling */
514: /* low-rank compression), it will use it's own preconditioned GMRES. Here we can disable */
515: /* the outer iterative solver, as PETSc uses STRUMPACK from within a KSP. */
516: PetscStackCall("STRUMPACK_set_Krylov_solver", STRUMPACK_set_Krylov_solver(*S, STRUMPACK_DIRECT));
518: PetscStackCall("STRUMPACK_Krylov_solver",itcurrent = STRUMPACK_Krylov_solver(*S));
519: PetscOptionsEnum("-mat_strumpack_iterative_solver","Select iterative solver from STRUMPACK","None",SolverTypes,(PetscEnum)itcurrent,(PetscEnum*)&itsolver,&set);
520: if (set) PetscStackCall("STRUMPACK_set_Krylov_solver",STRUMPACK_set_Krylov_solver(*S,itsolver));
522: PetscOptionsEnd();
524: *F = B;
525: return 0;
526: }
528: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_STRUMPACK(void)
529: {
530: MatSolverTypeRegister(MATSOLVERSTRUMPACK,MATMPIAIJ,MAT_FACTOR_LU,MatGetFactor_aij_strumpack);
531: MatSolverTypeRegister(MATSOLVERSTRUMPACK,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_aij_strumpack);
532: MatSolverTypeRegister(MATSOLVERSTRUMPACK,MATMPIAIJ,MAT_FACTOR_ILU,MatGetFactor_aij_strumpack);
533: MatSolverTypeRegister(MATSOLVERSTRUMPACK,MATSEQAIJ,MAT_FACTOR_ILU,MatGetFactor_aij_strumpack);
534: return 0;
535: }