Actual source code: strumpack.c

petsc-3.9.4 2018-09-11
Report Typos and Errors
  1:  #include <../src/mat/impls/aij/seq/aij.h>
  2:  #include <../src/mat/impls/aij/mpi/mpiaij.h>
  3: #include <StrumpackSparseSolver.h>

  7: static PetscErrorCode MatGetDiagonal_STRUMPACK(Mat A,Vec v)
  8: {
 10:   SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Mat type: STRUMPACK factor");
 11:   return(0);
 12: }

 16: static PetscErrorCode MatDestroy_STRUMPACK(Mat A)
 17: {
 18:   STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)A->spptr;
 19:   PetscErrorCode         ierr;
 20:   PetscBool              flg;

 23:   /* Deallocate STRUMPACK storage */
 24:   PetscStackCall("STRUMPACK_destroy",STRUMPACK_destroy(S));
 25:   PetscFree(A->spptr);
 26:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
 27:   if (flg) {
 28:     MatDestroy_SeqAIJ(A);
 29:   } else {
 30:     MatDestroy_MPIAIJ(A);
 31:   }

 33:   /* clear composed functions */
 34:   PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL);
 35:   PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetReordering_C",NULL);
 36:   PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetColPerm_C",NULL);
 37:   PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSRelTol_C",NULL);
 38:   PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSAbsTol_C",NULL);
 39:   PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSMaxRank_C",NULL);
 40:   PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSLeafSize_C",NULL);
 41:   PetscObjectComposeFunction((PetscObject)A,"MatSTRUMPACKSetHSSMinSepSize_C",NULL);

 43:   return(0);
 44: }


 49: static PetscErrorCode MatSTRUMPACKSetReordering_STRUMPACK(Mat F,MatSTRUMPACKReordering reordering)
 50: {
 51:   STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;

 54:   PetscStackCall("STRUMPACK_reordering_method",STRUMPACK_set_reordering_method(*S,(STRUMPACK_REORDERING_STRATEGY)reordering));
 55:   return(0);
 56: }
 59: /*@
 60:   MatSTRUMPACKSetReordering - Set STRUMPACK fill-reducing reordering

 62:    Input Parameters:
 63: +  F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
 64: -  reordering - the code to be used to find the fill-reducing reordering
 65:       Possible values: NATURAL=0 METIS=1 PARMETIS=2 SCOTCH=3 PTSCOTCH=4 RCM=5

 67:   Options Database:
 68: .   -mat_strumpack_reordering <METIS>  - Sparsity reducing matrix reordering (choose one of) NATURAL METIS PARMETIS SCOTCH PTSCOTCH RCM (None)

 70:    Level: beginner

 72:    References:
 73: .      STRUMPACK manual

 75: .seealso: MatGetFactor()
 76: @*/
 77: PetscErrorCode MatSTRUMPACKSetReordering(Mat F,MatSTRUMPACKReordering reordering)
 78: {

 84:   PetscTryMethod(F,"MatSTRUMPACKSetReordering_C",(Mat,MatSTRUMPACKReordering),(F,reordering));
 85:   return(0);
 86: }


 91: static PetscErrorCode MatSTRUMPACKSetColPerm_STRUMPACK(Mat F,PetscBool cperm)
 92: {
 93:   STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;

 96:   PetscStackCall("STRUMPACK_set_mc64job",STRUMPACK_set_mc64job(*S,cperm ? 5 : 0));
 97:   return(0);
 98: }
101: /*@
102:   MatSTRUMPACKSetColPerm - Set whether STRUMPACK should try to permute the columns of the matrix in order to get a nonzero diagonal
103:    Logically Collective on Mat

105:    Input Parameters:
106: +  F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
107: -  cperm - whether or not to permute (internally) the columns of the matrix

109:   Options Database:
110: .   -mat_strumpack_colperm <cperm>

112:    Level: beginner

114:    References:
115: .      STRUMPACK manual

117: .seealso: MatGetFactor()
118: @*/
119: PetscErrorCode MatSTRUMPACKSetColPerm(Mat F,PetscBool cperm)
120: {

126:   PetscTryMethod(F,"MatSTRUMPACKSetColPerm_C",(Mat,PetscBool),(F,cperm));
127:   return(0);
128: }

