Actual source code: strumpack.c

petsc-master 2019-09-23
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>

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