Actual source code: ilu.c

      1: /*
      2:    Defines a ILU factorization preconditioner for any Mat implementation
      3: */
      4: #include <../src/ksp/pc/impls/factor/ilu/ilu.h>

      6: static PetscErrorCode PCFactorReorderForNonzeroDiagonal_ILU(PC pc, PetscReal z)
      7: {
      8:   PC_ILU *ilu = (PC_ILU *)pc->data;

     10:   PetscFunctionBegin;
     11:   ilu->nonzerosalongdiagonal = PETSC_TRUE;
     12:   if (z == (PetscReal)PETSC_DECIDE) ilu->nonzerosalongdiagonaltol = 1.e-10;
     13:   else ilu->nonzerosalongdiagonaltol = z;
     14:   PetscFunctionReturn(PETSC_SUCCESS);
     15: }

     17: static PetscErrorCode PCReset_ILU(PC pc)
     18: {
     19:   PC_ILU *ilu = (PC_ILU *)pc->data;

     21:   PetscFunctionBegin;
     22:   if (!ilu->hdr.inplace) PetscCall(MatDestroy(&((PC_Factor *)ilu)->fact));
     23:   if (ilu->row && ilu->col && ilu->row != ilu->col) PetscCall(ISDestroy(&ilu->row));
     24:   PetscCall(ISDestroy(&ilu->col));
     25:   PetscFunctionReturn(PETSC_SUCCESS);
     26: }

     28: PetscErrorCode PCFactorSetDropTolerance_ILU(PC pc, PetscReal dt, PetscReal dtcol, PetscInt dtcount)
     29: {
     30:   PC_ILU *ilu = (PC_ILU *)pc->data;

     32:   PetscFunctionBegin;
     33:   if (pc->setupcalled && (((PC_Factor *)ilu)->info.dt != dt || ((PC_Factor *)ilu)->info.dtcol != dtcol || ((PC_Factor *)ilu)->info.dtcount != dtcount)) {
     34:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot change drop tolerance after using PC");
     35:   }
     36:   ((PC_Factor *)ilu)->info.dt      = dt;
     37:   ((PC_Factor *)ilu)->info.dtcol   = dtcol;
     38:   ((PC_Factor *)ilu)->info.dtcount = dtcount;
     39:   ((PC_Factor *)ilu)->info.usedt   = 1.0;
     40:   PetscFunctionReturn(PETSC_SUCCESS);
     41: }

     43: static PetscErrorCode PCSetFromOptions_ILU(PC pc, PetscOptionItems *PetscOptionsObject)
     44: {
     45:   PetscInt  itmp;
     46:   PetscBool flg, set;
     47:   PC_ILU   *ilu = (PC_ILU *)pc->data;
     48:   PetscReal tol;

     50:   PetscFunctionBegin;
     51:   PetscOptionsHeadBegin(PetscOptionsObject, "ILU Options");
     52:   PetscCall(PCSetFromOptions_Factor(pc, PetscOptionsObject));

     54:   PetscCall(PetscOptionsInt("-pc_factor_levels", "levels of fill", "PCFactorSetLevels", (PetscInt)((PC_Factor *)ilu)->info.levels, &itmp, &flg));
     55:   if (flg) ((PC_Factor *)ilu)->info.levels = itmp;

     57:   PetscCall(PetscOptionsBool("-pc_factor_diagonal_fill", "Allow fill into empty diagonal entry", "PCFactorSetAllowDiagonalFill", ((PC_Factor *)ilu)->info.diagonal_fill ? PETSC_TRUE : PETSC_FALSE, &flg, &set));
     58:   if (set) ((PC_Factor *)ilu)->info.diagonal_fill = (PetscReal)flg;
     59:   PetscCall(PetscOptionsName("-pc_factor_nonzeros_along_diagonal", "Reorder to remove zeros from diagonal", "PCFactorReorderForNonzeroDiagonal", &flg));
     60:   if (flg) {
     61:     tol = PETSC_DECIDE;
     62:     PetscCall(PetscOptionsReal("-pc_factor_nonzeros_along_diagonal", "Reorder to remove zeros from diagonal", "PCFactorReorderForNonzeroDiagonal", ilu->nonzerosalongdiagonaltol, &tol, NULL));
     63:     PetscCall(PCFactorReorderForNonzeroDiagonal(pc, tol));
     64:   }

     66:   PetscOptionsHeadEnd();
     67:   PetscFunctionReturn(PETSC_SUCCESS);
     68: }

     70: static PetscErrorCode PCSetUp_ILU(PC pc)
     71: {
     72:   PC_ILU        *ilu = (PC_ILU *)pc->data;
     73:   MatInfo        info;
     74:   PetscBool      flg;
     75:   MatSolverType  stype;
     76:   MatFactorError err;
     77:   const char    *prefix;

     79:   PetscFunctionBegin;
     80:   pc->failedreason = PC_NOERROR;
     81:   /* ugly hack to change default, since it is not support by some matrix types */
     82:   if (((PC_Factor *)ilu)->info.shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
     83:     PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQAIJ, &flg));
     84:     if (!flg) {
     85:       PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATMPIAIJ, &flg));
     86:       if (!flg) {
     87:         ((PC_Factor *)ilu)->info.shifttype = (PetscReal)MAT_SHIFT_INBLOCKS;
     88:         PetscCall(PetscInfo(pc, "Changing shift type from NONZERO to INBLOCKS because block matrices do not support NONZERO\n"));
     89:       }
     90:     }
     91:   }

     93:   PetscCall(PCGetOptionsPrefix(pc, &prefix));
     94:   PetscCall(MatSetOptionsPrefixFactor(pc->pmat, prefix));

     96:   PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
     97:   if (ilu->hdr.inplace) {
     98:     if (!pc->setupcalled) {
     99:       /* In-place factorization only makes sense with the natural ordering,
    100:          so we only need to get the ordering once, even if nonzero structure changes */
    101:       /* Should not get the ordering if the factorization routine does not use it, but do not yet have access to the factor matrix */
    102:       PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
    103:       PetscCall(MatDestroy(&((PC_Factor *)ilu)->fact));
    104:       PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)ilu)->ordering, &ilu->row, &ilu->col));
    105:     }

    107:     /* In place ILU only makes sense with fill factor of 1.0 because
    108:        cannot have levels of fill */
    109:     ((PC_Factor *)ilu)->info.fill          = 1.0;
    110:     ((PC_Factor *)ilu)->info.diagonal_fill = 0.0;

    112:     PetscCall(MatILUFactor(pc->pmat, ilu->row, ilu->col, &((PC_Factor *)ilu)->info));
    113:     PetscCall(MatFactorGetError(pc->pmat, &err));
    114:     if (err) { /* Factor() fails */
    115:       pc->failedreason = (PCFailedReason)err;
    116:       PetscFunctionReturn(PETSC_SUCCESS);
    117:     }

    119:     ((PC_Factor *)ilu)->fact = pc->pmat;
    120:     /* must update the pc record of the matrix state or the PC will attempt to run PCSetUp() yet again */
    121:     PetscCall(PetscObjectStateGet((PetscObject)pc->pmat, &pc->matstate));
    122:   } else {
    123:     if (!pc->setupcalled) {
    124:       /* first time in so compute reordering and symbolic factorization */
    125:       PetscBool canuseordering;

    127:       PetscCall(PCFactorSetUpMatSolverType(pc));
    128:       PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)ilu)->fact, &canuseordering));
    129:       if (canuseordering) {
    130:         PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
    131:         PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)ilu)->ordering, &ilu->row, &ilu->col));
    132:         /*  Remove zeros along diagonal?     */
    133:         if (ilu->nonzerosalongdiagonal) PetscCall(MatReorderForNonzeroDiagonal(pc->pmat, ilu->nonzerosalongdiagonaltol, ilu->row, ilu->col));
    134:       }
    135:       PetscCall(MatILUFactorSymbolic(((PC_Factor *)ilu)->fact, pc->pmat, ilu->row, ilu->col, &((PC_Factor *)ilu)->info));
    136:       PetscCall(MatGetInfo(((PC_Factor *)ilu)->fact, MAT_LOCAL, &info));
    137:       ilu->hdr.actualfill = info.fill_ratio_needed;
    138:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
    139:       if (!ilu->hdr.reuseordering) {
    140:         PetscBool canuseordering;

    142:         PetscCall(MatDestroy(&((PC_Factor *)ilu)->fact));
    143:         PetscCall(PCFactorSetUpMatSolverType(pc));
    144:         PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)ilu)->fact, &canuseordering));
    145:         if (canuseordering) {
    146:           /* compute a new ordering for the ILU */
    147:           PetscCall(ISDestroy(&ilu->row));
    148:           PetscCall(ISDestroy(&ilu->col));
    149:           PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
    150:           PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)ilu)->ordering, &ilu->row, &ilu->col));
    151:           /*  Remove zeros along diagonal?     */
    152:           if (ilu->nonzerosalongdiagonal) PetscCall(MatReorderForNonzeroDiagonal(pc->pmat, ilu->nonzerosalongdiagonaltol, ilu->row, ilu->col));
    153:         }
    154:       }
    155:       PetscCall(MatILUFactorSymbolic(((PC_Factor *)ilu)->fact, pc->pmat, ilu->row, ilu->col, &((PC_Factor *)ilu)->info));
    156:       PetscCall(MatGetInfo(((PC_Factor *)ilu)->fact, MAT_LOCAL, &info));
    157:       ilu->hdr.actualfill = info.fill_ratio_needed;
    158:     }
    159:     PetscCall(MatFactorGetError(((PC_Factor *)ilu)->fact, &err));
    160:     if (err) { /* FactorSymbolic() fails */
    161:       pc->failedreason = (PCFailedReason)err;
    162:       PetscFunctionReturn(PETSC_SUCCESS);
    163:     }

    165:     PetscCall(MatLUFactorNumeric(((PC_Factor *)ilu)->fact, pc->pmat, &((PC_Factor *)ilu)->info));
    166:     PetscCall(MatFactorGetError(((PC_Factor *)ilu)->fact, &err));
    167:     if (err) { /* FactorNumeric() fails */
    168:       pc->failedreason = (PCFailedReason)err;
    169:     }
    170:   }

    172:   PetscCall(PCFactorGetMatSolverType(pc, &stype));
    173:   if (!stype) {
    174:     MatSolverType solverpackage;
    175:     PetscCall(MatFactorGetSolverType(((PC_Factor *)ilu)->fact, &solverpackage));
    176:     PetscCall(PCFactorSetMatSolverType(pc, solverpackage));
    177:   }
    178:   PetscFunctionReturn(PETSC_SUCCESS);
    179: }

    181: static PetscErrorCode PCDestroy_ILU(PC pc)
    182: {
    183:   PC_ILU *ilu = (PC_ILU *)pc->data;

    185:   PetscFunctionBegin;
    186:   PetscCall(PCReset_ILU(pc));
    187:   PetscCall(PetscFree(((PC_Factor *)ilu)->solvertype));
    188:   PetscCall(PetscFree(((PC_Factor *)ilu)->ordering));
    189:   PetscCall(PetscFree(pc->data));
    190:   PetscCall(PCFactorClearComposedFunctions(pc));
    191:   PetscFunctionReturn(PETSC_SUCCESS);
    192: }

    194: static PetscErrorCode PCApply_ILU(PC pc, Vec x, Vec y)
    195: {
    196:   PC_ILU *ilu = (PC_ILU *)pc->data;

    198:   PetscFunctionBegin;
    199:   PetscCall(MatSolve(((PC_Factor *)ilu)->fact, x, y));
    200:   PetscFunctionReturn(PETSC_SUCCESS);
    201: }

    203: static PetscErrorCode PCMatApply_ILU(PC pc, Mat X, Mat Y)
    204: {
    205:   PC_ILU *ilu = (PC_ILU *)pc->data;

    207:   PetscFunctionBegin;
    208:   PetscCall(MatMatSolve(((PC_Factor *)ilu)->fact, X, Y));
    209:   PetscFunctionReturn(PETSC_SUCCESS);
    210: }

    212: static PetscErrorCode PCApplyTranspose_ILU(PC pc, Vec x, Vec y)
    213: {
    214:   PC_ILU *ilu = (PC_ILU *)pc->data;

    216:   PetscFunctionBegin;
    217:   PetscCall(MatSolveTranspose(((PC_Factor *)ilu)->fact, x, y));
    218:   PetscFunctionReturn(PETSC_SUCCESS);
    219: }

    221: static PetscErrorCode PCApplySymmetricLeft_ILU(PC pc, Vec x, Vec y)
    222: {
    223:   PC_ILU *icc = (PC_ILU *)pc->data;

    225:   PetscFunctionBegin;
    226:   PetscCall(MatForwardSolve(((PC_Factor *)icc)->fact, x, y));
    227:   PetscFunctionReturn(PETSC_SUCCESS);
    228: }

    230: static PetscErrorCode PCApplySymmetricRight_ILU(PC pc, Vec x, Vec y)
    231: {
    232:   PC_ILU *icc = (PC_ILU *)pc->data;

    234:   PetscFunctionBegin;
    235:   PetscCall(MatBackwardSolve(((PC_Factor *)icc)->fact, x, y));
    236:   PetscFunctionReturn(PETSC_SUCCESS);
    237: }

    239: /*MC
    240:      PCILU - Incomplete factorization preconditioners {cite}`dupont1968approximate`, {cite}`oliphant1961implicit`, {cite}`chan1997approximate`

    242:    Options Database Keys:
    243: +  -pc_factor_levels <k>                                 - number of levels of fill for ILU(k)
    244: .  -pc_factor_in_place                                   - only for ILU(0) with natural ordering, reuses the space of the matrix for
    245:                                                          its factorization (overwrites original matrix)
    246: .  -pc_factor_diagonal_fill                              - fill in a zero diagonal even if levels of fill indicate it wouldn't be fill
    247: .  -pc_factor_reuse_ordering                             - reuse ordering of factorized matrix from previous factorization
    248: .  -pc_factor_fill <nfill>                               - expected amount of fill in factored matrix compared to original matrix, nfill > 1
    249: .  -pc_factor_nonzeros_along_diagonal                    - reorder the matrix before factorization to remove zeros from the diagonal,
    250:                                                          this decreases the chance of getting a zero pivot
    251: .  -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
    252: -  -pc_factor_pivot_in_blocks                            - for block ILU(k) factorization, i.e. with `MATBAIJ` matrices with block size larger
    253:                                                          than 1 the diagonal blocks are factored with partial pivoting (this increases the
    254:                                                          stability of the ILU factorization

    256:    Level: beginner

    258:    Notes:
    259:    Only implemented for some matrix format and sequential. For parallel see `PCHYPRE` for hypre's ILU

    261:    For `MATSEQBAIJ` matrices this implements a point block ILU

    263:    The "symmetric" application of this preconditioner is not actually symmetric since L is not transpose(U)
    264:    even when the matrix is not symmetric since the U stores the diagonals of the factorization.

    266:    If you are using `MATSEQAIJCUSPARSE` matrices (or `MATMPIAIJCUSPARSE` matrices with block Jacobi), factorization
    267:    is never done on the GPU).

    269: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSOR`, `MatOrderingType`, `PCLU`, `PCICC`, `PCCHOLESKY`,
    270:           `PCFactorSetZeroPivot()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`,
    271:           `PCFactorSetDropTolerance()`, `PCFactorSetFill()`, `PCFactorSetMatOrderingType()`, `PCFactorSetReuseOrdering()`,
    272:           `PCFactorSetLevels()`, `PCFactorSetUseInPlace()`, `PCFactorSetAllowDiagonalFill()`, `PCFactorSetPivotInBlocks()`,
    273:           `PCFactorGetAllowDiagonalFill()`, `PCFactorGetUseInPlace()`
    274: M*/

    276: PETSC_EXTERN PetscErrorCode PCCreate_ILU(PC pc)
    277: {
    278:   PC_ILU *ilu;

    280:   PetscFunctionBegin;
    281:   PetscCall(PetscNew(&ilu));
    282:   pc->data = (void *)ilu;
    283:   PetscCall(PCFactorInitialize(pc, MAT_FACTOR_ILU));

    285:   ((PC_Factor *)ilu)->info.levels  = 0.;
    286:   ((PC_Factor *)ilu)->info.fill    = 1.0;
    287:   ilu->col                         = NULL;
    288:   ilu->row                         = NULL;
    289:   ((PC_Factor *)ilu)->info.dt      = PETSC_DEFAULT;
    290:   ((PC_Factor *)ilu)->info.dtcount = PETSC_DEFAULT;
    291:   ((PC_Factor *)ilu)->info.dtcol   = PETSC_DEFAULT;

    293:   pc->ops->reset               = PCReset_ILU;
    294:   pc->ops->destroy             = PCDestroy_ILU;
    295:   pc->ops->apply               = PCApply_ILU;
    296:   pc->ops->matapply            = PCMatApply_ILU;
    297:   pc->ops->applytranspose      = PCApplyTranspose_ILU;
    298:   pc->ops->setup               = PCSetUp_ILU;
    299:   pc->ops->setfromoptions      = PCSetFromOptions_ILU;
    300:   pc->ops->view                = PCView_Factor;
    301:   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_ILU;
    302:   pc->ops->applysymmetricright = PCApplySymmetricRight_ILU;
    303:   pc->ops->applyrichardson     = NULL;
    304:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorSetDropTolerance_C", PCFactorSetDropTolerance_ILU));
    305:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorReorderForNonzeroDiagonal_C", PCFactorReorderForNonzeroDiagonal_ILU));
    306:   PetscFunctionReturn(PETSC_SUCCESS);
    307: }