132: static PetscErrorCode MatSTRUMPACKSetHSSRelTol_STRUMPACK(Mat F,PetscReal rtol)
133: {
134:   STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;

137:   PetscStackCall("STRUMPACK_set_HSS_rel_tol", STRUMPACK_set_HSS_rel_tol(*S,rtol));
138:   return(0);
139: }
142: /*@
143:   MatSTRUMPACKSetHSSRelTol - Set STRUMPACK relative tolerance for HSS compression
144:    Logically Collective on Mat

146:    Input Parameters:
147: +  F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
148: -  rtol - relative compression tolerance

150:   Options Database:
151: .   -mat_strumpack_hss_rel_tol <1e-2>         - Relative compression tolerance (None)

153:    Level: beginner

155:    References:
156: .      STRUMPACK manual

158: .seealso: MatGetFactor()
159: @*/
160: PetscErrorCode MatSTRUMPACKSetHSSRelTol(Mat F,PetscReal rtol)
161: {

167:   PetscTryMethod(F,"MatSTRUMPACKSetHSSRelTol_C",(Mat,PetscReal),(F,rtol));
168:   return(0);
169: }


174: static PetscErrorCode MatSTRUMPACKSetHSSAbsTol_STRUMPACK(Mat F,PetscReal atol)
175: {
176:   STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;

179:   PetscStackCall("STRUMPACK_set_HSS_abs_tol", STRUMPACK_set_HSS_abs_tol(*S,atol));
180:   return(0);
181: }
184: /*@
185:   MatSTRUMPACKSetHSSAbsTol - Set STRUMPACK absolute tolerance for HSS compression
186:    Logically Collective on Mat

188:    Input Parameters:
189: +  F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
190: -  atol - absolute compression tolerance

192:   Options Database:
193: .   -mat_strumpack_hss_abs_tol <1e-8>         - Absolute compression tolerance (None)

195:    Level: beginner

197:    References:
198: .      STRUMPACK manual

200: .seealso: MatGetFactor()
201: @*/
202: PetscErrorCode MatSTRUMPACKSetHSSAbsTol(Mat F,PetscReal atol)
203: {

209:   PetscTryMethod(F,"MatSTRUMPACKSetHSSAbsTol_C",(Mat,PetscReal),(F,atol));
210:   return(0);
211: }


216: static PetscErrorCode MatSTRUMPACKSetHSSMaxRank_STRUMPACK(Mat F,PetscInt hssmaxrank)
217: {
218:   STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;

221:   PetscStackCall("STRUMPACK_set_HSS_max_rank", STRUMPACK_set_HSS_max_rank(*S,hssmaxrank));
222:   return(0);
223: }
226: /*@
227:   MatSTRUMPACKSetHSSMaxRank - Set STRUMPACK maximum HSS rank
228:    Logically Collective on Mat

230:    Input Parameters:
231: +  F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
232: -  hssmaxrank - maximum rank used in low-rank approximation

234:   Options Database:
235: .   -mat_strumpack_max_rank    - Maximum rank in HSS compression, when using pctype ilu (None)

237:    Level: beginner

239:    References:
240: .      STRUMPACK manual

242: .seealso: MatGetFactor()
243: @*/
244: PetscErrorCode MatSTRUMPACKSetHSSMaxRank(Mat F,PetscInt hssmaxrank)
245: {

251:   PetscTryMethod(F,"MatSTRUMPACKSetHSSMaxRank_C",(Mat,PetscInt),(F,hssmaxrank));
252:   return(0);
253: }


