Actual source code: klu.c

      1: /*
      2:    Provides an interface to the KLUv1.2 sparse solver

      4:    When build with PETSC_USE_64BIT_INDICES this will use SuiteSparse_long as the
      5:    integer type in KLU, otherwise it will use int. This means
      6:    all integers in this file are simply declared as PetscInt. Also it means
      7:    that KLU SuiteSparse_long version MUST be built with 64-bit integers when used.

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

     12: #if defined(PETSC_USE_64BIT_INDICES)
     13:   #define klu_K_defaults                        klu_l_defaults
     14:   #define klu_K_analyze(a, b, c, d)             klu_l_analyze((SuiteSparse_long)a, (SuiteSparse_long *)b, (SuiteSparse_long *)c, d)
     15:   #define klu_K_analyze_given(a, b, c, d, e, f) klu_l_analyze_given((SuiteSparse_long)a, (SuiteSparse_long *)b, (SuiteSparse_long *)c, (SuiteSparse_long *)d, (SuiteSparse_long *)e, f)
     16:   #define klu_K_free_symbolic                   klu_l_free_symbolic
     17:   #define klu_K_free_numeric                    klu_l_free_numeric
     18:   #define klu_K_common                          klu_l_common
     19:   #define klu_K_symbolic                        klu_l_symbolic
     20:   #define klu_K_numeric                         klu_l_numeric
     21:   #if defined(PETSC_USE_COMPLEX)
     22:     #define klu_K_factor(a, b, c, d, e)       klu_zl_factor((SuiteSparse_long *)a, (SuiteSparse_long *)b, c, d, e);
     23:     #define klu_K_solve                       klu_zl_solve
     24:     #define klu_K_tsolve(a, b, c, d, e, f, g) klu_zl_tsolve((a), (b), (c), (d), (PetscReal *)(e), (f), (g))
     25:     #define klu_K_refactor                    klu_zl_refactor
     26:     #define klu_K_sort                        klu_zl_sort
     27:     #define klu_K_flops                       klu_zl_flops
     28:     #define klu_K_rgrowth                     klu_zl_rgrowth
     29:     #define klu_K_condest                     klu_zl_condest
     30:     #define klu_K_rcond                       klu_zl_rcond
     31:     #define klu_K_scale                       klu_zl_scale
     32:   #else
     33:     #define klu_K_factor(a, b, c, d, e)       klu_l_factor((SuiteSparse_long *)a, (SuiteSparse_long *)b, c, d, e);
     34:     #define klu_K_solve                       klu_l_solve
     35:     #define klu_K_tsolve(a, b, c, d, e, f, g) klu_l_tsolve((a), (b), (c), (d), (e), (g))
     36:     #define klu_K_refactor                    klu_l_refactor
     37:     #define klu_K_sort                        klu_l_sort
     38:     #define klu_K_flops                       klu_l_flops
     39:     #define klu_K_rgrowth                     klu_l_rgrowth
     40:     #define klu_K_condest                     klu_l_condest
     41:     #define klu_K_rcond                       klu_l_rcond
     42:     #define klu_K_scale                       klu_l_scale
     43:   #endif
     44: #else
     45:   #define klu_K_defaults      klu_defaults
     46:   #define klu_K_analyze       klu_analyze
     47:   #define klu_K_analyze_given klu_analyze_given
     48:   #define klu_K_free_symbolic klu_free_symbolic
     49:   #define klu_K_free_numeric  klu_free_numeric
     50:   #define klu_K_common        klu_common
     51:   #define klu_K_symbolic      klu_symbolic
     52:   #define klu_K_numeric       klu_numeric
     53:   #if defined(PETSC_USE_COMPLEX)
     54:     #define klu_K_factor                      klu_z_factor
     55:     #define klu_K_solve                       klu_z_solve
     56:     #define klu_K_tsolve(a, b, c, d, e, f, g) klu_z_tsolve((a), (b), (c), (d), (PetscReal *)(e), (f), (g))
     57:     #define klu_K_refactor                    klu_z_refactor
     58:     #define klu_K_sort                        klu_z_sort
     59:     #define klu_K_flops                       klu_z_flops
     60:     #define klu_K_rgrowth                     klu_z_rgrowth
     61:     #define klu_K_condest                     klu_z_condest
     62:     #define klu_K_rcond                       klu_z_rcond
     63:     #define klu_K_scale                       klu_z_scale
     64:   #else
     65:     #define klu_K_factor                      klu_factor
     66:     #define klu_K_solve                       klu_solve
     67:     #define klu_K_tsolve(a, b, c, d, e, f, g) klu_tsolve((a), (b), (c), (d), (e), (g))
     68:     #define klu_K_refactor                    klu_refactor
     69:     #define klu_K_sort                        klu_sort
     70:     #define klu_K_flops                       klu_flops
     71:     #define klu_K_rgrowth                     klu_rgrowth
     72:     #define klu_K_condest                     klu_condest
     73:     #define klu_K_rcond                       klu_rcond
     74:     #define klu_K_scale                       klu_scale
     75:   #endif
     76: #endif

     78: EXTERN_C_BEGIN
     79: #include <klu.h>
     80: EXTERN_C_END

     82: static const char *KluOrderingTypes[] = {"AMD", "COLAMD"};
     83: static const char *scale[]            = {"NONE", "SUM", "MAX"};

     85: typedef struct {
     86:   klu_K_common    Common;
     87:   klu_K_symbolic *Symbolic;
     88:   klu_K_numeric  *Numeric;
     89:   PetscInt       *perm_c, *perm_r;
     90:   MatStructure    flg;
     91:   PetscBool       PetscMatOrdering;
     92:   PetscBool       CleanUpKLU;
     93: } Mat_KLU;

     95: static PetscErrorCode MatDestroy_KLU(Mat A)
     96: {
     97:   Mat_KLU *lu = (Mat_KLU *)A->data;

     99:   PetscFunctionBegin;
    100:   if (lu->CleanUpKLU) {
    101:     klu_K_free_symbolic(&lu->Symbolic, &lu->Common);
    102:     klu_K_free_numeric(&lu->Numeric, &lu->Common);
    103:     PetscCall(PetscFree2(lu->perm_r, lu->perm_c));
    104:   }
    105:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL));
    106:   PetscCall(PetscFree(A->data));
    107:   PetscFunctionReturn(PETSC_SUCCESS);
    108: }

    110: static PetscErrorCode MatSolveTranspose_KLU(Mat A, Vec b, Vec x)
    111: {
    112:   Mat_KLU     *lu = (Mat_KLU *)A->data;
    113:   PetscScalar *xa;
    114:   PetscInt     status;

    116:   PetscFunctionBegin;
    117:   /* KLU uses a column major format, solve Ax = b by klu_*_solve */
    118:   PetscCall(VecCopy(b, x)); /* klu_solve stores the solution in rhs */
    119:   PetscCall(VecGetArray(x, &xa));
    120:   status = klu_K_solve(lu->Symbolic, lu->Numeric, A->rmap->n, 1, (PetscReal *)xa, &lu->Common);
    121:   PetscCheck(status == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "KLU Solve failed");
    122:   PetscCall(VecRestoreArray(x, &xa));
    123:   PetscFunctionReturn(PETSC_SUCCESS);
    124: }

    126: static PetscErrorCode MatSolve_KLU(Mat A, Vec b, Vec x)
    127: {
    128:   Mat_KLU     *lu = (Mat_KLU *)A->data;
    129:   PetscScalar *xa;

    131:   PetscFunctionBegin;
    132:   /* KLU uses a column major format, solve A^Tx = b by klu_*_tsolve */
    133:   PetscCall(VecCopy(b, x)); /* klu_solve stores the solution in rhs */
    134:   PetscCall(VecGetArray(x, &xa));
    135:   PetscCheck(klu_K_tsolve(lu->Symbolic, lu->Numeric, A->rmap->n, 1, xa, 1, &lu->Common), PETSC_COMM_SELF, PETSC_ERR_LIB, "KLU Solve failed"); /* conjugate solve */
    136:   PetscCall(VecRestoreArray(x, &xa));
    137:   PetscFunctionReturn(PETSC_SUCCESS);
    138: }

    140: static PetscErrorCode MatLUFactorNumeric_KLU(Mat F, Mat A, const MatFactorInfo *info)
    141: {
    142:   Mat_KLU     *lu = (Mat_KLU *)F->data;
    143:   Mat_SeqAIJ  *a  = (Mat_SeqAIJ *)A->data;
    144:   PetscInt    *ai = a->i, *aj = a->j;
    145:   PetscScalar *av = a->a;

    147:   PetscFunctionBegin;
    148:   /* numeric factorization of A' */

    150:   if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric) klu_K_free_numeric(&lu->Numeric, &lu->Common);
    151:   lu->Numeric = klu_K_factor(ai, aj, (PetscReal *)av, lu->Symbolic, &lu->Common);
    152:   PetscCheck(lu->Numeric, PETSC_COMM_SELF, PETSC_ERR_LIB, "KLU Numeric factorization failed");

    154:   lu->flg                = SAME_NONZERO_PATTERN;
    155:   lu->CleanUpKLU         = PETSC_TRUE;
    156:   F->ops->solve          = MatSolve_KLU;
    157:   F->ops->solvetranspose = MatSolveTranspose_KLU;
    158:   PetscFunctionReturn(PETSC_SUCCESS);
    159: }

    161: static PetscErrorCode MatLUFactorSymbolic_KLU(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
    162: {
    163:   Mat_SeqAIJ     *a  = (Mat_SeqAIJ *)A->data;
    164:   Mat_KLU        *lu = (Mat_KLU *)F->data;
    165:   PetscInt        i, *ai = a->i, *aj = a->j, m = A->rmap->n, n = A->cmap->n;
    166:   const PetscInt *ra, *ca;

    168:   PetscFunctionBegin;
    169:   if (lu->PetscMatOrdering) {
    170:     PetscCall(ISGetIndices(r, &ra));
    171:     PetscCall(ISGetIndices(c, &ca));
    172:     PetscCall(PetscMalloc2(m, &lu->perm_r, n, &lu->perm_c));
    173:     /* we cannot simply memcpy on 64-bit archs */
    174:     for (i = 0; i < m; i++) lu->perm_r[i] = ra[i];
    175:     for (i = 0; i < n; i++) lu->perm_c[i] = ca[i];
    176:     PetscCall(ISRestoreIndices(r, &ra));
    177:     PetscCall(ISRestoreIndices(c, &ca));
    178:   }

    180:   /* symbolic factorization of A' */
    181:   if (r) {
    182:     lu->PetscMatOrdering = PETSC_TRUE;
    183:     lu->Symbolic         = klu_K_analyze_given(n, ai, aj, lu->perm_c, lu->perm_r, &lu->Common);
    184:   } else { /* use klu internal ordering */
    185:     lu->Symbolic = klu_K_analyze(n, ai, aj, &lu->Common);
    186:   }
    187:   PetscCheck(lu->Symbolic, PETSC_COMM_SELF, PETSC_ERR_LIB, "KLU Symbolic Factorization failed");

    189:   lu->flg                 = DIFFERENT_NONZERO_PATTERN;
    190:   lu->CleanUpKLU          = PETSC_TRUE;
    191:   F->ops->lufactornumeric = MatLUFactorNumeric_KLU;
    192:   PetscFunctionReturn(PETSC_SUCCESS);
    193: }

    195: static PetscErrorCode MatView_Info_KLU(Mat A, PetscViewer viewer)
    196: {
    197:   Mat_KLU       *lu      = (Mat_KLU *)A->data;
    198:   klu_K_numeric *Numeric = (klu_K_numeric *)lu->Numeric;

    200:   PetscFunctionBegin;
    201:   PetscCall(PetscViewerASCIIPrintf(viewer, "KLU stats:\n"));
    202:   PetscCall(PetscViewerASCIIPrintf(viewer, "  Number of diagonal blocks: %" PetscInt_FMT "\n", (PetscInt)Numeric->nblocks));
    203:   PetscCall(PetscViewerASCIIPrintf(viewer, "  Total nonzeros=%" PetscInt_FMT "\n", (PetscInt)(Numeric->lnz + Numeric->unz)));
    204:   PetscCall(PetscViewerASCIIPrintf(viewer, "KLU runtime parameters:\n"));
    205:   /* Control parameters used by numeric factorization */
    206:   PetscCall(PetscViewerASCIIPrintf(viewer, "  Partial pivoting tolerance: %g\n", lu->Common.tol));
    207:   /* BTF preordering */
    208:   PetscCall(PetscViewerASCIIPrintf(viewer, "  BTF preordering enabled: %" PetscInt_FMT "\n", (PetscInt)lu->Common.btf));
    209:   /* mat ordering */
    210:   if (!lu->PetscMatOrdering) PetscCall(PetscViewerASCIIPrintf(viewer, "  Ordering: %s (not using the PETSc ordering)\n", KluOrderingTypes[(int)lu->Common.ordering]));
    211:   /* matrix row scaling */
    212:   PetscCall(PetscViewerASCIIPrintf(viewer, "  Matrix row scaling: %s\n", scale[(int)lu->Common.scale]));
    213:   PetscFunctionReturn(PETSC_SUCCESS);
    214: }

    216: static PetscErrorCode MatView_KLU(Mat A, PetscViewer viewer)
    217: {
    218:   PetscBool         isascii;
    219:   PetscViewerFormat format;

    221:   PetscFunctionBegin;
    222:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    223:   if (isascii) {
    224:     PetscCall(PetscViewerGetFormat(viewer, &format));
    225:     if (format == PETSC_VIEWER_ASCII_INFO) PetscCall(MatView_Info_KLU(A, viewer));
    226:   }
    227:   PetscFunctionReturn(PETSC_SUCCESS);
    228: }

    230: static PetscErrorCode MatFactorGetSolverType_seqaij_klu(Mat A, MatSolverType *type)
    231: {
    232:   PetscFunctionBegin;
    233:   *type = MATSOLVERKLU;
    234:   PetscFunctionReturn(PETSC_SUCCESS);
    235: }

    237: /*MC
    238:   MATSOLVERKLU = "klu" - A matrix type providing direct solvers, LU, for sequential matrices
    239:   via the external package KLU.

    241:   `./configure --download-suitesparse` to install PETSc to use KLU

    243:   Use `-pc_type lu` `-pc_factor_mat_solver_type klu` to use this direct solver

    245:   Consult KLU documentation for more information on the options database keys below.

    247:   Options Database Keys:
    248: + -mat_klu_pivot_tol 0.001 - Partial pivoting tolerance
    249: . -mat_klu_use_btf 1       - Use BTF preordering
    250: . -mat_klu_ordering AMD    - KLU reordering scheme to reduce fill-in (choose one of) `AMD`, `COLAMD`, `PETSC`
    251: - -mat_klu_row_scale NONE  - Matrix row scaling (choose one of) `NONE`, `SUM`, `MAX`

    253:    Level: beginner

    255:    Note:
    256:    KLU is part of SuiteSparse <http://faculty.cse.tamu.edu/davis/suitesparse.html>

    258: .seealso: [](ch_matrices), `Mat`, `PCLU`, `MATSOLVERUMFPACK`, `MATSOLVERCHOLMOD`, `PCFactorSetMatSolverType()`, `MatSolverType`
    259: M*/

    261: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat A, MatFactorType ftype, Mat *F)
    262: {
    263:   Mat       B;
    264:   Mat_KLU  *lu;
    265:   PetscInt  m = A->rmap->n, n = A->cmap->n, idx = 0, status;
    266:   PetscBool flg;

    268:   PetscFunctionBegin;
    269:   /* Create the factorization matrix F */
    270:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
    271:   PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n));
    272:   PetscCall(PetscStrallocpy("klu", &((PetscObject)B)->type_name));
    273:   PetscCall(MatSetUp(B));

    275:   PetscCall(PetscNew(&lu));

    277:   B->data                  = lu;
    278:   B->ops->getinfo          = MatGetInfo_External;
    279:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_KLU;
    280:   B->ops->destroy          = MatDestroy_KLU;
    281:   B->ops->view             = MatView_KLU;

    283:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_klu));

    285:   B->factortype   = MAT_FACTOR_LU;
    286:   B->assembled    = PETSC_TRUE; /* required by -ksp_view */
    287:   B->preallocated = PETSC_TRUE;

    289:   PetscCall(PetscFree(B->solvertype));
    290:   PetscCall(PetscStrallocpy(MATSOLVERKLU, &B->solvertype));
    291:   B->canuseordering = PETSC_TRUE;
    292:   PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&B->preferredordering[MAT_FACTOR_LU]));

    294:   /* initializations */
    295:   /* get the default control parameters */
    296:   status = klu_K_defaults(&lu->Common);
    297:   PetscCheck(status > 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "KLU Initialization failed");

    299:   lu->Common.scale = 0; /* No row scaling */

    301:   PetscOptionsBegin(PetscObjectComm((PetscObject)B), ((PetscObject)B)->prefix, "KLU Options", "Mat");
    302:   /* Partial pivoting tolerance */
    303:   PetscCall(PetscOptionsReal("-mat_klu_pivot_tol", "Partial pivoting tolerance", "None", lu->Common.tol, &lu->Common.tol, NULL));
    304:   /* BTF pre-ordering */
    305:   PetscCall(PetscOptionsInt("-mat_klu_use_btf", "Enable BTF preordering", "None", (PetscInt)lu->Common.btf, (PetscInt *)&lu->Common.btf, NULL));
    306:   /* Matrix reordering */
    307:   PetscCall(PetscOptionsEList("-mat_klu_ordering", "Internal ordering method", "None", KluOrderingTypes, PETSC_STATIC_ARRAY_LENGTH(KluOrderingTypes), KluOrderingTypes[0], &idx, &flg));
    308:   lu->Common.ordering = (int)idx;
    309:   /* Matrix row scaling */
    310:   PetscCall(PetscOptionsEList("-mat_klu_row_scale", "Matrix row scaling", "None", scale, 3, scale[0], &idx, &flg));
    311:   PetscOptionsEnd();
    312:   *F = B;
    313:   PetscFunctionReturn(PETSC_SUCCESS);
    314: }