1 /*
2 Defines a ILU factorization preconditioner for any Mat implementation
3 */
4 #include <../src/ksp/pc/impls/factor/ilu/ilu.h> /*I "petscpc.h" I*/
5
PCFactorReorderForNonzeroDiagonal_ILU(PC pc,PetscReal z)6 static PetscErrorCode PCFactorReorderForNonzeroDiagonal_ILU(PC pc, PetscReal z)
7 {
8 PC_ILU *ilu = (PC_ILU *)pc->data;
9
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 }
16
PCReset_ILU(PC pc)17 static PetscErrorCode PCReset_ILU(PC pc)
18 {
19 PC_ILU *ilu = (PC_ILU *)pc->data;
20
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 }
27
PCFactorSetDropTolerance_ILU(PC pc,PetscReal dt,PetscReal dtcol,PetscInt dtcount)28 static PetscErrorCode PCFactorSetDropTolerance_ILU(PC pc, PetscReal dt, PetscReal dtcol, PetscInt dtcount)
29 {
30 PC_ILU *ilu = (PC_ILU *)pc->data;
31
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 }
40
PCSetFromOptions_ILU(PC pc,PetscOptionItems PetscOptionsObject)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;
47
48 PetscFunctionBegin;
49 PetscOptionsHeadBegin(PetscOptionsObject, "ILU Options");
50 PetscCall(PCSetFromOptions_Factor(pc, PetscOptionsObject));
51
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;
54
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 }
63
64 PetscOptionsHeadEnd();
65 PetscFunctionReturn(PETSC_SUCCESS);
66 }
67
PCSetUp_ILU(PC pc)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;
76
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 }
90
91 PetscCall(PCGetOptionsPrefix(pc, &prefix));
92 PetscCall(MatSetOptionsPrefixFactor(pc->pmat, prefix));
93
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 }
104
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;
109
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 }
116
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;
124
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 PetscCall(MatDestroy(&((PC_Factor *)ilu)->fact));
138 PetscCall(PCFactorSetUpMatSolverType(pc));
139 if (!ilu->hdr.reuseordering) {
140 PetscBool canuseordering;
141
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 }
162
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 }
169
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 }
178
PCDestroy_ILU(PC pc)179 static PetscErrorCode PCDestroy_ILU(PC pc)
180 {
181 PC_ILU *ilu = (PC_ILU *)pc->data;
182
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 }
191
PCApply_ILU(PC pc,Vec x,Vec y)192 static PetscErrorCode PCApply_ILU(PC pc, Vec x, Vec y)
193 {
194 PC_ILU *ilu = (PC_ILU *)pc->data;
195
196 PetscFunctionBegin;
197 PetscCall(MatSolve(((PC_Factor *)ilu)->fact, x, y));
198 PetscFunctionReturn(PETSC_SUCCESS);
199 }
200
PCMatApply_ILU(PC pc,Mat X,Mat Y)201 static PetscErrorCode PCMatApply_ILU(PC pc, Mat X, Mat Y)
202 {
203 PC_ILU *ilu = (PC_ILU *)pc->data;
204
205 PetscFunctionBegin;
206 PetscCall(MatMatSolve(((PC_Factor *)ilu)->fact, X, Y));
207 PetscFunctionReturn(PETSC_SUCCESS);
208 }
209
PCApplyTranspose_ILU(PC pc,Vec x,Vec y)210 static PetscErrorCode PCApplyTranspose_ILU(PC pc, Vec x, Vec y)
211 {
212 PC_ILU *ilu = (PC_ILU *)pc->data;
213
214 PetscFunctionBegin;
215 PetscCall(MatSolveTranspose(((PC_Factor *)ilu)->fact, x, y));
216 PetscFunctionReturn(PETSC_SUCCESS);
217 }
218
PCMatApplyTranspose_ILU(PC pc,Mat X,Mat Y)219 static PetscErrorCode PCMatApplyTranspose_ILU(PC pc, Mat X, Mat Y)
220 {
221 PC_ILU *ilu = (PC_ILU *)pc->data;
222
223 PetscFunctionBegin;
224 PetscCall(MatMatSolveTranspose(((PC_Factor *)ilu)->fact, X, Y));
225 PetscFunctionReturn(PETSC_SUCCESS);
226 }
227
PCApplySymmetricLeft_ILU(PC pc,Vec x,Vec y)228 static PetscErrorCode PCApplySymmetricLeft_ILU(PC pc, Vec x, Vec y)
229 {
230 PC_ILU *icc = (PC_ILU *)pc->data;
231
232 PetscFunctionBegin;
233 PetscCall(MatForwardSolve(((PC_Factor *)icc)->fact, x, y));
234 PetscFunctionReturn(PETSC_SUCCESS);
235 }
236
PCApplySymmetricRight_ILU(PC pc,Vec x,Vec y)237 static PetscErrorCode PCApplySymmetricRight_ILU(PC pc, Vec x, Vec y)
238 {
239 PC_ILU *icc = (PC_ILU *)pc->data;
240
241 PetscFunctionBegin;
242 PetscCall(MatBackwardSolve(((PC_Factor *)icc)->fact, x, y));
243 PetscFunctionReturn(PETSC_SUCCESS);
244 }
245
246 /*MC
247 PCILU - Incomplete factorization preconditioners {cite}`dupont1968approximate`, {cite}`oliphant1961implicit`, {cite}`chan1997approximate`
248
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
262
263 Level: beginner
264
265 Notes:
266 Only implemented for some matrix format and sequential. For parallel see `PCHYPRE` for hypre's ILU
267
268 For `MATSEQBAIJ` matrices this implements a point block ILU
269
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.
272
273 If you are using `MATSEQAIJCUSPARSE` matrices (or `MATMPIAIJCUSPARSE` matrices with block Jacobi), factorization
274 is never done on the GPU).
275
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*/
282
PCCreate_ILU(PC pc)283 PETSC_EXTERN PetscErrorCode PCCreate_ILU(PC pc)
284 {
285 PC_ILU *ilu;
286
287 PetscFunctionBegin;
288 PetscCall(PetscNew(&ilu));
289 pc->data = (void *)ilu;
290 PetscCall(PCFactorInitialize(pc, MAT_FACTOR_ILU));
291
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;
299
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 }
316