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;
12: ilu->nonzerosalongdiagonal = PETSC_TRUE;
13: if (z == PETSC_DECIDE) ilu->nonzerosalongdiagonaltol = 1.e-10;
14: else ilu->nonzerosalongdiagonaltol = z;
15: return(0);
16: }
18: PetscErrorCode PCReset_ILU(PC pc)
19: {
20: PC_ILU *ilu = (PC_ILU*)pc->data;
24: if (!ilu->hdr.inplace) {MatDestroy(&((PC_Factor*)ilu)->fact);}
25: if (ilu->row && ilu->col && ilu->row != ilu->col) {ISDestroy(&ilu->row);}
26: ISDestroy(&ilu->col);
27: return(0);
28: }
30: PetscErrorCode PCFactorSetDropTolerance_ILU(PC pc,PetscReal dt,PetscReal dtcol,PetscInt dtcount)
31: {
32: PC_ILU *ilu = (PC_ILU*)pc->data;
35: if (pc->setupcalled && (((PC_Factor*)ilu)->info.dt != dt || ((PC_Factor*)ilu)->info.dtcol != dtcol || ((PC_Factor*)ilu)->info.dtcount != dtcount)) {
36: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot change drop tolerance after using PC");
37: }
38: ((PC_Factor*)ilu)->info.dt = dt;
39: ((PC_Factor*)ilu)->info.dtcol = dtcol;
40: ((PC_Factor*)ilu)->info.dtcount = dtcount;
41: ((PC_Factor*)ilu)->info.usedt = 1.0;
42: return(0);
43: }
45: static PetscErrorCode PCSetFromOptions_ILU(PetscOptionItems *PetscOptionsObject,PC pc)
46: {
48: PetscInt itmp;
49: PetscBool flg,set;
50: PC_ILU *ilu = (PC_ILU*)pc->data;
51: PetscReal tol;
54: PetscOptionsHead(PetscOptionsObject,"ILU Options");
55: PCSetFromOptions_Factor(PetscOptionsObject,pc);
57: PetscOptionsInt("-pc_factor_levels","levels of fill","PCFactorSetLevels",(PetscInt)((PC_Factor*)ilu)->info.levels,&itmp,&flg);
58: if (flg) ((PC_Factor*)ilu)->info.levels = itmp;
60: PetscOptionsBool("-pc_factor_diagonal_fill","Allow fill into empty diagonal entry","PCFactorSetAllowDiagonalFill",((PC_Factor*)ilu)->info.diagonal_fill ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
61: if (set) ((PC_Factor*)ilu)->info.diagonal_fill = (PetscReal) flg;
62: PetscOptionsName("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",&flg);
63: if (flg) {
64: tol = PETSC_DECIDE;
65: PetscOptionsReal("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",ilu->nonzerosalongdiagonaltol,&tol,NULL);
66: PCFactorReorderForNonzeroDiagonal(pc,tol);
67: }
69: PetscOptionsTail();
70: return(0);
71: }
73: static PetscErrorCode PCSetUp_ILU(PC pc)
74: {
75: PetscErrorCode ierr;
76: PC_ILU *ilu = (PC_ILU*)pc->data;
77: MatInfo info;
78: PetscBool flg;
79: MatSolverType stype;
80: MatFactorError err;
83: pc->failedreason = PC_NOERROR;
84: /* ugly hack to change default, since it is not support by some matrix types */
85: if (((PC_Factor*)ilu)->info.shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
86: PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJ,&flg);
87: if (!flg) {
88: PetscObjectTypeCompare((PetscObject)pc->pmat,MATMPIAIJ,&flg);
89: if (!flg) {
90: ((PC_Factor*)ilu)->info.shifttype = (PetscReal)MAT_SHIFT_INBLOCKS;
91: PetscInfo(pc,"Changing shift type from NONZERO to INBLOCKS because block matrices do not support NONZERO\n");
92: }
93: }
94: }
96: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
97: if (ilu->hdr.inplace) {
98: 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: PCFactorSetDefaultOrdering_Factor(pc);
104: MatDestroy(&((PC_Factor*)ilu)->fact);
105: MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);
106: if (ilu->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->row);}
107: if (ilu->col) {PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->col);}
108: }
110: /* In place ILU only makes sense with fill factor of 1.0 because
111: cannot have levels of fill */
112: ((PC_Factor*)ilu)->info.fill = 1.0;
113: ((PC_Factor*)ilu)->info.diagonal_fill = 0.0;
115: MatILUFactor(pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);
116: MatFactorGetError(pc->pmat,&err);
117: if (err) { /* Factor() fails */
118: pc->failedreason = (PCFailedReason)err;
119: return(0);
120: }
122: ((PC_Factor*)ilu)->fact = pc->pmat;
123: /* must update the pc record of the matrix state or the PC will attempt to run PCSetUp() yet again */
124: PetscObjectStateGet((PetscObject)pc->pmat,&pc->matstate);
125: } else {
126: if (!pc->setupcalled) {
127: /* first time in so compute reordering and symbolic factorization */
128: PetscBool canuseordering;
129: if (!((PC_Factor*)ilu)->fact) {
130: MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);
131: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)ilu)->fact);
132: }
133: MatFactorGetCanUseOrdering(((PC_Factor*)ilu)->fact,&canuseordering);
134: if (canuseordering) {
135: PCFactorSetDefaultOrdering_Factor(pc);
136: MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);
137: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->row);
138: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->col);
139: /* Remove zeros along diagonal? */
140: if (ilu->nonzerosalongdiagonal) {
141: MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);
142: }
143: }
144: MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);
145: MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);
146: ilu->hdr.actualfill = info.fill_ratio_needed;
147: } else if (pc->flag != SAME_NONZERO_PATTERN) {
148: if (!ilu->hdr.reuseordering) {
149: PetscBool canuseordering;
150: MatDestroy(&((PC_Factor*)ilu)->fact);
151: MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);
152: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)ilu)->fact);
153: MatFactorGetCanUseOrdering(((PC_Factor*)ilu)->fact,&canuseordering);
154: if (canuseordering) {
155: /* compute a new ordering for the ILU */
156: ISDestroy(&ilu->row);
157: ISDestroy(&ilu->col);
158: PCFactorSetDefaultOrdering_Factor(pc);
159: MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);
160: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->row);
161: PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->col);
162: /* Remove zeros along diagonal? */
163: if (ilu->nonzerosalongdiagonal) {
164: MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);
165: }
166: }
167: }
168: MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);
169: MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);
170: ilu->hdr.actualfill = info.fill_ratio_needed;
171: }
172: MatFactorGetError(((PC_Factor*)ilu)->fact,&err);
173: if (err) { /* FactorSymbolic() fails */
174: pc->failedreason = (PCFailedReason)err;
175: return(0);
176: }
178: MatLUFactorNumeric(((PC_Factor*)ilu)->fact,pc->pmat,&((PC_Factor*)ilu)->info);
179: MatFactorGetError(((PC_Factor*)ilu)->fact,&err);
180: if (err) { /* FactorNumeric() fails */
181: pc->failedreason = (PCFailedReason)err;
182: }
183: }
185: PCFactorGetMatSolverType(pc,&stype);
186: if (!stype) {
187: MatSolverType solverpackage;
188: MatFactorGetSolverType(((PC_Factor*)ilu)->fact,&solverpackage);
189: PCFactorSetMatSolverType(pc,solverpackage);
190: }
191: return(0);
192: }
194: static PetscErrorCode PCDestroy_ILU(PC pc)
195: {
196: PC_ILU *ilu = (PC_ILU*)pc->data;
200: PCReset_ILU(pc);
201: PetscFree(((PC_Factor*)ilu)->solvertype);
202: PetscFree(((PC_Factor*)ilu)->ordering);
203: PetscFree(pc->data);
204: return(0);
205: }
207: static PetscErrorCode PCApply_ILU(PC pc,Vec x,Vec y)
208: {
209: PC_ILU *ilu = (PC_ILU*)pc->data;
213: MatSolve(((PC_Factor*)ilu)->fact,x,y);
214: return(0);
215: }
217: static PetscErrorCode PCMatApply_ILU(PC pc,Mat X,Mat Y)
218: {
219: PC_ILU *ilu = (PC_ILU*)pc->data;
223: MatMatSolve(((PC_Factor*)ilu)->fact,X,Y);
224: return(0);
225: }
227: static PetscErrorCode PCApplyTranspose_ILU(PC pc,Vec x,Vec y)
228: {
229: PC_ILU *ilu = (PC_ILU*)pc->data;
233: MatSolveTranspose(((PC_Factor*)ilu)->fact,x,y);
234: return(0);
235: }
237: static PetscErrorCode PCApplySymmetricLeft_ILU(PC pc,Vec x,Vec y)
238: {
240: PC_ILU *icc = (PC_ILU*)pc->data;
243: MatForwardSolve(((PC_Factor*)icc)->fact,x,y);
244: return(0);
245: }
247: static PetscErrorCode PCApplySymmetricRight_ILU(PC pc,Vec x,Vec y)
248: {
250: PC_ILU *icc = (PC_ILU*)pc->data;
253: MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);
254: return(0);
255: }
257: /*MC
258: PCILU - Incomplete factorization preconditioners.
260: Options Database Keys:
261: + -pc_factor_levels <k> - number of levels of fill for ILU(k)
262: . -pc_factor_in_place - only for ILU(0) with natural ordering, reuses the space of the matrix for
263: its factorization (overwrites original matrix)
264: . -pc_factor_diagonal_fill - fill in a zero diagonal even if levels of fill indicate it wouldn't be fill
265: . -pc_factor_reuse_ordering - reuse ordering of factorized matrix from previous factorization
266: . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
267: . -pc_factor_nonzeros_along_diagonal - reorder the matrix before factorization to remove zeros from the diagonal,
268: this decreases the chance of getting a zero pivot
269: . -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
270: - -pc_factor_pivot_in_blocks - for block ILU(k) factorization, i.e. with BAIJ matrices with block size larger
271: than 1 the diagonal blocks are factored with partial pivoting (this increases the
272: stability of the ILU factorization
274: Level: beginner
276: Notes:
277: Only implemented for some matrix formats. (for parallel see PCHYPRE for hypre's ILU)
279: For BAIJ matrices this implements a point block ILU
281: The "symmetric" application of this preconditioner is not actually symmetric since L is not transpose(U)
282: even when the matrix is not symmetric since the U stores the diagonals of the factorization.
284: If you are using MATSEQAIJCUSPARSE matrices (or MATMPIAIJCUSPARSE matrices with block Jacobi), factorization
285: is never done on the GPU).
287: References:
288: + 1. - T. Dupont, R. Kendall, and H. Rachford. An approximate factorization procedure for solving
289: self adjoint elliptic difference equations. SIAM J. Numer. Anal., 5, 1968.
290: . 2. - T.A. Oliphant. An implicit numerical method for solving two dimensional timedependent diffusion problems. Quart. Appl. Math., 19, 1961.
291: - 3. - TONY F. CHAN AND HENK A. VAN DER VORST, APPROXIMATE AND INCOMPLETE FACTORIZATIONS,
292: Chapter in Parallel Numerical
293: Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
294: Science and Engineering, Kluwer.
296: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
297: PCFactorSetZeroPivot(), PCFactorSetShiftSetType(), PCFactorSetAmount(),
298: PCFactorSetDropTolerance(),PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(),
299: PCFactorSetLevels(), PCFactorSetUseInPlace(), PCFactorSetAllowDiagonalFill(), PCFactorSetPivotInBlocks(),
300: PCFactorGetAllowDiagonalFill(), PCFactorGetUseInPlace()
302: M*/
304: PETSC_EXTERN PetscErrorCode PCCreate_ILU(PC pc)
305: {
307: PC_ILU *ilu;
310: PetscNewLog(pc,&ilu);
311: pc->data = (void*)ilu;
312: PCFactorInitialize(pc,MAT_FACTOR_ILU);
314: ((PC_Factor*)ilu)->info.levels = 0.;
315: ((PC_Factor*)ilu)->info.fill = 1.0;
316: ilu->col = NULL;
317: ilu->row = NULL;
318: ((PC_Factor*)ilu)->info.dt = PETSC_DEFAULT;
319: ((PC_Factor*)ilu)->info.dtcount = PETSC_DEFAULT;
320: ((PC_Factor*)ilu)->info.dtcol = PETSC_DEFAULT;
322: pc->ops->reset = PCReset_ILU;
323: pc->ops->destroy = PCDestroy_ILU;
324: pc->ops->apply = PCApply_ILU;
325: pc->ops->matapply = PCMatApply_ILU;
326: pc->ops->applytranspose = PCApplyTranspose_ILU;
327: pc->ops->setup = PCSetUp_ILU;
328: pc->ops->setfromoptions = PCSetFromOptions_ILU;
329: pc->ops->view = PCView_Factor;
330: pc->ops->applysymmetricleft = PCApplySymmetricLeft_ILU;
331: pc->ops->applysymmetricright = PCApplySymmetricRight_ILU;
332: pc->ops->applyrichardson = NULL;
333: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetDropTolerance_C",PCFactorSetDropTolerance_ILU);
334: PetscObjectComposeFunction((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C",PCFactorReorderForNonzeroDiagonal_ILU);
335: return(0);
336: }