Actual source code: cholesky.c

      1: /*
      2:    Defines a direct factorization preconditioner for any Mat implementation
      3:    Note: this need not be considered a preconditioner since it supplies
      4:          a direct solver.
      5: */
      6: #include <../src/ksp/pc/impls/factor/factor.h>

      8: typedef struct {
      9:   PC_Factor hdr;
     10:   IS        row, col; /* index sets used for reordering */
     11: } PC_Cholesky;

     13: static PetscErrorCode PCSetFromOptions_Cholesky(PC pc, PetscOptionItems PetscOptionsObject)
     14: {
     15:   PetscFunctionBegin;
     16:   PetscOptionsHeadBegin(PetscOptionsObject, "Cholesky options");
     17:   PetscCall(PCSetFromOptions_Factor(pc, PetscOptionsObject));
     18:   PetscOptionsHeadEnd();
     19:   PetscFunctionReturn(PETSC_SUCCESS);
     20: }

     22: static PetscErrorCode PCSetUp_Cholesky(PC pc)
     23: {
     24:   PetscBool      flg;
     25:   PC_Cholesky   *dir = (PC_Cholesky *)pc->data;
     26:   MatSolverType  stype;
     27:   MatFactorError err;
     28:   const char    *prefix;

     30:   PetscFunctionBegin;
     31:   pc->failedreason = PC_NOERROR;
     32:   if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor *)dir)->info.fill = dir->hdr.actualfill;

     34:   PetscCall(PCGetOptionsPrefix(pc, &prefix));
     35:   PetscCall(MatSetOptionsPrefixFactor(pc->pmat, prefix));

     37:   PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
     38:   if (dir->hdr.inplace) {
     39:     MatFactorType ftype;

     41:     PetscCall(MatGetFactorType(pc->pmat, &ftype));
     42:     if (ftype == MAT_FACTOR_NONE) {
     43:       if (dir->row && dir->col && (dir->row != dir->col)) PetscCall(ISDestroy(&dir->row));
     44:       PetscCall(ISDestroy(&dir->col));
     45:       /* This should only get the ordering if needed, but since MatGetFactor() is not called we can't know if it is needed */
     46:       PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
     47:       PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col));
     48:       if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */
     49:         PetscCall(ISDestroy(&dir->col));
     50:       }
     51:       PetscCall(MatCholeskyFactor(pc->pmat, dir->row, &((PC_Factor *)dir)->info));
     52:       PetscCall(MatFactorGetError(pc->pmat, &err));
     53:       if (err) { /* Factor() fails */
     54:         pc->failedreason = (PCFailedReason)err;
     55:         PetscFunctionReturn(PETSC_SUCCESS);
     56:       }
     57:     }
     58:     ((PC_Factor *)dir)->fact = pc->pmat;
     59:   } else {
     60:     MatInfo info;

     62:     if (!pc->setupcalled) {
     63:       PetscBool canuseordering;

     65:       PetscCall(PCFactorSetUpMatSolverType(pc));
     66:       PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering));
     67:       if (canuseordering) {
     68:         PetscBool external;

     70:         PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
     71:         PetscCall(PetscStrcmp(((PC_Factor *)dir)->ordering, MATORDERINGEXTERNAL, &external));
     72:         if (!external) {
     73:           PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col));
     74:           /* check if dir->row == dir->col */
     75:           if (dir->row) {
     76:             PetscCall(ISEqual(dir->row, dir->col, &flg));
     77:             PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "row and column permutations must be equal");
     78:           }
     79:           PetscCall(ISDestroy(&dir->col)); /* only pass one ordering into CholeskyFactor */

     81:           flg = PETSC_FALSE;
     82:           PetscCall(PetscOptionsHasName(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &flg));
     83:           if (flg) {
     84:             PetscReal tol = 1.e-10;
     85:             PetscCall(PetscOptionsGetReal(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &tol, NULL));
     86:             PetscCall(MatReorderForNonzeroDiagonal(pc->pmat, tol, dir->row, dir->row));
     87:           }
     88:         }
     89:       }
     90:       PetscCall(MatCholeskyFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->row, &((PC_Factor *)dir)->info));
     91:       PetscCall(MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info));
     92:       dir->hdr.actualfill = info.fill_ratio_needed;
     93:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
     94:       PetscCall(MatDestroy(&((PC_Factor *)dir)->fact));
     95:       PetscCall(PCFactorSetUpMatSolverType(pc));
     96:       if (!dir->hdr.reuseordering) {
     97:         PetscBool canuseordering;

     99:         PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering));
    100:         if (canuseordering) {
    101:           PetscBool external;

    103:           PetscCall(ISDestroy(&dir->row));
    104:           PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
    105:           PetscCall(PetscStrcmp(((PC_Factor *)dir)->ordering, MATORDERINGEXTERNAL, &external));
    106:           if (!external) {
    107:             PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col));
    108:             PetscCall(ISDestroy(&dir->col)); /* only use dir->row ordering in CholeskyFactor */

    110:             flg = PETSC_FALSE;
    111:             PetscCall(PetscOptionsHasName(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &flg));
    112:             if (flg) {
    113:               PetscReal tol = 1.e-10;
    114:               PetscCall(PetscOptionsGetReal(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &tol, NULL));
    115:               PetscCall(MatReorderForNonzeroDiagonal(pc->pmat, tol, dir->row, dir->row));
    116:             }
    117:           }
    118:         }
    119:       }
    120:       PetscCall(MatCholeskyFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->row, &((PC_Factor *)dir)->info));
    121:       PetscCall(MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info));
    122:       dir->hdr.actualfill = info.fill_ratio_needed;
    123:     } else {
    124:       PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
    125:       if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) {
    126:         PetscCall(MatFactorClearError(((PC_Factor *)dir)->fact));
    127:         pc->failedreason = PC_NOERROR;
    128:       }
    129:     }
    130:     PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
    131:     if (err) { /* FactorSymbolic() fails */
    132:       pc->failedreason = (PCFailedReason)err;
    133:       PetscFunctionReturn(PETSC_SUCCESS);
    134:     }

    136:     PetscCall(MatCholeskyFactorNumeric(((PC_Factor *)dir)->fact, pc->pmat, &((PC_Factor *)dir)->info));
    137:     PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
    138:     if (err) { /* FactorNumeric() fails */
    139:       pc->failedreason = (PCFailedReason)err;
    140:     }
    141:   }

    143:   PetscCall(PCFactorGetMatSolverType(pc, &stype));
    144:   if (!stype) {
    145:     MatSolverType solverpackage;
    146:     PetscCall(MatFactorGetSolverType(((PC_Factor *)dir)->fact, &solverpackage));
    147:     PetscCall(PCFactorSetMatSolverType(pc, solverpackage));
    148:   }
    149:   PetscFunctionReturn(PETSC_SUCCESS);
    150: }

    152: static PetscErrorCode PCReset_Cholesky(PC pc)
    153: {
    154:   PC_Cholesky *dir = (PC_Cholesky *)pc->data;

    156:   PetscFunctionBegin;
    157:   if (!dir->hdr.inplace && ((PC_Factor *)dir)->fact) PetscCall(MatDestroy(&((PC_Factor *)dir)->fact));
    158:   PetscCall(ISDestroy(&dir->row));
    159:   PetscCall(ISDestroy(&dir->col));
    160:   PetscFunctionReturn(PETSC_SUCCESS);
    161: }

    163: static PetscErrorCode PCDestroy_Cholesky(PC pc)
    164: {
    165:   PC_Cholesky *dir = (PC_Cholesky *)pc->data;

    167:   PetscFunctionBegin;
    168:   PetscCall(PCReset_Cholesky(pc));
    169:   PetscCall(PetscFree(((PC_Factor *)dir)->ordering));
    170:   PetscCall(PetscFree(((PC_Factor *)dir)->solvertype));
    171:   PetscCall(PCFactorClearComposedFunctions(pc));
    172:   PetscCall(PetscFree(pc->data));
    173:   PetscFunctionReturn(PETSC_SUCCESS);
    174: }

    176: static PetscErrorCode PCApply_Cholesky(PC pc, Vec x, Vec y)
    177: {
    178:   PC_Cholesky *dir = (PC_Cholesky *)pc->data;

    180:   PetscFunctionBegin;
    181:   if (dir->hdr.inplace) {
    182:     PetscCall(MatSolve(pc->pmat, x, y));
    183:   } else {
    184:     PetscCall(MatSolve(((PC_Factor *)dir)->fact, x, y));
    185:   }
    186:   PetscFunctionReturn(PETSC_SUCCESS);
    187: }

    189: static PetscErrorCode PCMatApply_Cholesky(PC pc, Mat X, Mat Y)
    190: {
    191:   PC_Cholesky *dir = (PC_Cholesky *)pc->data;

    193:   PetscFunctionBegin;
    194:   if (dir->hdr.inplace) {
    195:     PetscCall(MatMatSolve(pc->pmat, X, Y));
    196:   } else {
    197:     PetscCall(MatMatSolve(((PC_Factor *)dir)->fact, X, Y));
    198:   }
    199:   PetscFunctionReturn(PETSC_SUCCESS);
    200: }

    202: static PetscErrorCode PCApplySymmetricLeft_Cholesky(PC pc, Vec x, Vec y)
    203: {
    204:   PC_Cholesky *dir = (PC_Cholesky *)pc->data;

    206:   PetscFunctionBegin;
    207:   if (dir->hdr.inplace) {
    208:     PetscCall(MatForwardSolve(pc->pmat, x, y));
    209:   } else {
    210:     PetscCall(MatForwardSolve(((PC_Factor *)dir)->fact, x, y));
    211:   }
    212:   PetscFunctionReturn(PETSC_SUCCESS);
    213: }

    215: static PetscErrorCode PCApplySymmetricRight_Cholesky(PC pc, Vec x, Vec y)
    216: {
    217:   PC_Cholesky *dir = (PC_Cholesky *)pc->data;

    219:   PetscFunctionBegin;
    220:   if (dir->hdr.inplace) {
    221:     PetscCall(MatBackwardSolve(pc->pmat, x, y));
    222:   } else {
    223:     PetscCall(MatBackwardSolve(((PC_Factor *)dir)->fact, x, y));
    224:   }
    225:   PetscFunctionReturn(PETSC_SUCCESS);
    226: }

    228: static PetscErrorCode PCApplyTranspose_Cholesky(PC pc, Vec x, Vec y)
    229: {
    230:   PC_Cholesky *dir = (PC_Cholesky *)pc->data;

    232:   PetscFunctionBegin;
    233:   if (dir->hdr.inplace) {
    234:     PetscCall(MatSolveTranspose(pc->pmat, x, y));
    235:   } else {
    236:     PetscCall(MatSolveTranspose(((PC_Factor *)dir)->fact, x, y));
    237:   }
    238:   PetscFunctionReturn(PETSC_SUCCESS);
    239: }

    241: static PetscErrorCode PCMatApplyTranspose_Cholesky(PC pc, Mat X, Mat Y)
    242: {
    243:   PC_Cholesky *dir = (PC_Cholesky *)pc->data;

    245:   PetscFunctionBegin;
    246:   if (dir->hdr.inplace) {
    247:     PetscCall(MatMatSolveTranspose(pc->pmat, X, Y));
    248:   } else {
    249:     PetscCall(MatMatSolveTranspose(((PC_Factor *)dir)->fact, X, Y));
    250:   }
    251:   PetscFunctionReturn(PETSC_SUCCESS);
    252: }

    254: /*@
    255:   PCFactorSetReuseOrdering - When similar matrices are factored, this
    256:   causes the ordering computed in the first factor to be used for all
    257:   following factors.

    259:   Logically Collective

    261:   Input Parameters:
    262: + pc   - the preconditioner context
    263: - flag - `PETSC_TRUE` to reuse else `PETSC_FALSE`

    265:   Options Database Key:
    266: . -pc_factor_reuse_ordering - Activate `PCFactorSetReuseOrdering()`

    268:   Level: intermediate

    270: .seealso: [](ch_ksp), `PCLU`, `PCCHOLESKY`, `PCFactorSetReuseFill()`
    271: @*/
    272: PetscErrorCode PCFactorSetReuseOrdering(PC pc, PetscBool flag)
    273: {
    274:   PetscFunctionBegin;
    277:   PetscTryMethod(pc, "PCFactorSetReuseOrdering_C", (PC, PetscBool), (pc, flag));
    278:   PetscFunctionReturn(PETSC_SUCCESS);
    279: }

    281: /*MC
    282:    PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner

    284:    Options Database Keys:
    285: +  -pc_factor_reuse_ordering (true|false) - Activate `PCFactorSetReuseOrdering()`
    286: .  -pc_factor_mat_solver_type type        - Actives `PCFactorSetMatSolverType()` to choose the direct solver, like `superlu`
    287: .  -pc_factor_reuse_fill (true|false)     - Activates `PCFactorSetReuseFill()`
    288: .  -pc_factor_fill fill                   - Sets the explected fill amount
    289: .  -pc_factor_in_place (true|false)       - Activates in-place factorization
    290: .  -pc_factor_mat_ordering_type type      - Sets ordering routine used to determine the order the rows are used in the factorization to reduce fill, see `MatOrderingType`
    291: -  -mat_solvertype_optionname             - options for a specific solver package, for example `-mat_mumps_cntl_1`

    293:    Level: beginner

    295:    Notes:
    296:    The Cholesky factorization direct solver, `PCCHOLESKY` is only for symmetric positive-definite (SPD) matrices. For such
    297:    SPD matrices it is more efficient than using the LU factorization direct solver, `PCLU`.

    299:    Not all options work for all matrix formats

    301:    Usually this will compute an "exact" solution in one iteration and does
    302:    not need a Krylov method (i.e. you can use -ksp_type preonly, or
    303:    `KSPSetType`(ksp,`KSPPREONLY`) for the Krylov method

    305:    The options prefix of the factored matrix is set to be the same as the `PC` options prefix.

    307: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
    308:           `PCILU`, `PCLU`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`,
    309:           `PCFactorSetFill()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`,
    310:           `PCFactorSetUseInPlace()`, `PCFactorGetUseInPlace()`, `PCFactorSetMatOrderingType()`
    311: M*/

    313: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
    314: {
    315:   PC_Cholesky *dir;

    317:   PetscFunctionBegin;
    318:   PetscCall(PetscNew(&dir));
    319:   pc->data = (void *)dir;
    320:   PetscCall(PCFactorInitialize(pc, MAT_FACTOR_CHOLESKY));

    322:   ((PC_Factor *)dir)->info.fill = 5.0;

    324:   pc->ops->destroy             = PCDestroy_Cholesky;
    325:   pc->ops->reset               = PCReset_Cholesky;
    326:   pc->ops->apply               = PCApply_Cholesky;
    327:   pc->ops->matapply            = PCMatApply_Cholesky;
    328:   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_Cholesky;
    329:   pc->ops->applysymmetricright = PCApplySymmetricRight_Cholesky;
    330:   pc->ops->applytranspose      = PCApplyTranspose_Cholesky;
    331:   pc->ops->matapplytranspose   = PCMatApplyTranspose_Cholesky;
    332:   pc->ops->setup               = PCSetUp_Cholesky;
    333:   pc->ops->setfromoptions      = PCSetFromOptions_Cholesky;
    334:   pc->ops->view                = PCView_Factor;
    335:   pc->ops->applyrichardson     = NULL;
    336:   PetscFunctionReturn(PETSC_SUCCESS);
    337: }