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: static PetscErrorCode PCFactorSetDropTolerance_ILU(PC pc, PetscReal dt, PetscReal dtcol, PetscInt dtcount)
 29: {
 30:   PC_ILU *ilu = (PC_ILU *)pc->data;

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

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

 48:   PetscFunctionBegin;
 49:   PetscOptionsHeadBegin(PetscOptionsObject, "ILU Options");
 50:   PetscCall(PCSetFromOptions_Factor(pc, PetscOptionsObject));

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

 55:   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));
 56:   if (set) ((PC_Factor *)ilu)->info.diagonal_fill = (PetscReal)flg;
 57:   PetscCall(PetscOptionsName("-pc_factor_nonzeros_along_diagonal", "Reorder to remove zeros from diagonal", "PCFactorReorderForNonzeroDiagonal", &flg));
 58:   if (flg) {
 59:     tol = PETSC_DECIDE;
 60:     PetscCall(PetscOptionsReal("-pc_factor_nonzeros_along_diagonal", "Reorder to remove zeros from diagonal", "PCFactorReorderForNonzeroDiagonal", ilu->nonzerosalongdiagonaltol, &tol, NULL));
 61:     PetscCall(PCFactorReorderForNonzeroDiagonal(pc, tol));
 62:   }

 64:   PetscOptionsHeadEnd();
 65:   PetscFunctionReturn(PETSC_SUCCESS);
 66: }

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

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

 91:   PetscCall(PCGetOptionsPrefix(pc, &prefix));
 92:   PetscCall(MatSetOptionsPrefixFactor(pc->pmat, prefix));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

219: static PetscErrorCode PCMatApplyTranspose_ILU(PC pc, Mat X, Mat Y)
220: {
221:   PC_ILU *ilu = (PC_ILU *)pc->data;

223:   PetscFunctionBegin;
224:   PetscCall(MatMatSolveTranspose(((PC_Factor *)ilu)->fact, X, Y));
225:   PetscFunctionReturn(PETSC_SUCCESS);
226: }

228: static PetscErrorCode PCApplySymmetricLeft_ILU(PC pc, Vec x, Vec y)
229: {
230:   PC_ILU *icc = (PC_ILU *)pc->data;

232:   PetscFunctionBegin;
233:   PetscCall(MatForwardSolve(((PC_Factor *)icc)->fact, x, y));
234:   PetscFunctionReturn(PETSC_SUCCESS);
235: }

237: static PetscErrorCode PCApplySymmetricRight_ILU(PC pc, Vec x, Vec y)
238: {
239:   PC_ILU *icc = (PC_ILU *)pc->data;

241:   PetscFunctionBegin;
242:   PetscCall(MatBackwardSolve(((PC_Factor *)icc)->fact, x, y));
243:   PetscFunctionReturn(PETSC_SUCCESS);
244: }

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

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

263:    Level: beginner

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

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

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

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

276: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSOR`, `MatOrderingType`, `PCLU`, `PCICC`, `PCCHOLESKY`,
277:           `PCFactorSetZeroPivot()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`,
278:           `PCFactorSetDropTolerance()`, `PCFactorSetFill()`, `PCFactorSetMatOrderingType()`, `PCFactorSetReuseOrdering()`,
279:           `PCFactorSetLevels()`, `PCFactorSetUseInPlace()`, `PCFactorSetAllowDiagonalFill()`, `PCFactorSetPivotInBlocks()`,
280:           `PCFactorGetAllowDiagonalFill()`, `PCFactorGetUseInPlace()`
281: M*/

283: PETSC_EXTERN PetscErrorCode PCCreate_ILU(PC pc)
284: {
285:   PC_ILU *ilu;

287:   PetscFunctionBegin;
288:   PetscCall(PetscNew(&ilu));
289:   pc->data = (void *)ilu;
290:   PetscCall(PCFactorInitialize(pc, MAT_FACTOR_ILU));

292:   ((PC_Factor *)ilu)->info.levels  = 0.;
293:   ((PC_Factor *)ilu)->info.fill    = 1.0;
294:   ilu->col                         = NULL;
295:   ilu->row                         = NULL;
296:   ((PC_Factor *)ilu)->info.dt      = PETSC_DEFAULT;
297:   ((PC_Factor *)ilu)->info.dtcount = PETSC_DEFAULT;
298:   ((PC_Factor *)ilu)->info.dtcol   = PETSC_DEFAULT;

300:   pc->ops->reset               = PCReset_ILU;
301:   pc->ops->destroy             = PCDestroy_ILU;
302:   pc->ops->apply               = PCApply_ILU;
303:   pc->ops->matapply            = PCMatApply_ILU;
304:   pc->ops->applytranspose      = PCApplyTranspose_ILU;
305:   pc->ops->matapplytranspose   = PCMatApplyTranspose_ILU;
306:   pc->ops->setup               = PCSetUp_ILU;
307:   pc->ops->setfromoptions      = PCSetFromOptions_ILU;
308:   pc->ops->view                = PCView_Factor;
309:   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_ILU;
310:   pc->ops->applysymmetricright = PCApplySymmetricRight_ILU;
311:   pc->ops->applyrichardson     = NULL;
312:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorSetDropTolerance_C", PCFactorSetDropTolerance_ILU));
313:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorReorderForNonzeroDiagonal_C", PCFactorReorderForNonzeroDiagonal_ILU));
314:   PetscFunctionReturn(PETSC_SUCCESS);
315: }