258: static PetscErrorCode MatSTRUMPACKSetHSSLeafSize_STRUMPACK(Mat F,PetscInt leaf_size)
259: {
260:   STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;

263:   PetscStackCall("STRUMPACK_set_HSS_leaf_size", STRUMPACK_set_HSS_leaf_size(*S,leaf_size));
264:   return(0);
265: }
268: /*@
269:   MatSTRUMPACKSetHSSLeafSize - Set STRUMPACK HSS leaf size
270:    Logically Collective on Mat

272:    Input Parameters:
273: +  F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
274: -  leaf_size - Size of diagonal blocks in HSS approximation

276:   Options Database:
277: .   -mat_strumpack_leaf_size    - Size of diagonal blocks in HSS approximation, when using pctype ilu (None)

279:    Level: beginner

281:    References:
282: .      STRUMPACK manual

284: .seealso: MatGetFactor()
285: @*/
286: PetscErrorCode MatSTRUMPACKSetHSSLeafSize(Mat F,PetscInt leaf_size)
287: {

293:   PetscTryMethod(F,"MatSTRUMPACKSetHSSLeafSize_C",(Mat,PetscInt),(F,leaf_size));
294:   return(0);
295: }



301: static PetscErrorCode MatSTRUMPACKSetHSSMinSepSize_STRUMPACK(Mat F,PetscInt hssminsize)
302: {
303:   STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;

306:   PetscStackCall("STRUMPACK_set_HSS_min_sep_size", STRUMPACK_set_HSS_min_sep_size(*S,hssminsize));
307:   return(0);
308: }
311: /*@
312:   MatSTRUMPACKSetHSSMinSepSize - Set STRUMPACK minimum separator size for low-rank approximation
313:    Logically Collective on Mat

315:    Input Parameters:
316: +  F - the factored matrix obtained by calling MatGetFactor() from PETSc-STRUMPACK interface
317: -  hssminsize - minimum dense matrix size for low-rank approximation

319:   Options Database:
320: .   -mat_strumpack_hss_min_sep_size <hssminsize>

322:    Level: beginner

324:    References:
325: .      STRUMPACK manual

327: .seealso: MatGetFactor()
328: @*/
329: PetscErrorCode MatSTRUMPACKSetHSSMinSepSize(Mat F,PetscInt hssminsize)
330: {

336:   PetscTryMethod(F,"MatSTRUMPACKSetHSSMinSepSize_C",(Mat,PetscInt),(F,hssminsize));
337:   return(0);
338: }


343: static PetscErrorCode MatSolve_STRUMPACK(Mat A,Vec b_mpi,Vec x)
344: {
345:   STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)A->spptr;
346:   STRUMPACK_RETURN_CODE  sp_err;
347:   PetscErrorCode         ierr;
348:   const PetscScalar      *bptr;
349:   PetscScalar            *xptr;

352:   VecGetArray(x,&xptr);
353:   VecGetArrayRead(b_mpi,&bptr);

355:   PetscStackCall("STRUMPACK_solve",sp_err = STRUMPACK_solve(*S,(PetscScalar*)bptr,xptr,0));
356:   switch (sp_err) {
357:   case STRUMPACK_SUCCESS: break;
358:   case STRUMPACK_MATRIX_NOT_SET:   { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix was not set"); break; }
359:   case STRUMPACK_REORDERING_ERROR: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix reordering failed"); break; }
360:   default:                           SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: solve failed");
361:   }
362:   VecRestoreArray(x,&xptr);
363:   VecRestoreArrayRead(b_mpi,&bptr);
364:   return(0);
365: }

