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: {
8: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Mat type: STRUMPACK factor");
9: return(0);
10: }
12: static PetscErrorCode MatDestroy_STRUMPACK(Mat A)
13: {
14: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)A->spptr;
15: PetscErrorCode ierr;
16: PetscBool flg;
19: /* Deallocate STRUMPACK storage */
20: PetscStackCall("STRUMPACK_destroy",STRUMPACK_destroy(S));
21: PetscFree(A->spptr);
22: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
23: if (flg) {
24: MatDestroy_SeqAIJ(A);
25: } else {
26: MatDestroy_MPIAIJ(A);
27: }
29: /* clear composed functions */
30: PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL);
31: PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetReordering_C",NULL);
32: PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetColPerm_C",NULL);
33: PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSRelTol_C",NULL);
34: PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSAbsTol_C",NULL);
35: PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSMaxRank_C",NULL);
36: PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSLeafSize_C",NULL);
37: PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSMinSepSize_C",NULL);
39: return(0);
40: }
43: static PetscErrorCode MatSTRUMPACKSetReordering_STRUMPACK(Mat F,MatSTRUMPACKReordering reordering)
44: {
45: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
48: PetscStackCall("STRUMPACK_reordering_method",STRUMPACK_set_reordering_method(*S,(STRUMPACK_REORDERING_STRATEGY)reordering));
49: return(0);
50: }
52: /*@
53: MatSTRUMPACKSetReordering - Set STRUMPACK fill-reducing reordering
55: Input Parameters:
56: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
57: - reordering - the code to be used to find the fill-reducing reordering
58: Possible values: NATURAL=0 METIS=1 PARMETIS=2 SCOTCH=3 PTSCOTCH=4 RCM=5
60: Options Database:
61: . -mat_strumpack_reordering <METIS> - Sparsity reducing matrix reordering (choose one of) NATURAL METIS PARMETIS SCOTCH PTSCOTCH RCM (None)
63: Level: beginner
65: References:
66: . STRUMPACK manual
68: .seealso: MatGetFactor()
69: @*/
70: PetscErrorCode MatSTRUMPACKSetReordering(Mat F,MatSTRUMPACKReordering reordering)
71: {
77: PetscTryMethod(F,"MatSTRUMPACKSetReordering_C",(Mat,MatSTRUMPACKReordering),(F,reordering));
78: return(0);
79: }
81: static PetscErrorCode MatSTRUMPACKSetColPerm_STRUMPACK(Mat F,PetscBool cperm)
82: {
83: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
86: PetscStackCall("STRUMPACK_set_mc64job",STRUMPACK_set_mc64job(*S,cperm ? 5 : 0));
87: return(0);
88: }
90: /*@
91: MatSTRUMPACKSetColPerm - Set whether STRUMPACK should try to permute the columns of the matrix in order to get a nonzero diagonal
92: Logically Collective on Mat
94: Input Parameters:
95: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
96: - cperm - whether or not to permute (internally) the columns of the matrix
98: Options Database:
99: . -mat_strumpack_colperm <cperm>
101: Level: beginner
103: References:
104: . STRUMPACK manual
106: .seealso: MatGetFactor()
107: @*/
108: PetscErrorCode MatSTRUMPACKSetColPerm(Mat F,PetscBool cperm)
109: {
115: PetscTryMethod(F,"MatSTRUMPACKSetColPerm_C",(Mat,PetscBool),(F,cperm));
116: return(0);
117: }
119: static PetscErrorCode MatSTRUMPACKSetHSSRelTol_STRUMPACK(Mat F,PetscReal rtol)
120: {
121: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
124: PetscStackCall("STRUMPACK_set_HSS_rel_tol", STRUMPACK_set_HSS_rel_tol(*S,rtol));
125: return(0);
126: }
128: /*@
129: MatSTRUMPACKSetHSSRelTol - Set STRUMPACK relative tolerance for HSS compression
130: Logically Collective on Mat
132: Input Parameters:
133: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
134: - rtol - relative compression tolerance
136: Options Database:
137: . -mat_strumpack_hss_rel_tol <1e-2> - Relative compression tolerance (None)
139: Level: beginner
141: References:
142: . STRUMPACK manual
144: .seealso: MatGetFactor()
145: @*/
146: PetscErrorCode MatSTRUMPACKSetHSSRelTol(Mat F,PetscReal rtol)
147: {
153: PetscTryMethod(F,"MatSTRUMPACKSetHSSRelTol_C",(Mat,PetscReal),(F,rtol));
154: return(0);
155: }
157: static PetscErrorCode MatSTRUMPACKSetHSSAbsTol_STRUMPACK(Mat F,PetscReal atol)
158: {
159: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
162: PetscStackCall("STRUMPACK_set_HSS_abs_tol", STRUMPACK_set_HSS_abs_tol(*S,atol));
163: return(0);
164: }
166: /*@
167: MatSTRUMPACKSetHSSAbsTol - Set STRUMPACK absolute tolerance for HSS compression
168: Logically Collective on Mat
170: Input Parameters:
171: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
172: - atol - absolute compression tolerance
174: Options Database:
175: . -mat_strumpack_hss_abs_tol <1e-8> - Absolute compression tolerance (None)
177: Level: beginner
179: References:
180: . STRUMPACK manual
182: .seealso: MatGetFactor()
183: @*/
184: PetscErrorCode MatSTRUMPACKSetHSSAbsTol(Mat F,PetscReal atol)
185: {
191: PetscTryMethod(F,"MatSTRUMPACKSetHSSAbsTol_C",(Mat,PetscReal),(F,atol));
192: return(0);
193: }
195: static PetscErrorCode MatSTRUMPACKSetHSSMaxRank_STRUMPACK(Mat F,PetscInt hssmaxrank)
196: {
197: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
200: PetscStackCall("STRUMPACK_set_HSS_max_rank", STRUMPACK_set_HSS_max_rank(*S,hssmaxrank));
201: return(0);
202: }
204: /*@
205: MatSTRUMPACKSetHSSMaxRank - Set STRUMPACK maximum HSS rank
206: Logically Collective on Mat
208: Input Parameters:
209: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
210: - hssmaxrank - maximum rank used in low-rank approximation
212: Options Database:
213: . -mat_strumpack_max_rank - Maximum rank in HSS compression, when using pctype ilu (None)
215: Level: beginner
217: References:
218: . STRUMPACK manual
220: .seealso: MatGetFactor()
221: @*/
222: PetscErrorCode MatSTRUMPACKSetHSSMaxRank(Mat F,PetscInt hssmaxrank)
223: {
229: PetscTryMethod(F,"MatSTRUMPACKSetHSSMaxRank_C",(Mat,PetscInt),(F,hssmaxrank));
230: return(0);
231: }
233: static PetscErrorCode MatSTRUMPACKSetHSSLeafSize_STRUMPACK(Mat F,PetscInt leaf_size)
234: {
235: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
238: PetscStackCall("STRUMPACK_set_HSS_leaf_size", STRUMPACK_set_HSS_leaf_size(*S,leaf_size));
239: return(0);
240: }
242: /*@
243: MatSTRUMPACKSetHSSLeafSize - Set STRUMPACK HSS leaf size
244: Logically Collective on Mat
246: Input Parameters:
247: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
248: - leaf_size - Size of diagonal blocks in HSS approximation
250: Options Database:
251: . -mat_strumpack_leaf_size - Size of diagonal blocks in HSS approximation, when using pctype ilu (None)
253: Level: beginner
255: References:
256: . STRUMPACK manual
258: .seealso: MatGetFactor()
259: @*/
260: PetscErrorCode MatSTRUMPACKSetHSSLeafSize(Mat F,PetscInt leaf_size)
261: {
267: PetscTryMethod(F,"MatSTRUMPACKSetHSSLeafSize_C",(Mat,PetscInt),(F,leaf_size));
268: return(0);
269: }
271: static PetscErrorCode MatSTRUMPACKSetHSSMinSepSize_STRUMPACK(Mat F,PetscInt hssminsize)
272: {
273: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
276: PetscStackCall("STRUMPACK_set_HSS_min_sep_size", STRUMPACK_set_HSS_min_sep_size(*S,hssminsize));
277: return(0);
278: }
280: /*@
281: MatSTRUMPACKSetHSSMinSepSize - Set STRUMPACK minimum separator size for low-rank approximation
282: Logically Collective on Mat
284: Input Parameters:
285: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
286: - hssminsize - minimum dense matrix size for low-rank approximation
288: Options Database:
289: . -mat_strumpack_hss_min_sep_size <hssminsize>
291: Level: beginner
293: References:
294: . STRUMPACK manual
296: .seealso: MatGetFactor()
297: @*/
298: PetscErrorCode MatSTRUMPACKSetHSSMinSepSize(Mat F,PetscInt hssminsize)
299: {
305: PetscTryMethod(F,"MatSTRUMPACKSetHSSMinSepSize_C",(Mat,PetscInt),(F,hssminsize));
306: return(0);
307: }
309: static PetscErrorCode MatSolve_STRUMPACK(Mat A,Vec b_mpi,Vec x)
310: {
311: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)A->spptr;
312: STRUMPACK_RETURN_CODE sp_err;
313: PetscErrorCode ierr;
314: const PetscScalar *bptr;
315: PetscScalar *xptr;
318: VecGetArray(x,&xptr);
319: VecGetArrayRead(b_mpi,&bptr);
321: PetscStackCall("STRUMPACK_solve",sp_err = STRUMPACK_solve(*S,(PetscScalar*)bptr,xptr,0));
322: switch (sp_err) {
323: case STRUMPACK_SUCCESS: break;
324: case STRUMPACK_MATRIX_NOT_SET: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix was not set"); break; }
325: case STRUMPACK_REORDERING_ERROR: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix reordering failed"); break; }
326: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: solve failed");
327: }
328: VecRestoreArray(x,&xptr);
329: VecRestoreArrayRead(b_mpi,&bptr);
330: return(0);
331: }
333: static PetscErrorCode MatMatSolve_STRUMPACK(Mat A,Mat B_mpi,Mat X)
334: {
335: PetscErrorCode ierr;
336: PetscBool flg;
339: PetscObjectTypeCompareAny((PetscObject)B_mpi,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
340: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
341: PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
342: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
343: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatSolve_STRUMPACK() is not implemented yet");
344: return(0);
345: }
347: static PetscErrorCode MatView_Info_STRUMPACK(Mat A,PetscViewer viewer)
348: {
349: PetscErrorCode ierr;
352: /* check if matrix is strumpack type */
353: if (A->ops->solve != MatSolve_STRUMPACK) return(0);
354: PetscViewerASCIIPrintf(viewer,"STRUMPACK sparse solver!\n");
355: return(0);
356: }
358: static PetscErrorCode MatView_STRUMPACK(Mat A,PetscViewer viewer)
359: {
360: PetscErrorCode ierr;
361: PetscBool iascii;
362: PetscViewerFormat format;
365: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
366: if (iascii) {
367: PetscViewerGetFormat(viewer,&format);
368: if (format == PETSC_VIEWER_ASCII_INFO) {
369: MatView_Info_STRUMPACK(A,viewer);
370: }
371: }
372: return(0);
373: }
375: static PetscErrorCode MatLUFactorNumeric_STRUMPACK(Mat F,Mat A,const MatFactorInfo *info)
376: {
377: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
378: STRUMPACK_RETURN_CODE sp_err;
379: Mat_SeqAIJ *A_d,*A_o;
380: Mat_MPIAIJ *mat;
381: PetscErrorCode ierr;
382: PetscInt M=A->rmap->N,m=A->rmap->n;
383: PetscBool flg;
386: PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);
387: if (flg) { /* A is MATMPIAIJ */
388: mat = (Mat_MPIAIJ*)A->data;
389: A_d = (Mat_SeqAIJ*)(mat->A)->data;
390: A_o = (Mat_SeqAIJ*)(mat->B)->data;
391: 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));
392: } else { /* A is MATSEQAIJ */
393: A_d = (Mat_SeqAIJ*)A->data;
394: PetscStackCall("STRUMPACK_set_csr_matrix",STRUMPACK_set_csr_matrix(*S,&M,A_d->i,A_d->j,A_d->a,0));
395: }
397: /* Reorder and Factor the matrix. */
398: /* TODO figure out how to avoid reorder if the matrix values changed, but the pattern remains the same. */
399: PetscStackCall("STRUMPACK_reorder",sp_err = STRUMPACK_reorder(*S));
400: PetscStackCall("STRUMPACK_factor",sp_err = STRUMPACK_factor(*S));
401: switch (sp_err) {
402: case STRUMPACK_SUCCESS: break;
403: case STRUMPACK_MATRIX_NOT_SET: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix was not set"); break; }
404: case STRUMPACK_REORDERING_ERROR: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix reordering failed"); break; }
405: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: factorization failed");
406: }
407: F->assembled = PETSC_TRUE;
408: F->preallocated = PETSC_TRUE;
409: return(0);
410: }
412: static PetscErrorCode MatLUFactorSymbolic_STRUMPACK(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
413: {
415: F->ops->lufactornumeric = MatLUFactorNumeric_STRUMPACK;
416: F->ops->solve = MatSolve_STRUMPACK;
417: F->ops->matsolve = MatMatSolve_STRUMPACK;
418: return(0);
419: }
421: static PetscErrorCode MatFactorGetSolverType_aij_strumpack(Mat A,MatSolverType *type)
422: {
424: *type = MATSOLVERSTRUMPACK;
425: return(0);
426: }
428: /*MC
429: MATSOLVERSSTRUMPACK = "strumpack" - A solver package providing a direct sparse solver (PCLU)
430: and a preconditioner (PCILU) using low-rank compression via the external package STRUMPACK.
432: Consult the STRUMPACK-sparse manual for more info.
434: Use
435: ./configure --download-strumpack
436: to have PETSc installed with STRUMPACK
438: Use
439: -pc_type lu -pc_factor_mat_solver_type strumpack
440: to use this as an exact (direct) solver, use
441: -pc_type ilu -pc_factor_mat_solver_type strumpack
442: to enable low-rank compression (i.e, use as a preconditioner).
444: Works with AIJ matrices
446: Options Database Keys:
447: + -mat_strumpack_verbose
448: . -mat_strumpack_hss_rel_tol <1e-2> - Relative compression tolerance (None)
449: . -mat_strumpack_hss_abs_tol <1e-8> - Absolute compression tolerance (None)
450: . -mat_strumpack_colperm <TRUE> - Permute matrix to make diagonal nonzeros (None)
451: . -mat_strumpack_hss_min_sep_size <256> - Minimum size of separator for HSS compression (None)
452: . -mat_strumpack_max_rank - Maximum rank in HSS compression, when using pctype ilu (None)
453: . -mat_strumpack_leaf_size - Size of diagonal blocks in HSS approximation, when using pctype ilu (None)
454: . -mat_strumpack_reordering <METIS> - Sparsity reducing matrix reordering (choose one of) NATURAL METIS PARMETIS SCOTCH PTSCOTCH RCM (None)
455: - -mat_strumpack_iterative_solver <DIRECT> - Select iterative solver from STRUMPACK (choose one of) AUTO DIRECT REFINE PREC_GMRES GMRES PREC_BICGSTAB BICGSTAB (None)
457: Level: beginner
459: .seealso: PCLU, PCILU, MATSOLVERSUPERLU_DIST, MATSOLVERMUMPS, PCFactorSetMatSolverType(), MatSolverType
460: M*/
461: static PetscErrorCode MatGetFactor_aij_strumpack(Mat A,MatFactorType ftype,Mat *F)
462: {
463: Mat B;
464: PetscErrorCode ierr;
465: PetscInt M=A->rmap->N,N=A->cmap->N;
466: PetscBool verb,flg,set;
467: PetscReal ctol;
468: PetscInt hssminsize,max_rank,leaf_size;
469: STRUMPACK_SparseSolver *S;
470: STRUMPACK_INTERFACE iface;
471: STRUMPACK_REORDERING_STRATEGY ndcurrent,ndvalue;
472: STRUMPACK_KRYLOV_SOLVER itcurrent,itsolver;
473: const STRUMPACK_PRECISION table[2][2][2] =
474: {{{STRUMPACK_FLOATCOMPLEX_64, STRUMPACK_DOUBLECOMPLEX_64},
475: {STRUMPACK_FLOAT_64, STRUMPACK_DOUBLE_64}},
476: {{STRUMPACK_FLOATCOMPLEX, STRUMPACK_DOUBLECOMPLEX},
477: {STRUMPACK_FLOAT, STRUMPACK_DOUBLE}}};
478: const STRUMPACK_PRECISION prec =
479: table[(sizeof(PetscInt)==8)?0:1]
480: [(PETSC_SCALAR==PETSC_COMPLEX)?0:1]
481: [(PETSC_REAL==PETSC_FLOAT)?0:1];
482: const char *const STRUMPACKNDTypes[] =
483: {"NATURAL","METIS","PARMETIS","SCOTCH","PTSCOTCH","RCM","STRUMPACKNDTypes","",0};
484: const char *const SolverTypes[] =
485: {"AUTO","NONE","REFINE","PREC_GMRES","GMRES","PREC_BICGSTAB","BICGSTAB","SolverTypes","",0};
488: /* Create the factorization matrix */
489: MatCreate(PetscObjectComm((PetscObject)A),&B);
490: MatSetSizes(B,A->rmap->n,A->cmap->n,M,N);
491: MatSetType(B,((PetscObject)A)->type_name);
492: MatSeqAIJSetPreallocation(B,0,NULL);
493: MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
494: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
495: B->ops->lufactorsymbolic = MatLUFactorSymbolic_STRUMPACK;
496: B->ops->ilufactorsymbolic = MatLUFactorSymbolic_STRUMPACK;
497: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
498: B->ops->view = MatView_STRUMPACK;
499: B->ops->destroy = MatDestroy_STRUMPACK;
500: B->ops->getdiagonal = MatGetDiagonal_STRUMPACK;
501: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_aij_strumpack);
502: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetReordering_C",MatSTRUMPACKSetReordering_STRUMPACK);
503: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetColPerm_C",MatSTRUMPACKSetColPerm_STRUMPACK);
504: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSRelTol_C",MatSTRUMPACKSetHSSRelTol_STRUMPACK);
505: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSAbsTol_C",MatSTRUMPACKSetHSSAbsTol_STRUMPACK);
506: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSMaxRank_C",MatSTRUMPACKSetHSSMaxRank_STRUMPACK);
507: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSLeafSize_C",MatSTRUMPACKSetHSSLeafSize_STRUMPACK);
508: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSMinSepSize_C",MatSTRUMPACKSetHSSMinSepSize_STRUMPACK);
509: B->factortype = ftype;
510: PetscNewLog(B,&S);
511: B->spptr = S;
513: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
514: iface = flg ? STRUMPACK_MT : STRUMPACK_MPI_DIST;
516: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"STRUMPACK Options","Mat");
518: verb = PetscLogPrintInfo ? PETSC_TRUE : PETSC_FALSE;
519: PetscOptionsBool("-mat_strumpack_verbose","Print STRUMPACK information","None",verb,&verb,NULL);
521: PetscStackCall("STRUMPACK_init",STRUMPACK_init(S,PetscObjectComm((PetscObject)A),prec,iface,0,NULL,verb));
523: PetscStackCall("STRUMPACK_HSS_rel_tol",ctol = (PetscReal)STRUMPACK_HSS_rel_tol(*S));
524: PetscOptionsReal("-mat_strumpack_hss_rel_tol","Relative compression tolerance","None",ctol,&ctol,&set);
525: if (set) PetscStackCall("STRUMPACK_set_HSS_rel_tol",STRUMPACK_set_HSS_rel_tol(*S,(double)ctol));
527: PetscStackCall("STRUMPACK_HSS_abs_tol",ctol = (PetscReal)STRUMPACK_HSS_abs_tol(*S));
528: PetscOptionsReal("-mat_strumpack_hss_abs_tol","Absolute compression tolerance","None",ctol,&ctol,&set);
529: if (set) PetscStackCall("STRUMPACK_set_HSS_abs_tol",STRUMPACK_set_HSS_abs_tol(*S,(double)ctol));
531: PetscStackCall("STRUMPACK_mc64job",flg = (STRUMPACK_mc64job(*S) == 0) ? PETSC_FALSE : PETSC_TRUE);
532: PetscOptionsBool("-mat_strumpack_colperm","Find a col perm to get nonzero diagonal","None",flg,&flg,&set);
533: if (set) PetscStackCall("STRUMPACK_set_mc64job",STRUMPACK_set_mc64job(*S,flg ? 5 : 0));
535: PetscStackCall("STRUMPACK_HSS_min_sep_size",hssminsize = (PetscInt)STRUMPACK_HSS_min_sep_size(*S));
536: PetscOptionsInt("-mat_strumpack_hss_min_sep_size","Minimum size of separator for HSS compression","None",hssminsize,&hssminsize,&set);
537: if (set) PetscStackCall("STRUMPACK_set_HSS_min_sep_size",STRUMPACK_set_HSS_min_sep_size(*S,(int)hssminsize));
539: PetscStackCall("STRUMPACK_HSS_max_rank",max_rank = (PetscInt)STRUMPACK_HSS_max_rank(*S));
540: PetscOptionsInt("-mat_strumpack_max_rank","Maximum rank in HSS approximation","None",max_rank,&max_rank,&set);
541: if (set) PetscStackCall("STRUMPACK_set_HSS_max_rank",STRUMPACK_set_HSS_max_rank(*S,(int)max_rank));
543: PetscStackCall("STRUMPACK_HSS_leaf_size",leaf_size = (PetscInt)STRUMPACK_HSS_leaf_size(*S));
544: PetscOptionsInt("-mat_strumpack_leaf_size","Size of diagonal blocks in HSS approximation","None",leaf_size,&leaf_size,&set);
545: if (set) PetscStackCall("STRUMPACK_set_HSS_leaf_size",STRUMPACK_set_HSS_leaf_size(*S,(int)leaf_size));
547: PetscStackCall("STRUMPACK_reordering_method",ndcurrent = STRUMPACK_reordering_method(*S));
548: PetscOptionsEnum("-mat_strumpack_reordering","Sparsity reducing matrix reordering","None",STRUMPACKNDTypes,(PetscEnum)ndcurrent,(PetscEnum*)&ndvalue,&set);
549: if (set) PetscStackCall("STRUMPACK_set_reordering_method",STRUMPACK_set_reordering_method(*S,ndvalue));
551: if (ftype == MAT_FACTOR_ILU) {
552: /* When enabling HSS compression, the STRUMPACK solver becomes an incomplete */
553: /* (or approximate) LU factorization. */
554: PetscStackCall("STRUMPACK_enable_HSS",STRUMPACK_enable_HSS(*S));
555: }
557: /* Disable the outer iterative solver from STRUMPACK. */
558: /* When STRUMPACK is used as a direct solver, it will by default do iterative refinement. */
559: /* When STRUMPACK is used as an approximate factorization preconditioner (by enabling */
560: /* low-rank compression), it will use it's own preconditioned GMRES. Here we can disable */
561: /* the outer iterative solver, as PETSc uses STRUMPACK from within a KSP. */
562: PetscStackCall("STRUMPACK_set_Krylov_solver", STRUMPACK_set_Krylov_solver(*S, STRUMPACK_DIRECT));
564: PetscStackCall("STRUMPACK_Krylov_solver",itcurrent = STRUMPACK_Krylov_solver(*S));
565: PetscOptionsEnum("-mat_strumpack_iterative_solver","Select iterative solver from STRUMPACK","None",SolverTypes,(PetscEnum)itcurrent,(PetscEnum*)&itsolver,&set);
566: if (set) PetscStackCall("STRUMPACK_set_Krylov_solver",STRUMPACK_set_Krylov_solver(*S,itsolver));
568: PetscOptionsEnd();
570: *F = B;
571: return(0);
572: }
574: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_STRUMPACK(void)
575: {
579: MatSolverTypeRegister(MATSOLVERSTRUMPACK,MATMPIAIJ,MAT_FACTOR_LU,MatGetFactor_aij_strumpack);
580: MatSolverTypeRegister(MATSOLVERSTRUMPACK,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_aij_strumpack);
581: MatSolverTypeRegister(MATSOLVERSTRUMPACK,MATMPIAIJ,MAT_FACTOR_ILU,MatGetFactor_aij_strumpack);
582: MatSolverTypeRegister(MATSOLVERSTRUMPACK,MATSEQAIJ,MAT_FACTOR_ILU,MatGetFactor_aij_strumpack);
583: return(0);
584: }