Actual source code: ilu.c


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

120:     ((PC_Factor *)ilu)->fact = pc->pmat;
121:     /* must update the pc record of the matrix state or the PC will attempt to run PCSetUp() yet again */
122:     PetscCall(PetscObjectStateGet((PetscObject)pc->pmat, &pc->matstate));
123:   } else {
124:     if (!pc->setupcalled) {
125:       /* first time in so compute reordering and symbolic factorization */
126:       PetscBool canuseordering;
127:       if (!((PC_Factor *)ilu)->fact) { PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)ilu)->solvertype, MAT_FACTOR_ILU, &((PC_Factor *)ilu)->fact)); }
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;
141:         PetscCall(MatDestroy(&((PC_Factor *)ilu)->fact));
142:         PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)ilu)->solvertype, MAT_FACTOR_ILU, &((PC_Factor *)ilu)->fact));
143:         PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)ilu)->fact, &canuseordering));
144:         if (canuseordering) {
145:           /* compute a new ordering for the ILU */
146:           PetscCall(ISDestroy(&ilu->row));
147:           PetscCall(ISDestroy(&ilu->col));
148:           PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
149:           PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)ilu)->ordering, &ilu->row, &ilu->col));
150:           /*  Remove zeros along diagonal?     */
151:           if (ilu->nonzerosalongdiagonal) PetscCall(MatReorderForNonzeroDiagonal(pc->pmat, ilu->nonzerosalongdiagonaltol, ilu->row, ilu->col));
152:         }
153:       }
154:       PetscCall(MatILUFactorSymbolic(((PC_Factor *)ilu)->fact, pc->pmat, ilu->row, ilu->col, &((PC_Factor *)ilu)->info));
155:       PetscCall(MatGetInfo(((PC_Factor *)ilu)->fact, MAT_LOCAL, &info));
156:       ilu->hdr.actualfill = info.fill_ratio_needed;
157:     }
158:     PetscCall(MatFactorGetError(((PC_Factor *)ilu)->fact, &err));
159:     if (err) { /* FactorSymbolic() fails */
160:       pc->failedreason = (PCFailedReason)err;
161:       PetscFunctionReturn(PETSC_SUCCESS);
162:     }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

238: /*MC
239:      PCILU - Incomplete factorization preconditioners.

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

255:    Level: beginner

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

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

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

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

268:    References:
269: +  * - T. Dupont, R. Kendall, and H. Rachford. An approximate factorization procedure for solving
270:    self adjoint elliptic difference equations. SIAM J. Numer. Anal., 5, 1968.
271: .  * -  T.A. Oliphant. An implicit numerical method for solving two dimensional timedependent diffusion problems. Quart. Appl. Math., 19, 1961.
272: -  * -  TONY F. CHAN AND HENK A. VAN DER VORST, APPROXIMATE AND INCOMPLETE FACTORIZATIONS,
273:       Chapter in Parallel Numerical
274:       Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
275:       Science and Engineering, Kluwer.

277: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSOR`, `MatOrderingType`, `PCLU`, `PCICC`, `PCCHOLESKY`,
278:           `PCFactorSetZeroPivot()`, `PCFactorSetShiftSetType()`, `PCFactorSetAmount()`,
279:           `PCFactorSetDropTolerance()`, `PCFactorSetFill()`, `PCFactorSetMatOrderingType()`, `PCFactorSetReuseOrdering()`,
280:           `PCFactorSetLevels()`, `PCFactorSetUseInPlace()`, `PCFactorSetAllowDiagonalFill()`, `PCFactorSetPivotInBlocks()`,
281:           `PCFactorGetAllowDiagonalFill()`, `PCFactorGetUseInPlace()`
282: M*/

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

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

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

301:   pc->ops->reset               = PCReset_ILU;
302:   pc->ops->destroy             = PCDestroy_ILU;
303:   pc->ops->apply               = PCApply_ILU;
304:   pc->ops->matapply            = PCMatApply_ILU;
305:   pc->ops->applytranspose      = PCApplyTranspose_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: }