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(PetscOptionItems *PetscOptionsObject,PC pc)
42: {
43: PetscInt itmp;
44: PetscBool flg,set;
45: PC_ILU *ilu = (PC_ILU*)pc->data;
46: PetscReal tol;
48: PetscOptionsHead(PetscOptionsObject,"ILU Options");
49: PCSetFromOptions_Factor(PetscOptionsObject,pc);
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: PetscOptionsTail();
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;
75: pc->failedreason = PC_NOERROR;
76: /* ugly hack to change default, since it is not support by some matrix types */
77: if (((PC_Factor*)ilu)->info.shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
78: PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJ,&flg);
79: if (!flg) {
80: PetscObjectTypeCompare((PetscObject)pc->pmat,MATMPIAIJ,&flg);
81: if (!flg) {
82: ((PC_Factor*)ilu)->info.shifttype = (PetscReal)MAT_SHIFT_INBLOCKS;
83: PetscInfo(pc,"Changing shift type from NONZERO to INBLOCKS because block matrices do not support NONZERO\n");
84: }
85: }
86: }
88: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
89: if (ilu->hdr.inplace) {
90: if (!pc->setupcalled) {
92: /* In-place factorization only makes sense with the natural ordering,
93: so we only need to get the ordering once, even if nonzero structure changes */
94: /* Should not get the ordering if the factorization routine does not use it, but do not yet have access to the factor matrix */
95: PCFactorSetDefaultOrdering_Factor(pc);
96: MatDestroy(&((PC_Factor*)ilu)->fact);
97: MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);
98: if (ilu->row) PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->row);
99: if (ilu->col) PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->col);
100: }
102: /* In place ILU only makes sense with fill factor of 1.0 because
103: cannot have levels of fill */
104: ((PC_Factor*)ilu)->info.fill = 1.0;
105: ((PC_Factor*)ilu)->info.diagonal_fill = 0.0;
107: MatILUFactor(pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);
108: MatFactorGetError(pc->pmat,&err);
109: if (err) { /* Factor() fails */
110: pc->failedreason = (PCFailedReason)err;
111: return 0;
112: }
114: ((PC_Factor*)ilu)->fact = pc->pmat;
115: /* must update the pc record of the matrix state or the PC will attempt to run PCSetUp() yet again */
116: PetscObjectStateGet((PetscObject)pc->pmat,&pc->matstate);
117: } else {
118: if (!pc->setupcalled) {
119: /* first time in so compute reordering and symbolic factorization */
120: PetscBool canuseordering;
121: if (!((PC_Factor*)ilu)->fact) {
122: MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);
123: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)ilu)->fact);
124: }
125: MatFactorGetCanUseOrdering(((PC_Factor*)ilu)->fact,&canuseordering);
126: if (canuseordering) {
127: PCFactorSetDefaultOrdering_Factor(pc);
128: MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);
129: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->row);
130: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->col);
131: /* Remove zeros along diagonal? */
132: if (ilu->nonzerosalongdiagonal) {
133: MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);
134: }
135: }
136: MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);
137: MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);
138: ilu->hdr.actualfill = info.fill_ratio_needed;
139: } else if (pc->flag != SAME_NONZERO_PATTERN) {
140: if (!ilu->hdr.reuseordering) {
141: PetscBool canuseordering;
142: MatDestroy(&((PC_Factor*)ilu)->fact);
143: MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);
144: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)ilu)->fact);
145: MatFactorGetCanUseOrdering(((PC_Factor*)ilu)->fact,&canuseordering);
146: if (canuseordering) {
147: /* compute a new ordering for the ILU */
148: ISDestroy(&ilu->row);
149: ISDestroy(&ilu->col);
150: PCFactorSetDefaultOrdering_Factor(pc);
151: MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);
152: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->row);
153: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->col);
154: /* Remove zeros along diagonal? */
155: if (ilu->nonzerosalongdiagonal) {
156: MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);
157: }
158: }
159: }
160: MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);
161: MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);
162: ilu->hdr.actualfill = info.fill_ratio_needed;
163: }
164: MatFactorGetError(((PC_Factor*)ilu)->fact,&err);
165: if (err) { /* FactorSymbolic() fails */
166: pc->failedreason = (PCFailedReason)err;
167: return 0;
168: }
170: MatLUFactorNumeric(((PC_Factor*)ilu)->fact,pc->pmat,&((PC_Factor*)ilu)->info);
171: MatFactorGetError(((PC_Factor*)ilu)->fact,&err);
172: if (err) { /* FactorNumeric() fails */
173: pc->failedreason = (PCFailedReason)err;
174: }
175: }
177: PCFactorGetMatSolverType(pc,&stype);
178: if (!stype) {
179: MatSolverType solverpackage;
180: MatFactorGetSolverType(((PC_Factor*)ilu)->fact,&solverpackage);
181: PCFactorSetMatSolverType(pc,solverpackage);
182: }
183: return 0;
184: }
186: static PetscErrorCode PCDestroy_ILU(PC pc)
187: {
188: PC_ILU *ilu = (PC_ILU*)pc->data;
190: PCReset_ILU(pc);
191: PetscFree(((PC_Factor*)ilu)->solvertype);
192: PetscFree(((PC_Factor*)ilu)->ordering);
193: PetscFree(pc->data);
194: return 0;
195: }
197: static PetscErrorCode PCApply_ILU(PC pc,Vec x,Vec y)
198: {
199: PC_ILU *ilu = (PC_ILU*)pc->data;
201: MatSolve(((PC_Factor*)ilu)->fact,x,y);
202: return 0;
203: }
205: static PetscErrorCode PCMatApply_ILU(PC pc,Mat X,Mat Y)
206: {
207: PC_ILU *ilu = (PC_ILU*)pc->data;
209: MatMatSolve(((PC_Factor*)ilu)->fact,X,Y);
210: return 0;
211: }
213: static PetscErrorCode PCApplyTranspose_ILU(PC pc,Vec x,Vec y)
214: {
215: PC_ILU *ilu = (PC_ILU*)pc->data;
217: MatSolveTranspose(((PC_Factor*)ilu)->fact,x,y);
218: return 0;
219: }
221: static PetscErrorCode PCApplySymmetricLeft_ILU(PC pc,Vec x,Vec y)
222: {
223: PC_ILU *icc = (PC_ILU*)pc->data;
225: MatForwardSolve(((PC_Factor*)icc)->fact,x,y);
226: return 0;
227: }
229: static PetscErrorCode PCApplySymmetricRight_ILU(PC pc,Vec x,Vec y)
230: {
231: PC_ILU *icc = (PC_ILU*)pc->data;
233: MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);
234: return 0;
235: }
237: /*MC
238: PCILU - Incomplete factorization preconditioners.
240: Options Database Keys:
241: + -pc_factor_levels <k> - number of levels of fill for ILU(k)
242: . -pc_factor_in_place - only for ILU(0) with natural ordering, reuses the space of the matrix for
243: its factorization (overwrites original matrix)
244: . -pc_factor_diagonal_fill - fill in a zero diagonal even if levels of fill indicate it wouldn't be fill
245: . -pc_factor_reuse_ordering - reuse ordering of factorized matrix from previous factorization
246: . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
247: . -pc_factor_nonzeros_along_diagonal - reorder the matrix before factorization to remove zeros from the diagonal,
248: this decreases the chance of getting a zero pivot
249: . -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
250: - -pc_factor_pivot_in_blocks - for block ILU(k) factorization, i.e. with BAIJ matrices with block size larger
251: than 1 the diagonal blocks are factored with partial pivoting (this increases the
252: stability of the ILU factorization
254: Level: beginner
256: Notes:
257: Only implemented for some matrix formats. (for parallel see PCHYPRE for hypre's ILU)
259: For BAIJ matrices this implements a point block ILU
261: The "symmetric" application of this preconditioner is not actually symmetric since L is not transpose(U)
262: even when the matrix is not symmetric since the U stores the diagonals of the factorization.
264: If you are using MATSEQAIJCUSPARSE matrices (or MATMPIAIJCUSPARSE matrices with block Jacobi), factorization
265: is never done on the GPU).
267: References:
268: + * - T. Dupont, R. Kendall, and H. Rachford. An approximate factorization procedure for solving
269: self adjoint elliptic difference equations. SIAM J. Numer. Anal., 5, 1968.
270: . * - T.A. Oliphant. An implicit numerical method for solving two dimensional timedependent diffusion problems. Quart. Appl. Math., 19, 1961.
271: - * - TONY F. CHAN AND HENK A. VAN DER VORST, APPROXIMATE AND INCOMPLETE FACTORIZATIONS,
272: Chapter in Parallel Numerical
273: Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
274: Science and Engineering, Kluwer.
276: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
277: PCFactorSetZeroPivot(), PCFactorSetShiftSetType(), PCFactorSetAmount(),
278: PCFactorSetDropTolerance(),PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(),
279: PCFactorSetLevels(), PCFactorSetUseInPlace(), PCFactorSetAllowDiagonalFill(), PCFactorSetPivotInBlocks(),
280: PCFactorGetAllowDiagonalFill(), PCFactorGetUseInPlace()
282: M*/
284: PETSC_EXTERN PetscErrorCode PCCreate_ILU(PC pc)
285: {
286: PC_ILU *ilu;
288: PetscNewLog(pc,&ilu);
289: pc->data = (void*)ilu;
290: 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->setup = PCSetUp_ILU;
306: pc->ops->setfromoptions = PCSetFromOptions_ILU;
307: pc->ops->view = PCView_Factor;
308: pc->ops->applysymmetricleft = PCApplySymmetricLeft_ILU;
309: pc->ops->applysymmetricright = PCApplySymmetricRight_ILU;
310: pc->ops->applyrichardson = NULL;
311: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetDropTolerance_C",PCFactorSetDropTolerance_ILU);
312: PetscObjectComposeFunction((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C",PCFactorReorderForNonzeroDiagonal_ILU);
313: return 0;
314: }