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: ilu->nonzerosalongdiagonal = PETSC_TRUE;
12: if (z == PETSC_DECIDE) ilu->nonzerosalongdiagonaltol = 1.e-10;
13: else ilu->nonzerosalongdiagonaltol = z;
14: return 0;
15: }
17: PetscErrorCode PCReset_ILU(PC pc)
18: {
19: PC_ILU *ilu = (PC_ILU *)pc->data;
21: if (!ilu->hdr.inplace) MatDestroy(&((PC_Factor *)ilu)->fact);
22: if (ilu->row && ilu->col && ilu->row != ilu->col) ISDestroy(&ilu->row);
23: ISDestroy(&ilu->col);
24: return 0;
25: }
27: PetscErrorCode PCFactorSetDropTolerance_ILU(PC pc, PetscReal dt, PetscReal dtcol, PetscInt dtcount)
28: {
29: PC_ILU *ilu = (PC_ILU *)pc->data;
31: if (pc->setupcalled && (((PC_Factor *)ilu)->info.dt != dt || ((PC_Factor *)ilu)->info.dtcol != dtcol || ((PC_Factor *)ilu)->info.dtcount != dtcount)) {
32: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot change drop tolerance after using PC");
33: }
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: return 0;
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: PetscOptionsHeadBegin(PetscOptionsObject, "ILU Options");
49: PCSetFromOptions_Factor(pc, PetscOptionsObject);
51: PetscOptionsInt("-pc_factor_levels", "levels of fill", "PCFactorSetLevels", (PetscInt)((PC_Factor *)ilu)->info.levels, &itmp, &flg);
52: if (flg) ((PC_Factor *)ilu)->info.levels = itmp;
54: PetscOptionsBool("-pc_factor_diagonal_fill", "Allow fill into empty diagonal entry", "PCFactorSetAllowDiagonalFill", ((PC_Factor *)ilu)->info.diagonal_fill ? PETSC_TRUE : PETSC_FALSE, &flg, &set);
55: if (set) ((PC_Factor *)ilu)->info.diagonal_fill = (PetscReal)flg;
56: PetscOptionsName("-pc_factor_nonzeros_along_diagonal", "Reorder to remove zeros from diagonal", "PCFactorReorderForNonzeroDiagonal", &flg);
57: if (flg) {
58: tol = PETSC_DECIDE;
59: PetscOptionsReal("-pc_factor_nonzeros_along_diagonal", "Reorder to remove zeros from diagonal", "PCFactorReorderForNonzeroDiagonal", ilu->nonzerosalongdiagonaltol, &tol, NULL);
60: PCFactorReorderForNonzeroDiagonal(pc, tol);
61: }
63: PetscOptionsHeadEnd();
64: return 0;
65: }
67: static PetscErrorCode PCSetUp_ILU(PC pc)
68: {
69: PC_ILU *ilu = (PC_ILU *)pc->data;
70: MatInfo info;
71: PetscBool flg;
72: MatSolverType stype;
73: MatFactorError err;
74: const char *prefix;
76: pc->failedreason = PC_NOERROR;
77: /* ugly hack to change default, since it is not support by some matrix types */
78: if (((PC_Factor *)ilu)->info.shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
79: PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQAIJ, &flg);
80: if (!flg) {
81: PetscObjectTypeCompare((PetscObject)pc->pmat, MATMPIAIJ, &flg);
82: if (!flg) {
83: ((PC_Factor *)ilu)->info.shifttype = (PetscReal)MAT_SHIFT_INBLOCKS;
84: PetscInfo(pc, "Changing shift type from NONZERO to INBLOCKS because block matrices do not support NONZERO\n");
85: }
86: }
87: }
89: PCGetOptionsPrefix(pc, &prefix);
90: MatSetOptionsPrefixFactor(pc->pmat, prefix);
92: MatSetErrorIfFailure(pc->pmat, pc->erroriffailure);
93: if (ilu->hdr.inplace) {
94: if (!pc->setupcalled) {
95: /* In-place factorization only makes sense with the natural ordering,
96: so we only need to get the ordering once, even if nonzero structure changes */
97: /* Should not get the ordering if the factorization routine does not use it, but do not yet have access to the factor matrix */
98: PCFactorSetDefaultOrdering_Factor(pc);
99: MatDestroy(&((PC_Factor *)ilu)->fact);
100: MatGetOrdering(pc->pmat, ((PC_Factor *)ilu)->ordering, &ilu->row, &ilu->col);
101: }
103: /* In place ILU only makes sense with fill factor of 1.0 because
104: cannot have levels of fill */
105: ((PC_Factor *)ilu)->info.fill = 1.0;
106: ((PC_Factor *)ilu)->info.diagonal_fill = 0.0;
108: MatILUFactor(pc->pmat, ilu->row, ilu->col, &((PC_Factor *)ilu)->info);
109: MatFactorGetError(pc->pmat, &err);
110: if (err) { /* Factor() fails */
111: pc->failedreason = (PCFailedReason)err;
112: return 0;
113: }
115: ((PC_Factor *)ilu)->fact = pc->pmat;
116: /* must update the pc record of the matrix state or the PC will attempt to run PCSetUp() yet again */
117: PetscObjectStateGet((PetscObject)pc->pmat, &pc->matstate);
118: } else {
119: if (!pc->setupcalled) {
120: /* first time in so compute reordering and symbolic factorization */
121: PetscBool canuseordering;
122: if (!((PC_Factor *)ilu)->fact) { MatGetFactor(pc->pmat, ((PC_Factor *)ilu)->solvertype, MAT_FACTOR_ILU, &((PC_Factor *)ilu)->fact); }
123: MatFactorGetCanUseOrdering(((PC_Factor *)ilu)->fact, &canuseordering);
124: if (canuseordering) {
125: PCFactorSetDefaultOrdering_Factor(pc);
126: MatGetOrdering(pc->pmat, ((PC_Factor *)ilu)->ordering, &ilu->row, &ilu->col);
127: /* Remove zeros along diagonal? */
128: if (ilu->nonzerosalongdiagonal) MatReorderForNonzeroDiagonal(pc->pmat, ilu->nonzerosalongdiagonaltol, ilu->row, ilu->col);
129: }
130: MatILUFactorSymbolic(((PC_Factor *)ilu)->fact, pc->pmat, ilu->row, ilu->col, &((PC_Factor *)ilu)->info);
131: MatGetInfo(((PC_Factor *)ilu)->fact, MAT_LOCAL, &info);
132: ilu->hdr.actualfill = info.fill_ratio_needed;
133: } else if (pc->flag != SAME_NONZERO_PATTERN) {
134: if (!ilu->hdr.reuseordering) {
135: PetscBool canuseordering;
136: MatDestroy(&((PC_Factor *)ilu)->fact);
137: MatGetFactor(pc->pmat, ((PC_Factor *)ilu)->solvertype, MAT_FACTOR_ILU, &((PC_Factor *)ilu)->fact);
138: MatFactorGetCanUseOrdering(((PC_Factor *)ilu)->fact, &canuseordering);
139: if (canuseordering) {
140: /* compute a new ordering for the ILU */
141: ISDestroy(&ilu->row);
142: ISDestroy(&ilu->col);
143: PCFactorSetDefaultOrdering_Factor(pc);
144: MatGetOrdering(pc->pmat, ((PC_Factor *)ilu)->ordering, &ilu->row, &ilu->col);
145: /* Remove zeros along diagonal? */
146: if (ilu->nonzerosalongdiagonal) MatReorderForNonzeroDiagonal(pc->pmat, ilu->nonzerosalongdiagonaltol, ilu->row, ilu->col);
147: }
148: }
149: MatILUFactorSymbolic(((PC_Factor *)ilu)->fact, pc->pmat, ilu->row, ilu->col, &((PC_Factor *)ilu)->info);
150: MatGetInfo(((PC_Factor *)ilu)->fact, MAT_LOCAL, &info);
151: ilu->hdr.actualfill = info.fill_ratio_needed;
152: }
153: MatFactorGetError(((PC_Factor *)ilu)->fact, &err);
154: if (err) { /* FactorSymbolic() fails */
155: pc->failedreason = (PCFailedReason)err;
156: return 0;
157: }
159: MatLUFactorNumeric(((PC_Factor *)ilu)->fact, pc->pmat, &((PC_Factor *)ilu)->info);
160: MatFactorGetError(((PC_Factor *)ilu)->fact, &err);
161: if (err) { /* FactorNumeric() fails */
162: pc->failedreason = (PCFailedReason)err;
163: }
164: }
166: PCFactorGetMatSolverType(pc, &stype);
167: if (!stype) {
168: MatSolverType solverpackage;
169: MatFactorGetSolverType(((PC_Factor *)ilu)->fact, &solverpackage);
170: PCFactorSetMatSolverType(pc, solverpackage);
171: }
172: return 0;
173: }
175: static PetscErrorCode PCDestroy_ILU(PC pc)
176: {
177: PC_ILU *ilu = (PC_ILU *)pc->data;
179: PCReset_ILU(pc);
180: PetscFree(((PC_Factor *)ilu)->solvertype);
181: PetscFree(((PC_Factor *)ilu)->ordering);
182: PetscFree(pc->data);
183: PCFactorClearComposedFunctions(pc);
184: return 0;
185: }
187: static PetscErrorCode PCApply_ILU(PC pc, Vec x, Vec y)
188: {
189: PC_ILU *ilu = (PC_ILU *)pc->data;
191: MatSolve(((PC_Factor *)ilu)->fact, x, y);
192: return 0;
193: }
195: static PetscErrorCode PCMatApply_ILU(PC pc, Mat X, Mat Y)
196: {
197: PC_ILU *ilu = (PC_ILU *)pc->data;
199: MatMatSolve(((PC_Factor *)ilu)->fact, X, Y);
200: return 0;
201: }
203: static PetscErrorCode PCApplyTranspose_ILU(PC pc, Vec x, Vec y)
204: {
205: PC_ILU *ilu = (PC_ILU *)pc->data;
207: MatSolveTranspose(((PC_Factor *)ilu)->fact, x, y);
208: return 0;
209: }
211: static PetscErrorCode PCApplySymmetricLeft_ILU(PC pc, Vec x, Vec y)
212: {
213: PC_ILU *icc = (PC_ILU *)pc->data;
215: MatForwardSolve(((PC_Factor *)icc)->fact, x, y);
216: return 0;
217: }
219: static PetscErrorCode PCApplySymmetricRight_ILU(PC pc, Vec x, Vec y)
220: {
221: PC_ILU *icc = (PC_ILU *)pc->data;
223: MatBackwardSolve(((PC_Factor *)icc)->fact, x, y);
224: return 0;
225: }
227: /*MC
228: PCILU - Incomplete factorization preconditioners.
230: Options Database Keys:
231: + -pc_factor_levels <k> - number of levels of fill for ILU(k)
232: . -pc_factor_in_place - only for ILU(0) with natural ordering, reuses the space of the matrix for
233: its factorization (overwrites original matrix)
234: . -pc_factor_diagonal_fill - fill in a zero diagonal even if levels of fill indicate it wouldn't be fill
235: . -pc_factor_reuse_ordering - reuse ordering of factorized matrix from previous factorization
236: . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
237: . -pc_factor_nonzeros_along_diagonal - reorder the matrix before factorization to remove zeros from the diagonal,
238: this decreases the chance of getting a zero pivot
239: . -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
240: - -pc_factor_pivot_in_blocks - for block ILU(k) factorization, i.e. with BAIJ matrices with block size larger
241: than 1 the diagonal blocks are factored with partial pivoting (this increases the
242: stability of the ILU factorization
244: Level: beginner
246: Notes:
247: Only implemented for some matrix format and sequential. For parallel see `PCHYPRE` for hypre's ILU
249: For `MATSEQBAIJ` matrices this implements a point block ILU
251: The "symmetric" application of this preconditioner is not actually symmetric since L is not transpose(U)
252: even when the matrix is not symmetric since the U stores the diagonals of the factorization.
254: If you are using `MATSEQAIJCUSPARSE` matrices (or `MATMPIAIJCUSPARSE` matrices with block Jacobi), factorization
255: is never done on the GPU).
257: References:
258: + * - T. Dupont, R. Kendall, and H. Rachford. An approximate factorization procedure for solving
259: self adjoint elliptic difference equations. SIAM J. Numer. Anal., 5, 1968.
260: . * - T.A. Oliphant. An implicit numerical method for solving two dimensional timedependent diffusion problems. Quart. Appl. Math., 19, 1961.
261: - * - TONY F. CHAN AND HENK A. VAN DER VORST, APPROXIMATE AND INCOMPLETE FACTORIZATIONS,
262: Chapter in Parallel Numerical
263: Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
264: Science and Engineering, Kluwer.
266: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSOR`, `MatOrderingType`, `PCLU`, `PCICC`, `PCCHOLESKY`,
267: `PCFactorSetZeroPivot()`, `PCFactorSetShiftSetType()`, `PCFactorSetAmount()`,
268: `PCFactorSetDropTolerance()`, `PCFactorSetFill()`, `PCFactorSetMatOrderingType()`, `PCFactorSetReuseOrdering()`,
269: `PCFactorSetLevels()`, `PCFactorSetUseInPlace()`, `PCFactorSetAllowDiagonalFill()`, `PCFactorSetPivotInBlocks()`,
270: `PCFactorGetAllowDiagonalFill()`, `PCFactorGetUseInPlace()`
271: M*/
273: PETSC_EXTERN PetscErrorCode PCCreate_ILU(PC pc)
274: {
275: PC_ILU *ilu;
277: PetscNew(&ilu);
278: pc->data = (void *)ilu;
279: PCFactorInitialize(pc, MAT_FACTOR_ILU);
281: ((PC_Factor *)ilu)->info.levels = 0.;
282: ((PC_Factor *)ilu)->info.fill = 1.0;
283: ilu->col = NULL;
284: ilu->row = NULL;
285: ((PC_Factor *)ilu)->info.dt = PETSC_DEFAULT;
286: ((PC_Factor *)ilu)->info.dtcount = PETSC_DEFAULT;
287: ((PC_Factor *)ilu)->info.dtcol = PETSC_DEFAULT;
289: pc->ops->reset = PCReset_ILU;
290: pc->ops->destroy = PCDestroy_ILU;
291: pc->ops->apply = PCApply_ILU;
292: pc->ops->matapply = PCMatApply_ILU;
293: pc->ops->applytranspose = PCApplyTranspose_ILU;
294: pc->ops->setup = PCSetUp_ILU;
295: pc->ops->setfromoptions = PCSetFromOptions_ILU;
296: pc->ops->view = PCView_Factor;
297: pc->ops->applysymmetricleft = PCApplySymmetricLeft_ILU;
298: pc->ops->applysymmetricright = PCApplySymmetricRight_ILU;
299: pc->ops->applyrichardson = NULL;
300: PetscObjectComposeFunction((PetscObject)pc, "PCFactorSetDropTolerance_C", PCFactorSetDropTolerance_ILU);
301: PetscObjectComposeFunction((PetscObject)pc, "PCFactorReorderForNonzeroDiagonal_C", PCFactorReorderForNonzeroDiagonal_ILU);
302: return 0;
303: }