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