369: static PetscErrorCode MatMatSolve_STRUMPACK(Mat A,Mat B_mpi,Mat X)
370: {
371:   PetscErrorCode   ierr;
372:   PetscBool        flg;

375:   PetscObjectTypeCompareAny((PetscObject)B_mpi,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
376:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
377:   PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
378:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
379:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatSolve_STRUMPACK() is not implemented yet");
380:   return(0);
381: }

385: static PetscErrorCode MatView_Info_STRUMPACK(Mat A,PetscViewer viewer)
386: {
387:   PetscErrorCode  ierr;

390:   /* check if matrix is strumpack type */
391:   if (A->ops->solve != MatSolve_STRUMPACK) return(0);
392:   PetscViewerASCIIPrintf(viewer,"STRUMPACK sparse solver!\n");
393:   return(0);
394: }

398: static PetscErrorCode MatView_STRUMPACK(Mat A,PetscViewer viewer)
399: {
400:   PetscErrorCode    ierr;
401:   PetscBool         iascii;
402:   PetscViewerFormat format;

405:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
406:   if (iascii) {
407:     PetscViewerGetFormat(viewer,&format);
408:     if (format == PETSC_VIEWER_ASCII_INFO) {
409:       MatView_Info_STRUMPACK(A,viewer);
410:     }
411:   }
412:   return(0);
413: }

417: static PetscErrorCode MatLUFactorNumeric_STRUMPACK(Mat F,Mat A,const MatFactorInfo *info)
418: {
419:   STRUMPACK_SparseSolver *S = (STRUMPACK_SparseSolver*)F->spptr;
420:   STRUMPACK_RETURN_CODE  sp_err;
421:   Mat_SeqAIJ             *A_d,*A_o;
422:   Mat_MPIAIJ             *mat;
423:   PetscErrorCode         ierr;
424:   PetscInt               M=A->rmap->N,m=A->rmap->n;
425:   PetscBool              flg;

428:   PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);
429:   if (flg) { /* A is MATMPIAIJ */
430:     mat = (Mat_MPIAIJ*)A->data;
431:     A_d = (Mat_SeqAIJ*)(mat->A)->data;
432:     A_o = (Mat_SeqAIJ*)(mat->B)->data;
433:     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));
434:   } else { /* A is MATSEQAIJ */
435:     A_d = (Mat_SeqAIJ*)A->data;
436:     PetscStackCall("STRUMPACK_set_csr_matrix",STRUMPACK_set_csr_matrix(*S,&M,A_d->i,A_d->j,A_d->a,0));
437:   }

439:   /* Reorder and Factor the matrix. */
440:   /* TODO figure out how to avoid reorder if the matrix values changed, but the pattern remains the same. */
441:   PetscStackCall("STRUMPACK_reorder",sp_err = STRUMPACK_reorder(*S));
442:   PetscStackCall("STRUMPACK_factor",sp_err = STRUMPACK_factor(*S));
443:   switch (sp_err) {
444:   case STRUMPACK_SUCCESS: break;
445:   case STRUMPACK_MATRIX_NOT_SET:   { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix was not set"); break; }
446:   case STRUMPACK_REORDERING_ERROR: { SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: matrix reordering failed"); break; }
447:   default:                           SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"STRUMPACK error: factorization failed");
448:   }
449:   F->assembled    = PETSC_TRUE;
450:   F->preallocated = PETSC_TRUE;
451:   return(0);
452: }

