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: }