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