456: static PetscErrorCode MatLUFactorSymbolic_STRUMPACK(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
457: {
459:   F->ops->lufactornumeric = MatLUFactorNumeric_STRUMPACK;
460:   F->ops->solve           = MatSolve_STRUMPACK;
461:   F->ops->matsolve        = MatMatSolve_STRUMPACK;
462:   return(0);
463: }

467: static PetscErrorCode MatFactorGetSolverType_aij_strumpack(Mat A,MatSolverType *type)
468: {
470:   *type = MATSOLVERSTRUMPACK;
471:   return(0);
472: }

474: /*MC
475:   MATSOLVERSSTRUMPACK = "strumpack" - A solver package providing a direct sparse solver (PCLU)
476:   and a preconditioner (PCILU) using low-rank compression via the external package STRUMPACK.

478:   Consult the STRUMPACK-sparse manual for more info.

480:   Use
481:      ./configure --download-strumpack
482:   to have PETSc installed with STRUMPACK

484:   Use
485:     -pc_type lu -pc_factor_mat_solver_type strumpack
486:   to use this as an exact (direct) solver, use
487:     -pc_type ilu -pc_factor_mat_solver_type strumpack
488:   to enable low-rank compression (i.e, use as a preconditioner).

490:   Works with AIJ matrices

492:   Options Database Keys:
493: + -mat_strumpack_verbose
494: . -mat_strumpack_hss_rel_tol <1e-2>         - Relative compression tolerance (None)
495: . -mat_strumpack_hss_abs_tol <1e-8>         - Absolute compression tolerance (None)
496: . -mat_strumpack_colperm <TRUE>             - Permute matrix to make diagonal nonzeros (None)
497: . -mat_strumpack_hss_min_sep_size <256>     - Minimum size of separator for HSS compression (None)
498: . -mat_strumpack_max_rank                   - Maximum rank in HSS compression, when using pctype ilu (None)
499: . -mat_strumpack_leaf_size                  - Size of diagonal blocks in HSS approximation, when using pctype ilu (None)
500: . -mat_strumpack_reordering <METIS>         - Sparsity reducing matrix reordering (choose one of) NATURAL METIS PARMETIS SCOTCH PTSCOTCH RCM (None)
501: . -mat_strumpack_iterative_solver <DIRECT>  - Select iterative solver from STRUMPACK (choose one of) AUTO DIRECT REFINE PREC_GMRES GMRES PREC_BICGSTAB BICGSTAB (None)

503:  Level: beginner

505: .seealso: PCLU, PCILU, MATSOLVERSUPERLU_DIST, MATSOLVERMUMPS, PCFactorSetMatSolverType(), MatSolverType
506: M*/
509: static PetscErrorCode MatGetFactor_aij_strumpack(Mat A,MatFactorType ftype,Mat *F)
510: {
511:   Mat                           B;
512:   PetscErrorCode                ierr;
513:   PetscInt                      M=A->rmap->N,N=A->cmap->N;
514:   PetscBool                     verb,flg,set;
515:   PetscReal                     ctol;
516:   PetscInt                      hssminsize,max_rank,leaf_size;
517:   STRUMPACK_SparseSolver        *S;
518:   STRUMPACK_INTERFACE           iface;
519:   STRUMPACK_REORDERING_STRATEGY ndcurrent,ndvalue;
520:   STRUMPACK_KRYLOV_SOLVER       itcurrent,itsolver;
521:   const STRUMPACK_PRECISION table[2][2][2] =
522:     {{{STRUMPACK_FLOATCOMPLEX_64, STRUMPACK_DOUBLECOMPLEX_64},
523:       {STRUMPACK_FLOAT_64,        STRUMPACK_DOUBLE_64}},
524:      {{STRUMPACK_FLOATCOMPLEX,    STRUMPACK_DOUBLECOMPLEX},
525:       {STRUMPACK_FLOAT,           STRUMPACK_DOUBLE}}};
526:   const STRUMPACK_PRECISION prec =
527:     table[(sizeof(PetscInt)==8)?0:1]
528:     [(PETSC_SCALAR==PETSC_COMPLEX)?0:1]
529:     [(PETSC_REAL==PETSC_FLOAT)?0:1];
530:   const char *const STRUMPACKNDTypes[] =
531:     {"NATURAL","METIS","PARMETIS","SCOTCH","PTSCOTCH","RCM","STRUMPACKNDTypes","",0};
532:   const char *const SolverTypes[] =
533:     {"AUTO","NONE","REFINE","PREC_GMRES","GMRES","PREC_BICGSTAB","BICGSTAB","SolverTypes","",0};

536:   /* Create the factorization matrix */
537:   MatCreate(PetscObjectComm((PetscObject)A),&B);
538:   MatSetSizes(B,A->rmap->n,A->cmap->n,M,N);
539:   MatSetType(B,((PetscObject)A)->type_name);
540:   MatSeqAIJSetPreallocation(B,0,NULL);
541:   MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);
542:   if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
543:     B->ops->lufactorsymbolic  = MatLUFactorSymbolic_STRUMPACK;
544:     B->ops->ilufactorsymbolic = MatLUFactorSymbolic_STRUMPACK;
545:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
546:   B->ops->view        = MatView_STRUMPACK;
547:   B->ops->destroy     = MatDestroy_STRUMPACK;
548:   B->ops->getdiagonal = MatGetDiagonal_STRUMPACK;
549:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_aij_strumpack);
550:   PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetReordering_C",MatSTRUMPACKSetReordering_STRUMPACK);
551:   PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetColPerm_C",MatSTRUMPACKSetColPerm_STRUMPACK);
552:   PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSRelTol_C",MatSTRUMPACKSetHSSRelTol_STRUMPACK);
553:   PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSAbsTol_C",MatSTRUMPACKSetHSSAbsTol_STRUMPACK);
554:   PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSMaxRank_C",MatSTRUMPACKSetHSSMaxRank_STRUMPACK);
555:   PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSLeafSize_C",MatSTRUMPACKSetHSSLeafSize_STRUMPACK);
556:   PetscObjectComposeFunction((PetscObject)B,"MatSTRUMPACKSetHSSMinSepSize_C",MatSTRUMPACKSetHSSMinSepSize_STRUMPACK);
557:   B->factortype = ftype;
558:   PetscNewLog(B,&S);
559:   B->spptr = S;

561:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
562:   iface = flg ? STRUMPACK_MT : STRUMPACK_MPI_DIST;

564:   PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"STRUMPACK Options","Mat");

566:   verb = PetscLogPrintInfo ? PETSC_TRUE : PETSC_FALSE;
567:   PetscOptionsBool("-mat_strumpack_verbose","Print STRUMPACK information","None",verb,&verb,NULL);

569:   PetscStackCall("STRUMPACK_init",STRUMPACK_init(S,PetscObjectComm((PetscObject)A),prec,iface,0,NULL,verb));

571:   PetscStackCall("STRUMPACK_HSS_rel_tol",ctol = (PetscReal)STRUMPACK_HSS_rel_tol(*S));
572:   PetscOptionsReal("-mat_strumpack_hss_rel_tol","Relative compression tolerance","None",ctol,&ctol,&set);
573:   if (set) PetscStackCall("STRUMPACK_set_HSS_rel_tol",STRUMPACK_set_HSS_rel_tol(*S,(double)ctol));

575:   PetscStackCall("STRUMPACK_HSS_abs_tol",ctol = (PetscReal)STRUMPACK_HSS_abs_tol(*S));
576:   PetscOptionsReal("-mat_strumpack_hss_abs_tol","Absolute compression tolerance","None",ctol,&ctol,&set);
577:   if (set) PetscStackCall("STRUMPACK_set_HSS_abs_tol",STRUMPACK_set_HSS_abs_tol(*S,(double)ctol));

579:   PetscStackCall("STRUMPACK_mc64job",flg = (STRUMPACK_mc64job(*S) == 0) ? PETSC_FALSE : PETSC_TRUE);
580:   PetscOptionsBool("-mat_strumpack_colperm","Find a col perm to get nonzero diagonal","None",flg,&flg,&set);
581:   if (set) PetscStackCall("STRUMPACK_set_mc64job",STRUMPACK_set_mc64job(*S,flg ? 5 : 0));

583:   PetscStackCall("STRUMPACK_HSS_min_sep_size",hssminsize = (PetscInt)STRUMPACK_HSS_min_sep_size(*S));
584:   PetscOptionsInt("-mat_strumpack_hss_min_sep_size","Minimum size of separator for HSS compression","None",hssminsize,&hssminsize,&set);
585:   if (set) PetscStackCall("STRUMPACK_set_HSS_min_sep_size",STRUMPACK_set_HSS_min_sep_size(*S,(int)hssminsize));

