Actual source code: ilu.c
petsc-3.4.5 2014-06-29
2: /*
3: Defines a ILU factorization preconditioner for any Mat implementation
4: */
5: #include <../src/ksp/pc/impls/factor/ilu/ilu.h> /*I "petscpc.h" I*/
7: /* ------------------------------------------------------------------------------------------*/
10: PetscErrorCode PCFactorSetReuseFill_ILU(PC pc,PetscBool flag)
11: {
12: PC_ILU *lu = (PC_ILU*)pc->data;
15: lu->reusefill = flag;
16: return(0);
17: }
21: PetscErrorCode PCFactorReorderForNonzeroDiagonal_ILU(PC pc,PetscReal z)
22: {
23: PC_ILU *ilu = (PC_ILU*)pc->data;
26: ilu->nonzerosalongdiagonal = PETSC_TRUE;
27: if (z == PETSC_DECIDE) ilu->nonzerosalongdiagonaltol = 1.e-10;
28: else ilu->nonzerosalongdiagonaltol = z;
29: return(0);
30: }
34: PetscErrorCode PCReset_ILU(PC pc)
35: {
36: PC_ILU *ilu = (PC_ILU*)pc->data;
40: if (!ilu->inplace) {MatDestroy(&((PC_Factor*)ilu)->fact);}
41: if (ilu->row && ilu->col && ilu->row != ilu->col) {ISDestroy(&ilu->row);}
42: ISDestroy(&ilu->col);
43: return(0);
44: }
48: PetscErrorCode PCFactorSetDropTolerance_ILU(PC pc,PetscReal dt,PetscReal dtcol,PetscInt dtcount)
49: {
50: PC_ILU *ilu = (PC_ILU*)pc->data;
53: if (pc->setupcalled && (((PC_Factor*)ilu)->info.dt != dt || ((PC_Factor*)ilu)->info.dtcol != dtcol || ((PC_Factor*)ilu)->info.dtcount != dtcount)) {
54: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot change drop tolerance after using PC");
55: }
56: ((PC_Factor*)ilu)->info.dt = dt;
57: ((PC_Factor*)ilu)->info.dtcol = dtcol;
58: ((PC_Factor*)ilu)->info.dtcount = dtcount;
59: ((PC_Factor*)ilu)->info.usedt = 1.0;
60: return(0);
61: }
65: PetscErrorCode PCFactorSetReuseOrdering_ILU(PC pc,PetscBool flag)
66: {
67: PC_ILU *ilu = (PC_ILU*)pc->data;
70: ilu->reuseordering = flag;
71: return(0);
72: }
76: PetscErrorCode PCFactorSetUseInPlace_ILU(PC pc)
77: {
78: PC_ILU *dir = (PC_ILU*)pc->data;
81: dir->inplace = PETSC_TRUE;
82: return(0);
83: }
87: static PetscErrorCode PCSetFromOptions_ILU(PC pc)
88: {
90: PetscInt itmp;
91: PetscBool flg;
92: PC_ILU *ilu = (PC_ILU*)pc->data;
93: PetscReal tol;
94: /* PetscReal dt[3]; */
97: PetscOptionsHead("ILU Options");
98: PCSetFromOptions_Factor(pc);
100: PetscOptionsInt("-pc_factor_levels","levels of fill","PCFactorSetLevels",(PetscInt)((PC_Factor*)ilu)->info.levels,&itmp,&flg);
101: if (flg) ((PC_Factor*)ilu)->info.levels = itmp;
103: flg = PETSC_FALSE;
104: PetscOptionsBool("-pc_factor_diagonal_fill","Allow fill into empty diagonal entry","PCFactorSetAllowDiagonalFill",flg,&flg,NULL);
105: ((PC_Factor*)ilu)->info.diagonal_fill = (double) flg;
106: /*
107: dt[0] = ((PC_Factor*)ilu)->info.dt;
108: dt[1] = ((PC_Factor*)ilu)->info.dtcol;
109: dt[2] = ((PC_Factor*)ilu)->info.dtcount;
111: PetscInt dtmax = 3;
112: PetscOptionsRealArray("-pc_factor_drop_tolerance,","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg);
113: if (flg) {
114: PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);
115: }
116: */
117: PetscOptionsName("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",&flg);
118: if (flg) {
119: tol = PETSC_DECIDE;
120: PetscOptionsReal("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",ilu->nonzerosalongdiagonaltol,&tol,0);
121: PCFactorReorderForNonzeroDiagonal(pc,tol);
122: }
124: PetscOptionsTail();
125: return(0);
126: }
130: static PetscErrorCode PCView_ILU(PC pc,PetscViewer viewer)
131: {
132: PC_ILU *ilu = (PC_ILU*)pc->data;
134: PetscBool iascii;
137: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
138: if (iascii) {
139: if (ilu->inplace) {
140: PetscViewerASCIIPrintf(viewer," ILU: in-place factorization\n");
141: } else {
142: PetscViewerASCIIPrintf(viewer," ILU: out-of-place factorization\n");
143: }
145: if (ilu->reusefill) {PetscViewerASCIIPrintf(viewer," ILU: Reusing fill from past factorization\n");}
146: if (ilu->reuseordering) {PetscViewerASCIIPrintf(viewer," ILU: Reusing reordering from past factorization\n");}
147: }
148: PCView_Factor(pc,viewer);
149: return(0);
150: }
154: static PetscErrorCode PCSetUp_ILU(PC pc)
155: {
157: PC_ILU *ilu = (PC_ILU*)pc->data;
158: MatInfo info;
159: PetscBool flg;
162: /* ugly hack to change default, since it is not support by some matrix types */
163: if (((PC_Factor*)ilu)->info.shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
164: PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJ,&flg);
165: if (!flg) {
166: PetscObjectTypeCompare((PetscObject)pc->pmat,MATMPIAIJ,&flg);
167: if (!flg) {
168: ((PC_Factor*)ilu)->info.shifttype = (PetscReal)MAT_SHIFT_INBLOCKS;
169: PetscInfo(pc,"Changing shift type from NONZERO to INBLOCKS because block matrices do not support NONZERO");
170: }
171: }
172: }
174: if (ilu->inplace) {
175: if (!pc->setupcalled) {
177: /* In-place factorization only makes sense with the natural ordering,
178: so we only need to get the ordering once, even if nonzero structure changes */
179: MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);
180: if (ilu->row) {PetscLogObjectParent(pc,ilu->row);}
181: if (ilu->col) {PetscLogObjectParent(pc,ilu->col);}
182: }
184: /* In place ILU only makes sense with fill factor of 1.0 because
185: cannot have levels of fill */
186: ((PC_Factor*)ilu)->info.fill = 1.0;
187: ((PC_Factor*)ilu)->info.diagonal_fill = 0.0;
189: MatILUFactor(pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);
190: ((PC_Factor*)ilu)->fact = pc->pmat;
191: } else {
192: if (!pc->setupcalled) {
193: /* first time in so compute reordering and symbolic factorization */
194: MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);
195: if (ilu->row) {PetscLogObjectParent(pc,ilu->row);}
196: if (ilu->col) {PetscLogObjectParent(pc,ilu->col);}
197: /* Remove zeros along diagonal? */
198: if (ilu->nonzerosalongdiagonal) {
199: MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);
200: }
201: if (!((PC_Factor*)ilu)->fact) {
202: MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);
203: }
204: MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);
205: MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);
207: ilu->actualfill = info.fill_ratio_needed;
209: PetscLogObjectParent(pc,((PC_Factor*)ilu)->fact);
210: } else if (pc->flag != SAME_NONZERO_PATTERN) {
211: if (!ilu->reuseordering) {
212: /* compute a new ordering for the ILU */
213: ISDestroy(&ilu->row);
214: ISDestroy(&ilu->col);
215: MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);
216: if (ilu->row) {PetscLogObjectParent(pc,ilu->row);}
217: if (ilu->col) {PetscLogObjectParent(pc,ilu->col);}
218: /* Remove zeros along diagonal? */
219: if (ilu->nonzerosalongdiagonal) {
220: MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);
221: }
222: }
223: MatDestroy(&((PC_Factor*)ilu)->fact);
224: MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);
225: MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);
226: MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);
228: ilu->actualfill = info.fill_ratio_needed;
230: PetscLogObjectParent(pc,((PC_Factor*)ilu)->fact);
231: }
232: MatLUFactorNumeric(((PC_Factor*)ilu)->fact,pc->pmat,&((PC_Factor*)ilu)->info);
233: }
234: return(0);
235: }
239: static PetscErrorCode PCDestroy_ILU(PC pc)
240: {
241: PC_ILU *ilu = (PC_ILU*)pc->data;
245: PCReset_ILU(pc);
246: PetscFree(((PC_Factor*)ilu)->solvertype);
247: PetscFree(((PC_Factor*)ilu)->ordering);
248: PetscFree(pc->data);
249: return(0);
250: }
254: static PetscErrorCode PCApply_ILU(PC pc,Vec x,Vec y)
255: {
256: PC_ILU *ilu = (PC_ILU*)pc->data;
260: MatSolve(((PC_Factor*)ilu)->fact,x,y);
261: return(0);
262: }
266: static PetscErrorCode PCApplyTranspose_ILU(PC pc,Vec x,Vec y)
267: {
268: PC_ILU *ilu = (PC_ILU*)pc->data;
272: MatSolveTranspose(((PC_Factor*)ilu)->fact,x,y);
273: return(0);
274: }
278: static PetscErrorCode PCApplySymmetricLeft_ILU(PC pc,Vec x,Vec y)
279: {
281: PC_ILU *icc = (PC_ILU*)pc->data;
284: MatForwardSolve(((PC_Factor*)icc)->fact,x,y);
285: return(0);
286: }
290: static PetscErrorCode PCApplySymmetricRight_ILU(PC pc,Vec x,Vec y)
291: {
293: PC_ILU *icc = (PC_ILU*)pc->data;
296: MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);
297: return(0);
298: }
300: /*MC
301: PCILU - Incomplete factorization preconditioners.
303: Options Database Keys:
304: + -pc_factor_levels <k> - number of levels of fill for ILU(k)
305: . -pc_factor_in_place - only for ILU(0) with natural ordering, reuses the space of the matrix for
306: its factorization (overwrites original matrix)
307: . -pc_factor_diagonal_fill - fill in a zero diagonal even if levels of fill indicate it wouldn't be fill
308: . -pc_factor_reuse_ordering - reuse ordering of factorized matrix from previous factorization
309: . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
310: . -pc_factor_nonzeros_along_diagonal - reorder the matrix before factorization to remove zeros from the diagonal,
311: this decreases the chance of getting a zero pivot
312: . -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
313: . -pc_factor_pivot_in_blocks - for block ILU(k) factorization, i.e. with BAIJ matrices with block size larger
314: than 1 the diagonal blocks are factored with partial pivoting (this increases the
315: stability of the ILU factorization
317: Level: beginner
319: Concepts: incomplete factorization
321: Notes: Only implemented for some matrix formats. (for parallel see PCHYPRE for hypre's ILU)
323: For BAIJ matrices this implements a point block ILU
325: The "symmetric" application of this preconditioner is not actually symmetric since L is not transpose(U)
326: even when the matrix is not symmetric since the U stores the diagonals of the factorization.
328: If you are using MATSEQAIJCUSP matrices (or MATMPIAIJCUSP matrices with block Jacobi) you must have ./configured
329: PETSc with also --download-txpetscgpu to have the triangular solves performed on the GPU (factorization is never
330: done on the GPU).
332: References:
333: T. Dupont, R. Kendall, and H. Rachford. An approximate factorization procedure for solving
334: self-adjoint elliptic difference equations. SIAM J. Numer. Anal., 5:559--573, 1968.
336: T.A. Oliphant. An implicit numerical method for solving two-dimensional time-dependent dif-
337: fusion problems. Quart. Appl. Math., 19:221--229, 1961.
339: Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS, TONY F. CHAN AND HENK A. VAN DER VORST
340: http://igitur-archive.library.uu.nl/math/2001-0621-115821/proc.pdf chapter in Parallel Numerical
341: Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
342: Science and Engineering, Kluwer, pp. 167--202.
345: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
346: PCFactorSetZeroPivot(), PCFactorSetShiftSetType(), PCFactorSetAmount(),
347: PCFactorSetDropTolerance(),PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(),
348: PCFactorSetLevels(), PCFactorSetUseInPlace(), PCFactorSetAllowDiagonalFill(), PCFactorSetPivotInBlocks()
350: M*/
354: PETSC_EXTERN PetscErrorCode PCCreate_ILU(PC pc)
355: {
357: PC_ILU *ilu;
360: PetscNewLog(pc,PC_ILU,&ilu);
362: ((PC_Factor*)ilu)->fact = 0;
363: MatFactorInfoInitialize(&((PC_Factor*)ilu)->info);
364: ((PC_Factor*)ilu)->factortype = MAT_FACTOR_ILU;
365: ((PC_Factor*)ilu)->info.levels = 0.;
366: ((PC_Factor*)ilu)->info.fill = 1.0;
367: ilu->col = 0;
368: ilu->row = 0;
369: ilu->inplace = PETSC_FALSE;
370: PetscStrallocpy(MATSOLVERPETSC,&((PC_Factor*)ilu)->solvertype);
371: PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)ilu)->ordering);
372: ilu->reuseordering = PETSC_FALSE;
373: ((PC_Factor*)ilu)->info.dt = PETSC_DEFAULT;
374: ((PC_Factor*)ilu)->info.dtcount = PETSC_DEFAULT;
375: ((PC_Factor*)ilu)->info.dtcol = PETSC_DEFAULT;
376: ((PC_Factor*)ilu)->info.shifttype = (PetscReal)MAT_SHIFT_INBLOCKS;
377: ((PC_Factor*)ilu)->info.shiftamount = 100.0*PETSC_MACHINE_EPSILON;
378: ((PC_Factor*)ilu)->info.zeropivot = 100.0*PETSC_MACHINE_EPSILON;
379: ((PC_Factor*)ilu)->info.pivotinblocks = 1.0;
380: ilu->reusefill = PETSC_FALSE;
381: ((PC_Factor*)ilu)->info.diagonal_fill = 0.0;
382: pc->data = (void*)ilu;
384: pc->ops->reset = PCReset_ILU;
385: pc->ops->destroy = PCDestroy_ILU;
386: pc->ops->apply = PCApply_ILU;
387: pc->ops->applytranspose = PCApplyTranspose_ILU;
388: pc->ops->setup = PCSetUp_ILU;
389: pc->ops->setfromoptions = PCSetFromOptions_ILU;
390: pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor;
391: pc->ops->view = PCView_ILU;
392: pc->ops->applysymmetricleft = PCApplySymmetricLeft_ILU;
393: pc->ops->applysymmetricright = PCApplySymmetricRight_ILU;
394: pc->ops->applyrichardson = 0;
396: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetZeroPivot_C",PCFactorSetZeroPivot_Factor);
397: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Factor);
398: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftAmount_C",PCFactorSetShiftAmount_Factor);
399: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetMatSolverPackage_C",PCFactorGetMatSolverPackage_Factor);
400: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatSolverPackage_C",PCFactorSetMatSolverPackage_Factor);
401: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C",PCFactorSetUpMatSolverPackage_Factor);
402: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetDropTolerance_C",PCFactorSetDropTolerance_ILU);
403: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetFill_C",PCFactorSetFill_Factor);
404: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatOrderingType_C",PCFactorSetMatOrderingType_Factor);
405: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C",PCFactorSetReuseOrdering_ILU);
406: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseFill_C",PCFactorSetReuseFill_ILU);
407: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetLevels_C",PCFactorSetLevels_Factor);
408: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetLevels_C",PCFactorGetLevels_Factor);
409: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUseInPlace_C",PCFactorSetUseInPlace_ILU);
410: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetAllowDiagonalFill_C",PCFactorSetAllowDiagonalFill_Factor);
411: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetPivotInBlocks_C",PCFactorSetPivotInBlocks_Factor);
412: PetscObjectComposeFunction((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C",PCFactorReorderForNonzeroDiagonal_ILU);
413: return(0);
414: }