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: }
42: static PetscErrorCode MatSTRUMPACKSetReordering_STRUMPACK(Mat F,MatSTRUMPACKReordering reordering)
43: {
44: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
47: PetscStackCall("STRUMPACK_reordering_method",STRUMPACK_set_reordering_method(*S,(STRUMPACK_REORDERING_STRATEGY)reordering));
48: return(0);
49: }
51: /*@
52: MatSTRUMPACKSetReordering - Set STRUMPACK fill-reducing reordering
54: Input Parameters:
55: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
56: - reordering - the code to be used to find the fill-reducing reordering
57: Possible values: NATURAL=0 METIS=1 PARMETIS=2 SCOTCH=3 PTSCOTCH=4 RCM=5
59: Options Database:
60: . -mat_strumpack_reordering <METIS> - Sparsity reducing matrix reordering (choose one of) NATURAL METIS PARMETIS SCOTCH PTSCOTCH RCM (None)
62: Level: beginner
64: References:
65: . STRUMPACK manual
67: .seealso: MatGetFactor()
68: @*/
69: PetscErrorCode MatSTRUMPACKSetReordering(Mat F,MatSTRUMPACKReordering reordering)
70: {
76: PetscTryMethod(F,"MatSTRUMPACKSetReordering_C",(Mat,MatSTRUMPACKReordering),(F,reordering));
77: return(0);
78: }
80: static PetscErrorCode MatSTRUMPACKSetColPerm_STRUMPACK(Mat F,PetscBool cperm)
81: {
82: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
85: PetscStackCall("STRUMPACK_set_mc64job",STRUMPACK_set_mc64job(*S,cperm ? 5 : 0));
86: return(0);
87: }
89: /*@
90: MatSTRUMPACKSetColPerm - Set whether STRUMPACK should try to permute the columns of the matrix in order to get a nonzero diagonal
91: Logically Collective on Mat
93: Input Parameters:
94: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
95: - cperm - whether or not to permute (internally) the columns of the matrix
97: Options Database:
98: . -mat_strumpack_colperm <cperm>
100: Level: beginner
102: References:
103: . STRUMPACK manual
105: .seealso: MatGetFactor()
106: @*/
107: PetscErrorCode MatSTRUMPACKSetColPerm(Mat F,PetscBool cperm)
108: {
114: PetscTryMethod(F,"MatSTRUMPACKSetColPerm_C",(Mat,PetscBool),(F,cperm));
115: return(0);
116: }
118: static PetscErrorCode MatSTRUMPACKSetHSSRelTol_STRUMPACK(Mat F,PetscReal rtol)
119: {
120: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
123: PetscStackCall("STRUMPACK_set_HSS_rel_tol", STRUMPACK_set_HSS_rel_tol(*S,rtol));
124: return(0);
125: }
127: /*@
128: MatSTRUMPACKSetHSSRelTol - Set STRUMPACK relative tolerance for HSS compression
129: Logically Collective on Mat
131: Input Parameters:
132: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
133: - rtol - relative compression tolerance
135: Options Database:
136: . -mat_strumpack_hss_rel_tol <1e-2> - Relative compression tolerance (None)
138: Level: beginner
140: References:
141: . STRUMPACK manual
143: .seealso: MatGetFactor()
144: @*/
145: PetscErrorCode MatSTRUMPACKSetHSSRelTol(Mat F,PetscReal rtol)
146: {
152: PetscTryMethod(F,"MatSTRUMPACKSetHSSRelTol_C",(Mat,PetscReal),(F,rtol));
153: return(0);
154: }
156: static PetscErrorCode MatSTRUMPACKSetHSSAbsTol_STRUMPACK(Mat F,PetscReal atol)
157: {
158: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
161: PetscStackCall("STRUMPACK_set_HSS_abs_tol", STRUMPACK_set_HSS_abs_tol(*S,atol));
162: return(0);
163: }
165: /*@
166: MatSTRUMPACKSetHSSAbsTol - Set STRUMPACK absolute tolerance for HSS compression
167: Logically Collective on Mat
169: Input Parameters:
170: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
171: - atol - absolute compression tolerance
173: Options Database:
174: . -mat_strumpack_hss_abs_tol <1e-8> - Absolute compression tolerance (None)
176: Level: beginner
178: References:
179: . STRUMPACK manual
181: .seealso: MatGetFactor()
182: @*/
183: PetscErrorCode MatSTRUMPACKSetHSSAbsTol(Mat F,PetscReal atol)
184: {
190: PetscTryMethod(F,"MatSTRUMPACKSetHSSAbsTol_C",(Mat,PetscReal),(F,atol));
191: return(0);
192: }
194: static PetscErrorCode MatSTRUMPACKSetHSSMaxRank_STRUMPACK(Mat F,PetscInt hssmaxrank)
195: {
196: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
199: PetscStackCall("STRUMPACK_set_HSS_max_rank", STRUMPACK_set_HSS_max_rank(*S,hssmaxrank));
200: return(0);
201: }
203: /*@
204: MatSTRUMPACKSetHSSMaxRank - Set STRUMPACK maximum HSS rank
205: Logically Collective on Mat
207: Input Parameters:
208: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
209: - hssmaxrank - maximum rank used in low-rank approximation
211: Options Database:
212: . -mat_strumpack_max_rank - Maximum rank in HSS compression, when using pctype ilu (None)
214: Level: beginner
216: References:
217: . STRUMPACK manual
219: .seealso: MatGetFactor()
220: @*/
221: PetscErrorCode MatSTRUMPACKSetHSSMaxRank(Mat F,PetscInt hssmaxrank)
222: {
228: PetscTryMethod(F,"MatSTRUMPACKSetHSSMaxRank_C",(Mat,PetscInt),(F,hssmaxrank));
229: return(0);
230: }
232: static PetscErrorCode MatSTRUMPACKSetHSSLeafSize_STRUMPACK(Mat F,PetscInt leaf_size)
233: {
234: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
237: PetscStackCall("STRUMPACK_set_HSS_leaf_size", STRUMPACK_set_HSS_leaf_size(*S,leaf_size));
238: return(0);
239: }
241: /*@
242: MatSTRUMPACKSetHSSLeafSize - Set STRUMPACK HSS leaf size
243: Logically Collective on Mat
245: Input Parameters:
246: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
247: - leaf_size - Size of diagonal blocks in HSS approximation
249: Options Database:
250: . -mat_strumpack_leaf_size - Size of diagonal blocks in HSS approximation, when using pctype ilu (None)
252: Level: beginner
254: References:
255: . STRUMPACK manual
257: .seealso: MatGetFactor()
258: @*/
259: PetscErrorCode MatSTRUMPACKSetHSSLeafSize(Mat F,PetscInt leaf_size)
260: {
266: PetscTryMethod(F,"MatSTRUMPACKSetHSSLeafSize_C",(Mat,PetscInt),(F,leaf_size));
267: return(0);
268: }
270: static PetscErrorCode MatSTRUMPACKSetHSSMinSepSize_STRUMPACK(Mat F,PetscInt hssminsize)
271: {
272: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
275: PetscStackCall("STRUMPACK_set_HSS_min_sep_size", STRUMPACK_set_HSS_min_sep_size(*S,hssminsize));
276: return(0);
277: }
279: /*@
280: MatSTRUMPACKSetHSSMinSepSize - Set STRUMPACK minimum separator size for low-rank approximation
281: Logically Collective on Mat
283: Input Parameters:
284: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
285: - hssminsize - minimum dense matrix size for low-rank approximation
287: Options Database:
288: . -mat_strumpack_hss_min_sep_size <hssminsize>
290: Level: beginner
292: References:
293: . STRUMPACK manual
295: .seealso: MatGetFactor()
296: @*/
297: PetscErrorCode MatSTRUMPACKSetHSSMinSepSize(Mat F,PetscInt hssminsize)
298: {
304: PetscTryMethod(F,"MatSTRUMPACKSetHSSMinSepSize_C",(Mat,PetscInt),(F,hssminsize));
305: return(0);
306: }
308: static PetscErrorCode MatSolve_STRUMPACK(Mat A,Vec b_mpi,Vec x)
309: {
310: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)A->spptr;
311: STRUMPACK_RETURN_CODE sp_err;
312: PetscErrorCode ierr;
313: const PetscScalar *bptr;
314: PetscScalar *xptr;
317: VecGetArray(x,&xptr);
318: VecGetArrayRead(b_mpi,&bptr);
320: PetscStackCall("STRUMPACK_solve",sp_err = STRUMPACK_solve(*S,(PetscScalar*)bptr,xptr,0));
321: switch (sp_err) {
322: case STRUMPACK_SUCCESS: break;
323: case STRUMPACK_MATRIX_NOT_SET: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix was not set"); break; }
324: case STRUMPACK_REORDERING_ERROR: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix reordering failed"); break; }
325: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: solve failed");
326: }
327: VecRestoreArray(x,&xptr);
328: VecRestoreArrayRead(b_mpi,&bptr);
329: return(0);
330: }
332: static PetscErrorCode MatMatSolve_STRUMPACK(Mat A,Mat B_mpi,Mat X)
333: {
334: PetscErrorCode ierr;
335: PetscBool flg;
338: PetscObjectTypeCompareAny((PetscObject)B_mpi,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
339: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
340: PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
341: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
342: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatSolve_STRUMPACK() is not implemented yet");
343: return(0);
344: }
346: static PetscErrorCode MatView_Info_STRUMPACK(Mat A,PetscViewer viewer)
347: {
348: PetscErrorCode ierr;
351: /* check if matrix is strumpack type */
352: if (A->ops->solve != MatSolve_STRUMPACK) return(0);
353: PetscViewerASCIIPrintf(viewer,"STRUMPACK sparse solver!\n");
354: return(0);
355: }
357: static PetscErrorCode MatView_STRUMPACK(Mat A,PetscViewer viewer)
358: {
359: PetscErrorCode ierr;
360: PetscBool iascii;
361: PetscViewerFormat format;
364: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
365: if (iascii) {
366: PetscViewerGetFormat(viewer,&format);
367: if (format == PETSC_VIEWER_ASCII_INFO) {
368: MatView_Info_STRUMPACK(A,viewer);
369: }
370: }
371: return(0);
372: }
374: static PetscErrorCode MatLUFactorNumeric_STRUMPACK(Mat F,Mat A,const MatFactorInfo *info)
375: {
376: STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
377: STRUMPACK_RETURN_CODE sp_err;
378: Mat_SeqAIJ *A_d,*A_o;
379: Mat_MPIAIJ *mat;
380: PetscErrorCode ierr;
381: PetscInt M=A->rmap->N,m=A->rmap->n;
382: PetscBool flg;
385: PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);
386: if (flg) { /* A is MATMPIAIJ */
387: mat = (Mat_MPIAIJ*)A->data;
388: A_d = (Mat_SeqAIJ*)(mat->A)->data;
389: A_o = (Mat_SeqAIJ*)(mat->B)->data;
390: 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));
391: } else { /* A is MATSEQAIJ */
392: A_d = (Mat_SeqAIJ*)A->data;
393: PetscStackCall("STRUMPACK_set_csr_matrix",STRUMPACK_set_csr_matrix(*S,&M,A_d->i,A_d->j,A_d->a,0));
394: }
396: /* Reorder and Factor the matrix. */
397: /* TODO figure out how to avoid reorder if the matrix values changed, but the pattern remains the same. */
398: PetscStackCall("STRUMPACK_reorder",sp_err = STRUMPACK_reorder(*S));
399: PetscStackCall("STRUMPACK_factor",sp_err = STRUMPACK_factor(*S));
400: switch (sp_err) {
401: case STRUMPACK_SUCCESS: break;
402: case STRUMPACK_MATRIX_NOT_SET: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix was not set"); break; }
403: case STRUMPACK_REORDERING_ERROR: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix reordering failed"); break; }
404: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: factorization failed");
405: }
406: F->assembled = PETSC_TRUE;
407: F->preallocated = PETSC_TRUE;
408: return(0);
409: }
411: static PetscErrorCode MatLUFactorSymbolic_STRUMPACK(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
412: {
414: F->ops->lufactornumeric = MatLUFactorNumeric_STRUMPACK;
415: F->ops->solve = MatSolve_STRUMPACK;
416: F->ops->matsolve = MatMatSolve_STRUMPACK;
417: return(0);
418: }
420: static PetscErrorCode MatFactorGetSolverType_aij_strumpack(Mat A,MatSolverType *type)
421: {
423: *type = MATSOLVERSTRUMPACK;
424: return(0);
425: }
427: /*MC
428: MATSOLVERSSTRUMPACK = "strumpack" - A solver package providing a direct sparse solver (PCLU)
429: and a preconditioner (PCILU) using low-rank compression via the external package STRUMPACK.
431: Consult the STRUMPACK-sparse manual for more info.
433: Use
434: ./configure --download-strumpack
435: to have PETSc installed with STRUMPACK
437: Use
438: -pc_type lu -pc_factor_mat_solver_type strumpack
439: to use this as an exact (direct) solver, use
440: -pc_type ilu -pc_factor_mat_solver_type strumpack
441: to enable low-rank compression (i.e, use as a preconditioner).
443: Works with AIJ matrices
445: Options Database Keys:
446: + -mat_strumpack_verbose
447: . -mat_strumpack_hss_rel_tol <1e-2> - Relative compression tolerance (None)
448: . -mat_strumpack_hss_abs_tol <1e-8> - Absolute compression tolerance (None)
449: . -mat_strumpack_colperm <TRUE> - Permute matrix to make diagonal nonzeros (None)
450: . -mat_strumpack_hss_min_sep_size <256> - Minimum size of separator for HSS compression (None)
451: . -mat_strumpack_max_rank - Maximum rank in HSS compression, when using pctype ilu (None)
452: . -mat_strumpack_leaf_size - Size of diagonal blocks in HSS approximation, when using pctype ilu (None)
453: . -mat_strumpack_reordering <METIS> - Sparsity reducing matrix reordering (choose one of) NATURAL METIS PARMETIS SCOTCH PTSCOTCH RCM (None)
454: - -mat_strumpack_iterative_solver <DIRECT> - Select iterative solver from STRUMPACK (choose one of) AUTO DIRECT REFINE PREC_GMRES GMRES PREC_BICGSTAB BICGSTAB (None)
456: Level: beginner
458: .seealso: PCLU, PCILU, MATSOLVERSUPERLU_DIST, MATSOLVERMUMPS, PCFactorSetMatSolverType(), MatSolverType
459: M*/
460: static PetscErrorCode MatGetFactor_aij_strumpack(Mat A,MatFactorType ftype,Mat *F)
461: {
462: Mat B;
463: PetscErrorCode ierr;
464: PetscInt M=A->rmap->N,N=A->cmap->N;
465: PetscBool verb,flg,set;
466: PetscReal ctol;
467: PetscInt hssminsize,max_rank,leaf_size;
468: STRUMPACK_SparseSolver *S;
469: STRUMPACK_INTERFACE iface;
470: STRUMPACK_REORDERING_STRATEGY ndcurrent,ndvalue;
471: STRUMPACK_KRYLOV_SOLVER itcurrent,itsolver;
472: const STRUMPACK_PRECISION table[2][2][2] =
473: {{{STRUMPACK_FLOATCOMPLEX_64, STRUMPACK_DOUBLECOMPLEX_64},
474: {STRUMPACK_FLOAT_64, STRUMPACK_DOUBLE_64}},
475: {{STRUMPACK_FLOATCOMPLEX, STRUMPACK_DOUBLECOMPLEX},
476: {STRUMPACK_FLOAT, STRUMPACK_DOUBLE}}};
477: const STRUMPACK_PRECISION prec = table[(sizeof(PetscInt)==8)?0:1][(PETSC_SCALAR==PETSC_COMPLEX)?0:1][(PETSC_REAL==PETSC_FLOAT)?0:1];
478: const char *const STRUMPACKNDTypes[] = {"NATURAL","METIS","PARMETIS","SCOTCH","PTSCOTCH","RCM","STRUMPACKNDTypes","",0};
479: const char *const SolverTypes[] = {"AUTO","NONE","REFINE","PREC_GMRES","GMRES","PREC_BICGSTAB","BICGSTAB","SolverTypes","",0};
482: /* Create the factorization matrix */
483: MatCreate(PetscObjectComm((PetscObject)A),&B);
484: MatSetSizes(B,A->rmap->n,A->cmap->n,M,N);
485: MatSetType(B,((PetscObject)A)->type_name);
486: MatSeqAIJSetPreallocation(B,0,NULL);
487: MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
488: B->trivialsymbolic = PETSC_TRUE;
489: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
490: B->ops->lufactorsymbolic = MatLUFactorSymbolic_STRUMPACK;
491: B->ops->ilufactorsymbolic = MatLUFactorSymbolic_STRUMPACK;
492: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
493: B->ops->view = MatView_STRUMPACK;
494: B->ops->destroy = MatDestroy_STRUMPACK;
495: B->ops->getdiagonal = MatGetDiagonal_STRUMPACK;
496: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_aij_strumpack);
497: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetReordering_C",MatSTRUMPACKSetReordering_STRUMPACK);
498: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetColPerm_C",MatSTRUMPACKSetColPerm_STRUMPACK);
499: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSRelTol_C",MatSTRUMPACKSetHSSRelTol_STRUMPACK);
500: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSAbsTol_C",MatSTRUMPACKSetHSSAbsTol_STRUMPACK);
501: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSMaxRank_C",MatSTRUMPACKSetHSSMaxRank_STRUMPACK);
502: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSLeafSize_C",MatSTRUMPACKSetHSSLeafSize_STRUMPACK);
503: PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSMinSepSize_C",MatSTRUMPACKSetHSSMinSepSize_STRUMPACK);
504: B->factortype = ftype;
505: PetscNewLog(B,&S);
506: B->spptr = S;
508: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
509: iface = flg ? STRUMPACK_MT : STRUMPACK_MPI_DIST;
511: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"STRUMPACK Options","Mat");
513: verb = PetscLogPrintInfo ? PETSC_TRUE : PETSC_FALSE;
514: PetscOptionsBool("-mat_strumpack_verbose","Print STRUMPACK information","None",verb,&verb,NULL);
516: PetscStackCall("STRUMPACK_init",STRUMPACK_init(S,PetscObjectComm((PetscObject)A),prec,iface,0,NULL,verb));
518: PetscStackCall("STRUMPACK_HSS_rel_tol",ctol = (PetscReal)STRUMPACK_HSS_rel_tol(*S));
519: PetscOptionsReal("-mat_strumpack_hss_rel_tol","Relative compression tolerance","None",ctol,&ctol,&set);
520: if (set) PetscStackCall("STRUMPACK_set_HSS_rel_tol",STRUMPACK_set_HSS_rel_tol(*S,(double)ctol));
522: PetscStackCall("STRUMPACK_HSS_abs_tol",ctol = (PetscReal)STRUMPACK_HSS_abs_tol(*S));
523: PetscOptionsReal("-mat_strumpack_hss_abs_tol","Absolute compression tolerance","None",ctol,&ctol,&set);
524: if (set) PetscStackCall("STRUMPACK_set_HSS_abs_tol",STRUMPACK_set_HSS_abs_tol(*S,(double)ctol));
526: PetscStackCall("STRUMPACK_mc64job",flg = (STRUMPACK_mc64job(*S) == 0) ? PETSC_FALSE : PETSC_TRUE);
527: PetscOptionsBool("-mat_strumpack_colperm","Find a col perm to get nonzero diagonal","None",flg,&flg,&set);
528: if (set) PetscStackCall("STRUMPACK_set_mc64job",STRUMPACK_set_mc64job(*S,flg ? 5 : 0));
530: PetscStackCall("STRUMPACK_HSS_min_sep_size",hssminsize = (PetscInt)STRUMPACK_HSS_min_sep_size(*S));
531: PetscOptionsInt("-mat_strumpack_hss_min_sep_size","Minimum size of separator for HSS compression","None",hssminsize,&hssminsize,&set);
532: if (set) PetscStackCall("STRUMPACK_set_HSS_min_sep_size",STRUMPACK_set_HSS_min_sep_size(*S,(int)hssminsize));
534: PetscStackCall("STRUMPACK_HSS_max_rank",max_rank = (PetscInt)STRUMPACK_HSS_max_rank(*S));
535: PetscOptionsInt("-mat_strumpack_max_rank","Maximum rank in HSS approximation","None",max_rank,&max_rank,&set);
536: if (set) PetscStackCall("STRUMPACK_set_HSS_max_rank",STRUMPACK_set_HSS_max_rank(*S,(int)max_rank));
538: PetscStackCall("STRUMPACK_HSS_leaf_size",leaf_size = (PetscInt)STRUMPACK_HSS_leaf_size(*S));
539: PetscOptionsInt("-mat_strumpack_leaf_size","Size of diagonal blocks in HSS approximation","None",leaf_size,&leaf_size,&set);
540: if (set) PetscStackCall("STRUMPACK_set_HSS_leaf_size",STRUMPACK_set_HSS_leaf_size(*S,(int)leaf_size));
542: PetscStackCall("STRUMPACK_reordering_method",ndcurrent = STRUMPACK_reordering_method(*S));
543: PetscOptionsEnum("-mat_strumpack_reordering","Sparsity reducing matrix reordering","None",STRUMPACKNDTypes,(PetscEnum)ndcurrent,(PetscEnum*)&ndvalue,&set);
544: if (set) PetscStackCall("STRUMPACK_set_reordering_method",STRUMPACK_set_reordering_method(*S,ndvalue));
546: if (ftype == MAT_FACTOR_ILU) {
547: /* When enabling HSS compression, the STRUMPACK solver becomes an incomplete */
548: /* (or approximate) LU factorization. */
549: PetscStackCall("STRUMPACK_enable_HSS",STRUMPACK_enable_HSS(*S));
550: }
552: /* Disable the outer iterative solver from STRUMPACK. */
553: /* When STRUMPACK is used as a direct solver, it will by default do iterative refinement. */
554: /* When STRUMPACK is used as an approximate factorization preconditioner (by enabling */
555: /* low-rank compression), it will use it's own preconditioned GMRES. Here we can disable */
556: /* the outer iterative solver, as PETSc uses STRUMPACK from within a KSP. */
557: PetscStackCall("STRUMPACK_set_Krylov_solver", STRUMPACK_set_Krylov_solver(*S, STRUMPACK_DIRECT));
559: PetscStackCall("STRUMPACK_Krylov_solver",itcurrent = STRUMPACK_Krylov_solver(*S));
560: PetscOptionsEnum("-mat_strumpack_iterative_solver","Select iterative solver from STRUMPACK","None",SolverTypes,(PetscEnum)itcurrent,(PetscEnum*)&itsolver,&set);
561: if (set) PetscStackCall("STRUMPACK_set_Krylov_solver",STRUMPACK_set_Krylov_solver(*S,itsolver));
563: PetscOptionsEnd();
565: *F = B;
566: return(0);
567: }
569: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_STRUMPACK(void)
570: {
574: MatSolverTypeRegister(MATSOLVERSTRUMPACK,MATMPIAIJ,MAT_FACTOR_LU,MatGetFactor_aij_strumpack);
575: MatSolverTypeRegister(MATSOLVERSTRUMPACK,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_aij_strumpack);
576: MatSolverTypeRegister(MATSOLVERSTRUMPACK,MATMPIAIJ,MAT_FACTOR_ILU,MatGetFactor_aij_strumpack);
577: MatSolverTypeRegister(MATSOLVERSTRUMPACK,MATSEQAIJ,MAT_FACTOR_ILU,MatGetFactor_aij_strumpack);
578: return(0);
579: }