Actual source code: lu.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: */

      7: #include <../src/ksp/pc/impls/factor/lu/lu.h>

      9: static PetscErrorCode PCFactorReorderForNonzeroDiagonal_LU(PC pc, PetscReal z)
     10: {
     11:   PC_LU *lu = (PC_LU *)pc->data;

     13:   PetscFunctionBegin;
     14:   lu->nonzerosalongdiagonal = PETSC_TRUE;
     15:   if (z == (PetscReal)PETSC_DECIDE) lu->nonzerosalongdiagonaltol = 1.e-10;
     16:   else lu->nonzerosalongdiagonaltol = z;
     17:   PetscFunctionReturn(PETSC_SUCCESS);
     18: }

     20: static PetscErrorCode PCSetFromOptions_LU(PC pc, PetscOptionItems PetscOptionsObject)
     21: {
     22:   PC_LU    *lu  = (PC_LU *)pc->data;
     23:   PetscBool flg = PETSC_FALSE;
     24:   PetscReal tol;

     26:   PetscFunctionBegin;
     27:   PetscOptionsHeadBegin(PetscOptionsObject, "LU options");
     28:   PetscCall(PCSetFromOptions_Factor(pc, PetscOptionsObject));

     30:   PetscCall(PetscOptionsName("-pc_factor_nonzeros_along_diagonal", "Reorder to remove zeros from diagonal", "PCFactorReorderForNonzeroDiagonal", &flg));
     31:   if (flg) {
     32:     tol = PETSC_DECIDE;
     33:     PetscCall(PetscOptionsReal("-pc_factor_nonzeros_along_diagonal", "Reorder to remove zeros from diagonal", "PCFactorReorderForNonzeroDiagonal", lu->nonzerosalongdiagonaltol, &tol, NULL));
     34:     PetscCall(PCFactorReorderForNonzeroDiagonal(pc, tol));
     35:   }
     36:   PetscOptionsHeadEnd();
     37:   PetscFunctionReturn(PETSC_SUCCESS);
     38: }

     40: static PetscErrorCode PCSetUp_LU(PC pc)
     41: {
     42:   PC_LU         *dir = (PC_LU *)pc->data;
     43:   MatSolverType  stype;
     44:   MatFactorError err;
     45:   const char    *prefix;

     47:   PetscFunctionBegin;
     48:   pc->failedreason = PC_NOERROR;
     49:   if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor *)dir)->info.fill = dir->hdr.actualfill;

     51:   PetscCall(PCGetOptionsPrefix(pc, &prefix));
     52:   PetscCall(MatSetOptionsPrefixFactor(pc->pmat, prefix));

     54:   PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
     55:   if (dir->hdr.inplace) {
     56:     MatFactorType ftype;

     58:     PetscCall(MatGetFactorType(pc->pmat, &ftype));
     59:     if (ftype == MAT_FACTOR_NONE) {
     60:       if (dir->row && dir->col && dir->row != dir->col) PetscCall(ISDestroy(&dir->row));
     61:       PetscCall(ISDestroy(&dir->col));
     62:       /* This should only get the ordering if needed, but since MatGetFactor() is not called we can't know if it is needed */
     63:       PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
     64:       PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col));
     65:       PetscCall(MatLUFactor(pc->pmat, dir->row, dir->col, &((PC_Factor *)dir)->info));
     66:       PetscCall(MatFactorGetError(pc->pmat, &err));
     67:       if (err) { /* Factor() fails */
     68:         pc->failedreason = (PCFailedReason)err;
     69:         PetscFunctionReturn(PETSC_SUCCESS);
     70:       }
     71:     }
     72:     ((PC_Factor *)dir)->fact = pc->pmat;
     73:   } else {
     74:     MatInfo info;

     76:     if (!pc->setupcalled) {
     77:       PetscBool canuseordering;

     79:       PetscCall(PCFactorSetUpMatSolverType(pc));
     80:       PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering));
     81:       if (canuseordering) {
     82:         PetscBool external;

     84:         PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
     85:         PetscCall(PetscStrcmp(((PC_Factor *)dir)->ordering, MATORDERINGEXTERNAL, &external));
     86:         if (!external) {
     87:           PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col));
     88:           if (dir->nonzerosalongdiagonal) PetscCall(MatReorderForNonzeroDiagonal(pc->pmat, dir->nonzerosalongdiagonaltol, dir->row, dir->col));
     89:         }
     90:       }
     91:       PetscCall(MatLUFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->row, dir->col, &((PC_Factor *)dir)->info));
     92:       PetscCall(MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info));
     93:       dir->hdr.actualfill = info.fill_ratio_needed;
     94:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
     95:       PetscBool canuseordering;

     97:       PetscCall(MatDestroy(&((PC_Factor *)dir)->fact));
     98:       PetscCall(PCFactorSetUpMatSolverType(pc));
     99:       if (!dir->hdr.reuseordering) {
    100:         PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering));
    101:         if (canuseordering) {
    102:           PetscBool external;

    104:           if (dir->row && dir->col && dir->row != dir->col) PetscCall(ISDestroy(&dir->row));
    105:           PetscCall(ISDestroy(&dir->col));
    106:           PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
    107:           PetscCall(PetscStrcmp(((PC_Factor *)dir)->ordering, MATORDERINGEXTERNAL, &external));
    108:           if (!external) {
    109:             PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col));
    110:             if (dir->nonzerosalongdiagonal) PetscCall(MatReorderForNonzeroDiagonal(pc->pmat, dir->nonzerosalongdiagonaltol, dir->row, dir->col));
    111:           }
    112:         }
    113:       }
    114:       PetscCall(MatLUFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->row, dir->col, &((PC_Factor *)dir)->info));
    115:       PetscCall(MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info));
    116:       dir->hdr.actualfill = info.fill_ratio_needed;
    117:     } else {
    118:       PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
    119:       if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) {
    120:         PetscCall(MatFactorClearError(((PC_Factor *)dir)->fact));
    121:         pc->failedreason = PC_NOERROR;
    122:       }
    123:     }
    124:     PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
    125:     if (err) { /* FactorSymbolic() fails */
    126:       pc->failedreason = (PCFailedReason)err;
    127:       PetscFunctionReturn(PETSC_SUCCESS);
    128:     }

    130:     PetscCall(MatLUFactorNumeric(((PC_Factor *)dir)->fact, pc->pmat, &((PC_Factor *)dir)->info));
    131:     PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
    132:     if (err) { /* FactorNumeric() fails */
    133:       pc->failedreason = (PCFailedReason)err;
    134:     }
    135:   }

    137:   PetscCall(PCFactorGetMatSolverType(pc, &stype));
    138:   if (!stype) {
    139:     MatSolverType solverpackage;
    140:     PetscCall(MatFactorGetSolverType(((PC_Factor *)dir)->fact, &solverpackage));
    141:     PetscCall(PCFactorSetMatSolverType(pc, solverpackage));
    142:   }
    143:   PetscFunctionReturn(PETSC_SUCCESS);
    144: }

    146: static PetscErrorCode PCReset_LU(PC pc)
    147: {
    148:   PC_LU *dir = (PC_LU *)pc->data;

    150:   PetscFunctionBegin;
    151:   if (!dir->hdr.inplace && ((PC_Factor *)dir)->fact) PetscCall(MatDestroy(&((PC_Factor *)dir)->fact));
    152:   if (dir->row && dir->col && dir->row != dir->col) PetscCall(ISDestroy(&dir->row));
    153:   PetscCall(ISDestroy(&dir->col));
    154:   PetscFunctionReturn(PETSC_SUCCESS);
    155: }

    157: static PetscErrorCode PCDestroy_LU(PC pc)
    158: {
    159:   PC_LU *dir = (PC_LU *)pc->data;

    161:   PetscFunctionBegin;
    162:   PetscCall(PCReset_LU(pc));
    163:   PetscCall(PetscFree(((PC_Factor *)dir)->ordering));
    164:   PetscCall(PetscFree(((PC_Factor *)dir)->solvertype));
    165:   PetscCall(PCFactorClearComposedFunctions(pc));
    166:   PetscCall(PetscFree(pc->data));
    167:   PetscFunctionReturn(PETSC_SUCCESS);
    168: }

    170: static PetscErrorCode PCApply_LU(PC pc, Vec x, Vec y)
    171: {
    172:   PC_LU *dir = (PC_LU *)pc->data;

    174:   PetscFunctionBegin;
    175:   if (dir->hdr.inplace) {
    176:     PetscCall(MatSolve(pc->pmat, x, y));
    177:   } else {
    178:     PetscCall(MatSolve(((PC_Factor *)dir)->fact, x, y));
    179:   }
    180:   PetscFunctionReturn(PETSC_SUCCESS);
    181: }

    183: static PetscErrorCode PCMatApply_LU(PC pc, Mat X, Mat Y)
    184: {
    185:   PC_LU *dir = (PC_LU *)pc->data;

    187:   PetscFunctionBegin;
    188:   if (dir->hdr.inplace) {
    189:     PetscCall(MatMatSolve(pc->pmat, X, Y));
    190:   } else {
    191:     PetscCall(MatMatSolve(((PC_Factor *)dir)->fact, X, Y));
    192:   }
    193:   PetscFunctionReturn(PETSC_SUCCESS);
    194: }

    196: static PetscErrorCode PCApplyTranspose_LU(PC pc, Vec x, Vec y)
    197: {
    198:   PC_LU *dir = (PC_LU *)pc->data;

    200:   PetscFunctionBegin;
    201:   if (dir->hdr.inplace) {
    202:     PetscCall(MatSolveTranspose(pc->pmat, x, y));
    203:   } else {
    204:     PetscCall(MatSolveTranspose(((PC_Factor *)dir)->fact, x, y));
    205:   }
    206:   PetscFunctionReturn(PETSC_SUCCESS);
    207: }

    209: static PetscErrorCode PCMatApplyTranspose_LU(PC pc, Mat X, Mat Y)
    210: {
    211:   PC_LU *dir = (PC_LU *)pc->data;

    213:   PetscFunctionBegin;
    214:   if (dir->hdr.inplace) {
    215:     PetscCall(MatMatSolveTranspose(pc->pmat, X, Y));
    216:   } else {
    217:     PetscCall(MatMatSolveTranspose(((PC_Factor *)dir)->fact, X, Y));
    218:   }
    219:   PetscFunctionReturn(PETSC_SUCCESS);
    220: }

    222: /*MC
    223:    PCLU - Uses a direct solver, based on LU factorization, as a preconditioner

    225:    Options Database Keys:
    226: +  -pc_factor_reuse_ordering - Activate `PCFactorSetReuseOrdering()`
    227: .  -pc_factor_mat_solver_type - Actives `PCFactorSetMatSolverType()` to choose the direct solver, like superlu
    228: .  -pc_factor_reuse_fill - Activates `PCFactorSetReuseFill()`
    229: .  -pc_factor_fill <fill> - Sets fill amount
    230: .  -pc_factor_in_place - Activates in-place factorization
    231: .  -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
    232: .  -pc_factor_pivot_in_blocks <true,false> - allow pivoting within the small blocks during factorization (may increase
    233:                                          stability of factorization.
    234: .  -pc_factor_shift_type <shifttype> - Sets shift type or -1 for the default; use '-help' for a list of available types
    235: .  -pc_factor_shift_amount <shiftamount> - Sets shift amount or -1 for the default
    236: .  -pc_factor_nonzeros_along_diagonal - permutes the rows and columns to try to put nonzero value along the diagonal.
    237: .  -pc_factor_mat_solver_type <packagename> - use an external package for the solve, see `MatSolverType` for possibilities
    238: -  -mat_solvertype_optionname - options for a specific solver package, for example -mat_mumps_cntl_1

    240:    Level: beginner

    242:    Notes:
    243:    Not all options work for all matrix formats

    245:    Run with `-help` to see additional options for particular matrix formats or factorization algorithms

    247:    The Cholesky factorization direct solver, `PCCHOLESKY` will be more efficient than `PCLU` for symmetric positive-definite (SPD) matrices

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

    253: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `MatSolverType`, `MatGetFactor()`, `PCQR`, `PCSVD`,
    254:           `PCILU`, `PCCHOLESKY`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`,
    255:           `PCFactorSetFill()`, `PCFactorSetUseInPlace()`, `PCFactorSetMatOrderingType()`, `PCFactorSetColumnPivot()`,
    256:           `PCFactorSetPivotInBlocks()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`
    257:           `PCFactorReorderForNonzeroDiagonal()`
    258: M*/

    260: PETSC_EXTERN PetscErrorCode PCCreate_LU(PC pc)
    261: {
    262:   PC_LU *dir;

    264:   PetscFunctionBegin;
    265:   PetscCall(PetscNew(&dir));
    266:   pc->data = (void *)dir;
    267:   PetscCall(PCFactorInitialize(pc, MAT_FACTOR_LU));
    268:   dir->nonzerosalongdiagonal = PETSC_FALSE;

    270:   ((PC_Factor *)dir)->info.fill      = 5.0;
    271:   ((PC_Factor *)dir)->info.dtcol     = 1.e-6; /* default to pivoting; this is only thing PETSc LU supports */
    272:   ((PC_Factor *)dir)->info.shifttype = (PetscReal)MAT_SHIFT_NONE;
    273:   dir->col                           = NULL;
    274:   dir->row                           = NULL;

    276:   pc->ops->reset             = PCReset_LU;
    277:   pc->ops->destroy           = PCDestroy_LU;
    278:   pc->ops->apply             = PCApply_LU;
    279:   pc->ops->matapply          = PCMatApply_LU;
    280:   pc->ops->applytranspose    = PCApplyTranspose_LU;
    281:   pc->ops->matapplytranspose = PCMatApplyTranspose_LU;
    282:   pc->ops->setup             = PCSetUp_LU;
    283:   pc->ops->setfromoptions    = PCSetFromOptions_LU;
    284:   pc->ops->view              = PCView_Factor;
    285:   pc->ops->applyrichardson   = NULL;
    286:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorReorderForNonzeroDiagonal_C", PCFactorReorderForNonzeroDiagonal_LU));
    287:   PetscFunctionReturn(PETSC_SUCCESS);
    288: }