587:   PetscStackCall("STRUMPACK_HSS_max_rank",max_rank = (PetscInt)STRUMPACK_HSS_max_rank(*S));
588:   PetscOptionsInt("-mat_strumpack_max_rank","Maximum rank in HSS approximation","None",max_rank,&max_rank,&set);
589:   if (set) PetscStackCall("STRUMPACK_set_HSS_max_rank",STRUMPACK_set_HSS_max_rank(*S,(int)max_rank));

591:   PetscStackCall("STRUMPACK_HSS_leaf_size",leaf_size = (PetscInt)STRUMPACK_HSS_leaf_size(*S));
592:   PetscOptionsInt("-mat_strumpack_leaf_size","Size of diagonal blocks in HSS approximation","None",leaf_size,&leaf_size,&set);
593:   if (set) PetscStackCall("STRUMPACK_set_HSS_leaf_size",STRUMPACK_set_HSS_leaf_size(*S,(int)leaf_size));

595:   PetscStackCall("STRUMPACK_reordering_method",ndcurrent = STRUMPACK_reordering_method(*S));
596:   PetscOptionsEnum("-mat_strumpack_reordering","Sparsity reducing matrix reordering","None",STRUMPACKNDTypes,(PetscEnum)ndcurrent,(PetscEnum*)&ndvalue,&set);
597:   if (set) PetscStackCall("STRUMPACK_set_reordering_method",STRUMPACK_set_reordering_method(*S,ndvalue));

599:   if (ftype == MAT_FACTOR_ILU) {
600:     /* When enabling HSS compression, the STRUMPACK solver becomes an incomplete              */
601:     /* (or approximate) LU factorization.                                                     */
602:     PetscStackCall("STRUMPACK_enable_HSS",STRUMPACK_enable_HSS(*S));
603:   }

605:   /* Disable the outer iterative solver from STRUMPACK.                                       */
606:   /* When STRUMPACK is used as a direct solver, it will by default do iterative refinement.   */
607:   /* When STRUMPACK is used as an approximate factorization preconditioner (by enabling       */
608:   /* low-rank compression), it will use it's own preconditioned GMRES. Here we can disable    */
609:   /* the outer iterative solver, as PETSc uses STRUMPACK from within a KSP.                   */
610:   PetscStackCall("STRUMPACK_set_Krylov_solver", STRUMPACK_set_Krylov_solver(*S, STRUMPACK_DIRECT));

612:   PetscStackCall("STRUMPACK_Krylov_solver",itcurrent = STRUMPACK_Krylov_solver(*S));
613:   PetscOptionsEnum("-mat_strumpack_iterative_solver","Select iterative solver from STRUMPACK","None",SolverTypes,(PetscEnum)itcurrent,(PetscEnum*)&itsolver,&set);
614:   if (set) PetscStackCall("STRUMPACK_set_Krylov_solver",STRUMPACK_set_Krylov_solver(*S,itsolver));

616:   PetscOptionsEnd();

618:   *F = B;
619:   return(0);
620: }

624: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_STRUMPACK(void)
625: {

629:   MatSolverTypeRegister(MATSOLVERSTRUMPACK,MATMPIAIJ,MAT_FACTOR_LU,MatGetFactor_aij_strumpack);
630:   MatSolverTypeRegister(MATSOLVERSTRUMPACK,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_aij_strumpack);
631:   MatSolverTypeRegister(MATSOLVERSTRUMPACK,MATMPIAIJ,MAT_FACTOR_ILU,MatGetFactor_aij_strumpack);
632:   MatSolverTypeRegister(MATSOLVERSTRUMPACK,MATSEQAIJ,MAT_FACTOR_ILU,MatGetFactor_aij_strumpack);
633:   return(0);
634: }