Actual source code: umfpack.c

      1: /*
      2:    Provides an interface to the UMFPACK sparse solver available through SuiteSparse version 4.2.1

      4:    When build with PETSC_USE_64BIT_INDICES this will use Suitesparse_long as the
      5:    integer type in UMFPACK, otherwise it will use int. This means
      6:    all integers in this file as simply declared as PetscInt. Also it means
      7:    that one cannot use 64BIT_INDICES on 32-bit pointer systems [as Suitesparse_long is 32-bit only]

      9: */
     10: #include <../src/mat/impls/aij/seq/aij.h>

     12: #if defined(PETSC_USE_64BIT_INDICES)
     13:   #if defined(PETSC_USE_COMPLEX)
     14:     #define umfpack_UMF_free_symbolic umfpack_zl_free_symbolic
     15:     #define umfpack_UMF_free_numeric  umfpack_zl_free_numeric
     16:     /* the type casts are needed because PetscInt is long long while SuiteSparse_long is long and compilers warn even when they are identical */
     17:     #define umfpack_UMF_wsolve(a, b, c, d, e, f, g, h, i, j, k, l, m, n) umfpack_zl_wsolve(a, (SuiteSparse_long *)b, (SuiteSparse_long *)c, d, e, f, g, h, i, (SuiteSparse_long *)j, k, l, (SuiteSparse_long *)m, n)
     18:     #define umfpack_UMF_numeric(a, b, c, d, e, f, g, h)                  umfpack_zl_numeric((SuiteSparse_long *)a, (SuiteSparse_long *)b, c, d, e, f, g, h)
     19:     #define umfpack_UMF_report_numeric                                   umfpack_zl_report_numeric
     20:     #define umfpack_UMF_report_control                                   umfpack_zl_report_control
     21:     #define umfpack_UMF_report_status                                    umfpack_zl_report_status
     22:     #define umfpack_UMF_report_info                                      umfpack_zl_report_info
     23:     #define umfpack_UMF_report_symbolic                                  umfpack_zl_report_symbolic
     24:     #define umfpack_UMF_qsymbolic(a, b, c, d, e, f, g, h, i, j)          umfpack_zl_qsymbolic(a, b, (SuiteSparse_long *)c, (SuiteSparse_long *)d, e, f, (SuiteSparse_long *)g, h, i, j)
     25:     #define umfpack_UMF_symbolic(a, b, c, d, e, f, g, h, i)              umfpack_zl_symbolic(a, b, (SuiteSparse_long *)c, (SuiteSparse_long *)d, e, f, g, h, i)
     26:     #define umfpack_UMF_defaults                                         umfpack_zl_defaults

     28:   #else
     29:     #define umfpack_UMF_free_symbolic                           umfpack_dl_free_symbolic
     30:     #define umfpack_UMF_free_numeric                            umfpack_dl_free_numeric
     31:     #define umfpack_UMF_wsolve(a, b, c, d, e, f, g, h, i, j, k) umfpack_dl_wsolve(a, (SuiteSparse_long *)b, (SuiteSparse_long *)c, d, e, f, g, h, i, (SuiteSparse_long *)j, k)
     32:     #define umfpack_UMF_numeric(a, b, c, d, e, f, g)            umfpack_dl_numeric((SuiteSparse_long *)a, (SuiteSparse_long *)b, c, d, e, f, g)
     33:     #define umfpack_UMF_report_numeric                          umfpack_dl_report_numeric
     34:     #define umfpack_UMF_report_control                          umfpack_dl_report_control
     35:     #define umfpack_UMF_report_status                           umfpack_dl_report_status
     36:     #define umfpack_UMF_report_info                             umfpack_dl_report_info
     37:     #define umfpack_UMF_report_symbolic                         umfpack_dl_report_symbolic
     38:     #define umfpack_UMF_qsymbolic(a, b, c, d, e, f, g, h, i)    umfpack_dl_qsymbolic(a, b, (SuiteSparse_long *)c, (SuiteSparse_long *)d, e, (SuiteSparse_long *)f, g, h, i)
     39:     #define umfpack_UMF_symbolic(a, b, c, d, e, f, g, h)        umfpack_dl_symbolic(a, b, (SuiteSparse_long *)c, (SuiteSparse_long *)d, e, f, g, h)
     40:     #define umfpack_UMF_defaults                                umfpack_dl_defaults
     41:   #endif

     43: #else
     44:   #if defined(PETSC_USE_COMPLEX)
     45:     #define umfpack_UMF_free_symbolic   umfpack_zi_free_symbolic
     46:     #define umfpack_UMF_free_numeric    umfpack_zi_free_numeric
     47:     #define umfpack_UMF_wsolve          umfpack_zi_wsolve
     48:     #define umfpack_UMF_numeric         umfpack_zi_numeric
     49:     #define umfpack_UMF_report_numeric  umfpack_zi_report_numeric
     50:     #define umfpack_UMF_report_control  umfpack_zi_report_control
     51:     #define umfpack_UMF_report_status   umfpack_zi_report_status
     52:     #define umfpack_UMF_report_info     umfpack_zi_report_info
     53:     #define umfpack_UMF_report_symbolic umfpack_zi_report_symbolic
     54:     #define umfpack_UMF_qsymbolic       umfpack_zi_qsymbolic
     55:     #define umfpack_UMF_symbolic        umfpack_zi_symbolic
     56:     #define umfpack_UMF_defaults        umfpack_zi_defaults

     58:   #else
     59:     #define umfpack_UMF_free_symbolic   umfpack_di_free_symbolic
     60:     #define umfpack_UMF_free_numeric    umfpack_di_free_numeric
     61:     #define umfpack_UMF_wsolve          umfpack_di_wsolve
     62:     #define umfpack_UMF_numeric         umfpack_di_numeric
     63:     #define umfpack_UMF_report_numeric  umfpack_di_report_numeric
     64:     #define umfpack_UMF_report_control  umfpack_di_report_control
     65:     #define umfpack_UMF_report_status   umfpack_di_report_status
     66:     #define umfpack_UMF_report_info     umfpack_di_report_info
     67:     #define umfpack_UMF_report_symbolic umfpack_di_report_symbolic
     68:     #define umfpack_UMF_qsymbolic       umfpack_di_qsymbolic
     69:     #define umfpack_UMF_symbolic        umfpack_di_symbolic
     70:     #define umfpack_UMF_defaults        umfpack_di_defaults
     71:   #endif
     72: #endif

     74: EXTERN_C_BEGIN
     75: #include <umfpack.h>
     76: EXTERN_C_END

     78: static const char *const UmfpackOrderingTypes[] = {"CHOLMOD", "AMD", "GIVEN", "METIS", "BEST", "NONE", "USER", "UmfpackOrderingTypes", "UMFPACK_ORDERING_", 0};

     80: typedef struct {
     81:   void        *Symbolic, *Numeric;
     82:   double       Info[UMFPACK_INFO], Control[UMFPACK_CONTROL], *W;
     83:   PetscInt    *Wi, *perm_c;
     84:   Mat          A; /* Matrix used for factorization */
     85:   MatStructure flg;

     87:   /* Flag to clean up UMFPACK objects during Destroy */
     88:   PetscBool CleanUpUMFPACK;
     89: } Mat_UMFPACK;

     91: static PetscErrorCode MatDestroy_UMFPACK(Mat A)
     92: {
     93:   Mat_UMFPACK *lu = (Mat_UMFPACK *)A->data;

     95:   PetscFunctionBegin;
     96:   if (lu->CleanUpUMFPACK) {
     97:     umfpack_UMF_free_symbolic(&lu->Symbolic);
     98:     umfpack_UMF_free_numeric(&lu->Numeric);
     99:     PetscCall(PetscFree(lu->Wi));
    100:     PetscCall(PetscFree(lu->W));
    101:     PetscCall(PetscFree(lu->perm_c));
    102:   }
    103:   PetscCall(MatDestroy(&lu->A));
    104:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL));
    105:   PetscCall(PetscFree(A->data));
    106:   PetscFunctionReturn(PETSC_SUCCESS);
    107: }

    109: static PetscErrorCode MatSolve_UMFPACK_Private(Mat A, Vec b, Vec x, int uflag)
    110: {
    111:   Mat_UMFPACK       *lu = (Mat_UMFPACK *)A->data;
    112:   Mat_SeqAIJ        *a  = (Mat_SeqAIJ *)lu->A->data;
    113:   PetscScalar       *av = a->a, *xa;
    114:   const PetscScalar *ba;
    115:   PetscInt          *ai = a->i, *aj = a->j;
    116:   int                status;
    117:   static PetscBool   cite = PETSC_FALSE;

    119:   PetscFunctionBegin;
    120:   if (!A->rmap->n) PetscFunctionReturn(PETSC_SUCCESS);
    121:   PetscCall(PetscCitationsRegister("@article{davis2004algorithm,\n  title={Algorithm 832: {UMFPACK} V4.3---An Unsymmetric-Pattern Multifrontal Method},\n  author={Davis, Timothy A},\n  journal={ACM Transactions on Mathematical Software (TOMS)},\n  "
    122:                                    "volume={30},\n  number={2},\n  pages={196--199},\n  year={2004},\n  publisher={ACM}\n}\n",
    123:                                    &cite));
    124:   /* solve Ax = b by umfpack_*_wsolve */

    126:   if (!lu->Wi) { /* first time, allocate working space for wsolve */
    127:     PetscCall(PetscMalloc1(A->rmap->n, &lu->Wi));
    128:     PetscCall(PetscMalloc1(5 * A->rmap->n, &lu->W));
    129:   }

    131:   PetscCall(VecGetArrayRead(b, &ba));
    132:   PetscCall(VecGetArray(x, &xa));
    133: #if defined(PETSC_USE_COMPLEX)
    134:   status = umfpack_UMF_wsolve(uflag, ai, aj, (PetscReal *)av, NULL, (PetscReal *)xa, NULL, (PetscReal *)ba, NULL, lu->Numeric, lu->Control, lu->Info, lu->Wi, lu->W);
    135: #else
    136:   status = umfpack_UMF_wsolve(uflag, ai, aj, av, xa, ba, lu->Numeric, lu->Control, lu->Info, lu->Wi, lu->W);
    137: #endif
    138:   umfpack_UMF_report_info(lu->Control, lu->Info);
    139:   if (status < 0) {
    140:     umfpack_UMF_report_status(lu->Control, status);
    141:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "umfpack_UMF_wsolve failed");
    142:   }

    144:   PetscCall(VecRestoreArrayRead(b, &ba));
    145:   PetscCall(VecRestoreArray(x, &xa));
    146:   PetscFunctionReturn(PETSC_SUCCESS);
    147: }

    149: static PetscErrorCode MatSolve_UMFPACK(Mat A, Vec b, Vec x)
    150: {
    151:   PetscFunctionBegin;
    152:   /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */
    153:   PetscCall(MatSolve_UMFPACK_Private(A, b, x, UMFPACK_Aat));
    154:   PetscFunctionReturn(PETSC_SUCCESS);
    155: }

    157: static PetscErrorCode MatSolveTranspose_UMFPACK(Mat A, Vec b, Vec x)
    158: {
    159:   PetscFunctionBegin;
    160:   /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */
    161:   PetscCall(MatSolve_UMFPACK_Private(A, b, x, UMFPACK_A));
    162:   PetscFunctionReturn(PETSC_SUCCESS);
    163: }

    165: static PetscErrorCode MatLUFactorNumeric_UMFPACK(Mat F, Mat A, const MatFactorInfo *info)
    166: {
    167:   Mat_UMFPACK *lu = (Mat_UMFPACK *)F->data;
    168:   Mat_SeqAIJ  *a  = (Mat_SeqAIJ *)A->data;
    169:   PetscInt    *ai = a->i, *aj = a->j;
    170:   int          status;
    171:   PetscScalar *av = a->a;

    173:   PetscFunctionBegin;
    174:   if (!A->rmap->n) PetscFunctionReturn(PETSC_SUCCESS);
    175:   /* numeric factorization of A' */

    177:   if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric) umfpack_UMF_free_numeric(&lu->Numeric);
    178: #if defined(PETSC_USE_COMPLEX)
    179:   status = umfpack_UMF_numeric(ai, aj, (double *)av, NULL, lu->Symbolic, &lu->Numeric, lu->Control, lu->Info);
    180: #else
    181:   status = umfpack_UMF_numeric(ai, aj, av, lu->Symbolic, &lu->Numeric, lu->Control, lu->Info);
    182: #endif
    183:   if (status < 0) {
    184:     umfpack_UMF_report_status(lu->Control, status);
    185:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "umfpack_UMF_numeric failed");
    186:   }
    187:   /* report numeric factorization of A' when Control[PRL] > 3 */
    188:   (void)umfpack_UMF_report_numeric(lu->Numeric, lu->Control);

    190:   PetscCall(PetscObjectReference((PetscObject)A));
    191:   PetscCall(MatDestroy(&lu->A));

    193:   lu->A                  = A;
    194:   lu->flg                = SAME_NONZERO_PATTERN;
    195:   lu->CleanUpUMFPACK     = PETSC_TRUE;
    196:   F->ops->solve          = MatSolve_UMFPACK;
    197:   F->ops->solvetranspose = MatSolveTranspose_UMFPACK;
    198:   PetscFunctionReturn(PETSC_SUCCESS);
    199: }

    201: static PetscErrorCode MatLUFactorSymbolic_UMFPACK(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
    202: {
    203:   Mat_SeqAIJ  *a  = (Mat_SeqAIJ *)A->data;
    204:   Mat_UMFPACK *lu = (Mat_UMFPACK *)F->data;
    205:   PetscInt     i, *ai = a->i, *aj = a->j, m = A->rmap->n, n = A->cmap->n, idx;
    206:   int          status;
    207: #if !defined(PETSC_USE_COMPLEX)
    208:   PetscScalar *av = a->a;
    209: #endif
    210:   const PetscInt *ra;
    211:   const char     *strategy[] = {"AUTO", "UNSYMMETRIC", "SYMMETRIC"};
    212:   const char     *scale[]    = {"NONE", "SUM", "MAX"};
    213:   PetscBool       flg;

    215:   PetscFunctionBegin;
    216:   F->ops->lufactornumeric = MatLUFactorNumeric_UMFPACK;
    217:   if (!n) PetscFunctionReturn(PETSC_SUCCESS);

    219:   /* Set options to F */
    220:   PetscOptionsBegin(PetscObjectComm((PetscObject)F), ((PetscObject)F)->prefix, "UMFPACK Options", "Mat");
    221:   /* Control parameters used by reporting routiones */
    222:   PetscCall(PetscOptionsReal("-mat_umfpack_prl", "Control[UMFPACK_PRL]", "None", lu->Control[UMFPACK_PRL], &lu->Control[UMFPACK_PRL], NULL));

    224:   /* Control parameters for symbolic factorization */
    225:   PetscCall(PetscOptionsEList("-mat_umfpack_strategy", "ordering and pivoting strategy", "None", strategy, 3, strategy[0], &idx, &flg));
    226:   if (flg) {
    227:     switch (idx) {
    228:     case 0:
    229:       lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_AUTO;
    230:       break;
    231:     case 1:
    232:       lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC;
    233:       break;
    234:     case 2:
    235:       lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_SYMMETRIC;
    236:       break;
    237:     }
    238:   }
    239:   PetscCall(PetscOptionsEList("-mat_umfpack_ordering", "Internal ordering method", "None", UmfpackOrderingTypes, PETSC_STATIC_ARRAY_LENGTH(UmfpackOrderingTypes), UmfpackOrderingTypes[(int)lu->Control[UMFPACK_ORDERING]], &idx, &flg));
    240:   if (flg) lu->Control[UMFPACK_ORDERING] = (int)idx;
    241:   PetscCall(PetscOptionsReal("-mat_umfpack_dense_col", "Control[UMFPACK_DENSE_COL]", "None", lu->Control[UMFPACK_DENSE_COL], &lu->Control[UMFPACK_DENSE_COL], NULL));
    242:   PetscCall(PetscOptionsReal("-mat_umfpack_dense_row", "Control[UMFPACK_DENSE_ROW]", "None", lu->Control[UMFPACK_DENSE_ROW], &lu->Control[UMFPACK_DENSE_ROW], NULL));
    243:   PetscCall(PetscOptionsReal("-mat_umfpack_amd_dense", "Control[UMFPACK_AMD_DENSE]", "None", lu->Control[UMFPACK_AMD_DENSE], &lu->Control[UMFPACK_AMD_DENSE], NULL));
    244:   PetscCall(PetscOptionsReal("-mat_umfpack_block_size", "Control[UMFPACK_BLOCK_SIZE]", "None", lu->Control[UMFPACK_BLOCK_SIZE], &lu->Control[UMFPACK_BLOCK_SIZE], NULL));
    245:   PetscCall(PetscOptionsReal("-mat_umfpack_fixq", "Control[UMFPACK_FIXQ]", "None", lu->Control[UMFPACK_FIXQ], &lu->Control[UMFPACK_FIXQ], NULL));
    246:   PetscCall(PetscOptionsReal("-mat_umfpack_aggressive", "Control[UMFPACK_AGGRESSIVE]", "None", lu->Control[UMFPACK_AGGRESSIVE], &lu->Control[UMFPACK_AGGRESSIVE], NULL));

    248:   /* Control parameters used by numeric factorization */
    249:   PetscCall(PetscOptionsReal("-mat_umfpack_pivot_tolerance", "Control[UMFPACK_PIVOT_TOLERANCE]", "None", lu->Control[UMFPACK_PIVOT_TOLERANCE], &lu->Control[UMFPACK_PIVOT_TOLERANCE], NULL));
    250:   PetscCall(PetscOptionsReal("-mat_umfpack_sym_pivot_tolerance", "Control[UMFPACK_SYM_PIVOT_TOLERANCE]", "None", lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE], &lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE], NULL));
    251:   PetscCall(PetscOptionsEList("-mat_umfpack_scale", "Control[UMFPACK_SCALE]", "None", scale, 3, scale[0], &idx, &flg));
    252:   if (flg) {
    253:     switch (idx) {
    254:     case 0:
    255:       lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_NONE;
    256:       break;
    257:     case 1:
    258:       lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_SUM;
    259:       break;
    260:     case 2:
    261:       lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_MAX;
    262:       break;
    263:     }
    264:   }
    265:   PetscCall(PetscOptionsReal("-mat_umfpack_alloc_init", "Control[UMFPACK_ALLOC_INIT]", "None", lu->Control[UMFPACK_ALLOC_INIT], &lu->Control[UMFPACK_ALLOC_INIT], NULL));
    266:   PetscCall(PetscOptionsReal("-mat_umfpack_front_alloc_init", "Control[UMFPACK_FRONT_ALLOC_INIT]", "None", lu->Control[UMFPACK_FRONT_ALLOC_INIT], &lu->Control[UMFPACK_ALLOC_INIT], NULL));
    267:   PetscCall(PetscOptionsReal("-mat_umfpack_droptol", "Control[UMFPACK_DROPTOL]", "None", lu->Control[UMFPACK_DROPTOL], &lu->Control[UMFPACK_DROPTOL], NULL));

    269:   /* Control parameters used by solve */
    270:   PetscCall(PetscOptionsReal("-mat_umfpack_irstep", "Control[UMFPACK_IRSTEP]", "None", lu->Control[UMFPACK_IRSTEP], &lu->Control[UMFPACK_IRSTEP], NULL));
    271:   PetscOptionsEnd();

    273:   if (r) {
    274:     PetscCall(ISGetIndices(r, &ra));
    275:     PetscCall(PetscMalloc1(m, &lu->perm_c));
    276:     /* we cannot simply memcpy on 64-bit archs */
    277:     for (i = 0; i < m; i++) lu->perm_c[i] = ra[i];
    278:     PetscCall(ISRestoreIndices(r, &ra));
    279:   }

    281:   /* print the control parameters */
    282:   if (lu->Control[UMFPACK_PRL] > 1) umfpack_UMF_report_control(lu->Control);

    284:   /* symbolic factorization of A' */
    285:   if (r) { /* use PETSc row ordering */
    286: #if !defined(PETSC_USE_COMPLEX)
    287:     status = umfpack_UMF_qsymbolic(n, m, ai, aj, av, lu->perm_c, &lu->Symbolic, lu->Control, lu->Info);
    288: #else
    289:     status = umfpack_UMF_qsymbolic(n, m, ai, aj, NULL, NULL, lu->perm_c, &lu->Symbolic, lu->Control, lu->Info);
    290: #endif
    291:   } else { /* use Umfpack col ordering */
    292: #if !defined(PETSC_USE_COMPLEX)
    293:     status = umfpack_UMF_symbolic(n, m, ai, aj, av, &lu->Symbolic, lu->Control, lu->Info);
    294: #else
    295:     status = umfpack_UMF_symbolic(n, m, ai, aj, NULL, NULL, &lu->Symbolic, lu->Control, lu->Info);
    296: #endif
    297:   }
    298:   if (status < 0) {
    299:     umfpack_UMF_report_info(lu->Control, lu->Info);
    300:     umfpack_UMF_report_status(lu->Control, status);
    301:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "umfpack_UMF_symbolic failed");
    302:   }
    303:   /* report sumbolic factorization of A' when Control[PRL] > 3 */
    304:   (void)umfpack_UMF_report_symbolic(lu->Symbolic, lu->Control);

    306:   lu->flg            = DIFFERENT_NONZERO_PATTERN;
    307:   lu->CleanUpUMFPACK = PETSC_TRUE;
    308:   PetscFunctionReturn(PETSC_SUCCESS);
    309: }

    311: static PetscErrorCode MatView_Info_UMFPACK(Mat A, PetscViewer viewer)
    312: {
    313:   Mat_UMFPACK *lu = (Mat_UMFPACK *)A->data;

    315:   PetscFunctionBegin;
    316:   /* check if matrix is UMFPACK type */
    317:   if (A->ops->solve != MatSolve_UMFPACK) PetscFunctionReturn(PETSC_SUCCESS);

    319:   PetscCall(PetscViewerASCIIPrintf(viewer, "UMFPACK run parameters:\n"));
    320:   /* Control parameters used by reporting routiones */
    321:   PetscCall(PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_PRL]: %g\n", lu->Control[UMFPACK_PRL]));

    323:   /* Control parameters used by symbolic factorization */
    324:   PetscCall(PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_STRATEGY]: %g\n", lu->Control[UMFPACK_STRATEGY]));
    325:   PetscCall(PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_DENSE_COL]: %g\n", lu->Control[UMFPACK_DENSE_COL]));
    326:   PetscCall(PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_DENSE_ROW]: %g\n", lu->Control[UMFPACK_DENSE_ROW]));
    327:   PetscCall(PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_AMD_DENSE]: %g\n", lu->Control[UMFPACK_AMD_DENSE]));
    328:   PetscCall(PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_BLOCK_SIZE]: %g\n", lu->Control[UMFPACK_BLOCK_SIZE]));
    329:   PetscCall(PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_FIXQ]: %g\n", lu->Control[UMFPACK_FIXQ]));
    330:   PetscCall(PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_AGGRESSIVE]: %g\n", lu->Control[UMFPACK_AGGRESSIVE]));

    332:   /* Control parameters used by numeric factorization */
    333:   PetscCall(PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_PIVOT_TOLERANCE]: %g\n", lu->Control[UMFPACK_PIVOT_TOLERANCE]));
    334:   PetscCall(PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_SYM_PIVOT_TOLERANCE]: %g\n", lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE]));
    335:   PetscCall(PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_SCALE]: %g\n", lu->Control[UMFPACK_SCALE]));
    336:   PetscCall(PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_ALLOC_INIT]: %g\n", lu->Control[UMFPACK_ALLOC_INIT]));
    337:   PetscCall(PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_DROPTOL]: %g\n", lu->Control[UMFPACK_DROPTOL]));

    339:   /* Control parameters used by solve */
    340:   PetscCall(PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_IRSTEP]: %g\n", lu->Control[UMFPACK_IRSTEP]));

    342:   /* mat ordering */
    343:   if (!lu->perm_c) PetscCall(PetscViewerASCIIPrintf(viewer, "  Control[UMFPACK_ORDERING]: %s (not using the PETSc ordering)\n", UmfpackOrderingTypes[(int)lu->Control[UMFPACK_ORDERING]]));
    344:   PetscFunctionReturn(PETSC_SUCCESS);
    345: }

    347: static PetscErrorCode MatView_UMFPACK(Mat A, PetscViewer viewer)
    348: {
    349:   PetscBool         isascii;
    350:   PetscViewerFormat format;

    352:   PetscFunctionBegin;
    353:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    354:   if (isascii) {
    355:     PetscCall(PetscViewerGetFormat(viewer, &format));
    356:     if (format == PETSC_VIEWER_ASCII_INFO) PetscCall(MatView_Info_UMFPACK(A, viewer));
    357:   }
    358:   PetscFunctionReturn(PETSC_SUCCESS);
    359: }

    361: static PetscErrorCode MatFactorGetSolverType_seqaij_umfpack(Mat A, MatSolverType *type)
    362: {
    363:   PetscFunctionBegin;
    364:   *type = MATSOLVERUMFPACK;
    365:   PetscFunctionReturn(PETSC_SUCCESS);
    366: }

    368: /*MC
    369:   MATSOLVERUMFPACK = "umfpack" - A matrix type providing direct solvers, LU, for sequential matrices
    370:   via the external package UMFPACK.

    372:   Use `./configure --download-suitesparse` to install PETSc to use UMFPACK

    374:   Use `-pc_type lu` `-pc_factor_mat_solver_type umfpack` to use this direct solver

    376:   Consult UMFPACK documentation for more information about the Control parameters
    377:   which correspond to the options database keys below.

    379:   Options Database Keys:
    380: + -mat_umfpack_ordering                - `CHOLMOD`, `AMD`, `GIVEN`, `METIS`, `BEST`, `NONE`
    381: . -mat_umfpack_prl                     - UMFPACK print level: Control[UMFPACK_PRL]
    382: . -mat_umfpack_strategy <AUTO>         - (choose one of) `AUTO`, `UNSYMMETRIC`, `SYMMETRIC 2BY2`
    383: . -mat_umfpack_dense_col <alpha_c>     - UMFPACK dense column threshold: Control[UMFPACK_DENSE_COL]
    384: . -mat_umfpack_dense_row <0.2>         - Control[UMFPACK_DENSE_ROW]
    385: . -mat_umfpack_amd_dense <10>          - Control[UMFPACK_AMD_DENSE]
    386: . -mat_umfpack_block_size <bs>         - UMFPACK block size for BLAS-Level 3 calls: Control[UMFPACK_BLOCK_SIZE]
    387: . -mat_umfpack_2by2_tolerance <0.01>   - Control[UMFPACK_2BY2_TOLERANCE]
    388: . -mat_umfpack_fixq <0>                - Control[UMFPACK_FIXQ]
    389: . -mat_umfpack_aggressive <1>          - Control[UMFPACK_AGGRESSIVE]
    390: . -mat_umfpack_pivot_tolerance <delta> - UMFPACK partial pivot tolerance: Control[UMFPACK_PIVOT_TOLERANCE]
    391: . -mat_umfpack_sym_pivot_tolerance <0.001> - Control[UMFPACK_SYM_PIVOT_TOLERANCE]
    392: . -mat_umfpack_scale <NONE>           - (choose one of) NONE SUM MAX
    393: . -mat_umfpack_alloc_init <delta>      - UMFPACK factorized matrix allocation modifier: Control[UMFPACK_ALLOC_INIT]
    394: . -mat_umfpack_droptol <0>            - Control[UMFPACK_DROPTOL]
    395: - -mat_umfpack_irstep <maxit>          - UMFPACK maximum number of iterative refinement steps: Control[UMFPACK_IRSTEP]

    397:    Level: beginner

    399:    Note:
    400:    UMFPACK is part of SuiteSparse <http://faculty.cse.tamu.edu/davis/suitesparse.html>

    402: .seealso: [](ch_matrices), `Mat`, `PCLU`, `MATSOLVERSUPERLU`, `MATSOLVERMUMPS`, `PCFactorSetMatSolverType()`, `MatSolverType`
    403: M*/

    405: static PetscErrorCode MatGetFactor_seqaij_umfpack(Mat A, MatFactorType ftype, Mat *F)
    406: {
    407:   Mat          B;
    408:   Mat_UMFPACK *lu;
    409:   PetscInt     m = A->rmap->n, n = A->cmap->n;

    411:   PetscFunctionBegin;
    412:   /* Create the factorization matrix F */
    413:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
    414:   PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n));
    415:   PetscCall(PetscStrallocpy("umfpack", &((PetscObject)B)->type_name));
    416:   PetscCall(MatSetUp(B));

    418:   PetscCall(PetscNew(&lu));

    420:   B->data                  = lu;
    421:   B->ops->getinfo          = MatGetInfo_External;
    422:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_UMFPACK;
    423:   B->ops->destroy          = MatDestroy_UMFPACK;
    424:   B->ops->view             = MatView_UMFPACK;
    425:   B->ops->matsolve         = NULL;

    427:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_umfpack));

    429:   B->factortype   = MAT_FACTOR_LU;
    430:   B->assembled    = PETSC_TRUE; /* required by -ksp_view */
    431:   B->preallocated = PETSC_TRUE;

    433:   PetscCall(PetscFree(B->solvertype));
    434:   PetscCall(PetscStrallocpy(MATSOLVERUMFPACK, &B->solvertype));
    435:   B->canuseordering = PETSC_TRUE;
    436:   PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&B->preferredordering[MAT_FACTOR_LU]));

    438:   /* initializations */
    439:   /* get the default control parameters */
    440:   umfpack_UMF_defaults(lu->Control);
    441:   lu->perm_c                  = NULL; /* use default UMFPACK col permutation */
    442:   lu->Control[UMFPACK_IRSTEP] = 0;    /* max num of iterative refinement steps to attempt */

    444:   *F = B;
    445:   PetscFunctionReturn(PETSC_SUCCESS);
    446: }

    448: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_cholmod(Mat, MatFactorType, Mat *);
    449: PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat, MatFactorType, Mat *);
    450: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat, MatFactorType, Mat *);
    451: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_spqr(Mat, MatFactorType, Mat *);

    453: PETSC_INTERN PetscErrorCode MatSolverTypeRegister_SuiteSparse(void)
    454: {
    455:   PetscFunctionBegin;
    456:   PetscCall(MatSolverTypeRegister(MATSOLVERUMFPACK, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_umfpack));
    457:   PetscCall(MatSolverTypeRegister(MATSOLVERCHOLMOD, MATSEQAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_seqaij_cholmod));
    458:   PetscCall(MatSolverTypeRegister(MATSOLVERCHOLMOD, MATSEQSBAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_seqsbaij_cholmod));
    459:   PetscCall(MatSolverTypeRegister(MATSOLVERKLU, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_klu));
    460:   PetscCall(MatSolverTypeRegister(MATSOLVERSPQR, MATSEQAIJ, MAT_FACTOR_QR, MatGetFactor_seqaij_spqr));
    461:   if (!PetscDefined(USE_COMPLEX)) PetscCall(MatSolverTypeRegister(MATSOLVERSPQR, MATNORMAL, MAT_FACTOR_QR, MatGetFactor_seqaij_spqr));
    462:   PetscCall(MatSolverTypeRegister(MATSOLVERSPQR, MATNORMALHERMITIAN, MAT_FACTOR_QR, MatGetFactor_seqaij_spqr));
    463:   PetscFunctionReturn(PETSC_SUCCESS);
    464: }