Actual source code: ilu.c
petsc-3.3-p7 2013-05-11
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: /* ------------------------------------------------------------------------------------------*/
8: EXTERN_C_BEGIN
11: PetscErrorCode PCFactorSetReuseFill_ILU(PC pc,PetscBool flag)
12: {
13: PC_ILU *lu = (PC_ILU*)pc->data;
14:
16: lu->reusefill = flag;
17: return(0);
18: }
19: EXTERN_C_END
21: EXTERN_C_BEGIN
24: PetscErrorCode PCFactorReorderForNonzeroDiagonal_ILU(PC pc,PetscReal z)
25: {
26: PC_ILU *ilu = (PC_ILU*)pc->data;
29: ilu->nonzerosalongdiagonal = PETSC_TRUE;
30: if (z == PETSC_DECIDE) {
31: ilu->nonzerosalongdiagonaltol = 1.e-10;
32: } else {
33: ilu->nonzerosalongdiagonaltol = z;
34: }
35: return(0);
36: }
37: EXTERN_C_END
41: PetscErrorCode PCReset_ILU(PC pc)
42: {
43: PC_ILU *ilu = (PC_ILU*)pc->data;
47: if (!ilu->inplace) {MatDestroy(&((PC_Factor*)ilu)->fact);}
48: if (ilu->row && ilu->col && ilu->row != ilu->col) {ISDestroy(&ilu->row);}
49: ISDestroy(&ilu->col);
50: return(0);
51: }
53: EXTERN_C_BEGIN
56: PetscErrorCode PCFactorSetDropTolerance_ILU(PC pc,PetscReal dt,PetscReal dtcol,PetscInt dtcount)
57: {
58: PC_ILU *ilu = (PC_ILU*)pc->data;
61: if (pc->setupcalled && (((PC_Factor*)ilu)->info.dt != dt || ((PC_Factor*)ilu)->info.dtcol != dtcol || ((PC_Factor*)ilu)->info.dtcount != dtcount)) {
62: SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Cannot change drop tolerance after using PC");
63: }
64: ((PC_Factor*)ilu)->info.dt = dt;
65: ((PC_Factor*)ilu)->info.dtcol = dtcol;
66: ((PC_Factor*)ilu)->info.dtcount = dtcount;
67: ((PC_Factor*)ilu)->info.usedt = 1.0;
68: return(0);
69: }
70: EXTERN_C_END
72: EXTERN_C_BEGIN
75: PetscErrorCode PCFactorSetReuseOrdering_ILU(PC pc,PetscBool flag)
76: {
77: PC_ILU *ilu = (PC_ILU*)pc->data;
80: ilu->reuseordering = flag;
81: return(0);
82: }
83: EXTERN_C_END
85: EXTERN_C_BEGIN
88: PetscErrorCode PCFactorSetUseInPlace_ILU(PC pc)
89: {
90: PC_ILU *dir = (PC_ILU*)pc->data;
93: dir->inplace = PETSC_TRUE;
94: return(0);
95: }
96: EXTERN_C_END
100: static PetscErrorCode PCSetFromOptions_ILU(PC pc)
101: {
103: PetscInt itmp;
104: PetscBool flg;
105: /* PetscReal dt[3]; */
106: PC_ILU *ilu = (PC_ILU*)pc->data;
107: PetscReal tol;
110: PetscOptionsHead("ILU Options");
111: PCSetFromOptions_Factor(pc);
113: PetscOptionsInt("-pc_factor_levels","levels of fill","PCFactorSetLevels",(PetscInt)((PC_Factor*)ilu)->info.levels,&itmp,&flg);
114: if (flg) ((PC_Factor*)ilu)->info.levels = itmp;
115: flg = PETSC_FALSE;
116: PetscOptionsBool("-pc_factor_diagonal_fill","Allow fill into empty diagonal entry","PCFactorSetAllowDiagonalFill",flg,&flg,PETSC_NULL);
117: ((PC_Factor*)ilu)->info.diagonal_fill = (double) flg;
118: /*
119: dt[0] = ((PC_Factor*)ilu)->info.dt;
120: dt[1] = ((PC_Factor*)ilu)->info.dtcol;
121: dt[2] = ((PC_Factor*)ilu)->info.dtcount;
122:
123: PetscInt dtmax = 3;
124: PetscOptionsRealArray("-pc_factor_drop_tolerance,","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg);
125: if (flg) {
126: PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);
127: }
128: */
129: PetscOptionsName("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",&flg);
130: if (flg) {
131: tol = PETSC_DECIDE;
132: PetscOptionsReal("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",ilu->nonzerosalongdiagonaltol,&tol,0);
133: PCFactorReorderForNonzeroDiagonal(pc,tol);
134: }
136: PetscOptionsTail();
137: return(0);
138: }
142: static PetscErrorCode PCView_ILU(PC pc,PetscViewer viewer)
143: {
144: PC_ILU *ilu = (PC_ILU*)pc->data;
146: PetscBool iascii;
149: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
150: if (iascii) {
151: if (ilu->inplace) {
152: PetscViewerASCIIPrintf(viewer," ILU: in-place factorization\n");
153: } else {
154: PetscViewerASCIIPrintf(viewer," ILU: out-of-place factorization\n");
155: }
156:
157: if (ilu->reusefill) {PetscViewerASCIIPrintf(viewer," ILU: Reusing fill from past factorization\n");}
158: if (ilu->reuseordering) {PetscViewerASCIIPrintf(viewer," ILU: Reusing reordering from past factorization\n");}
159: }
160: PCView_Factor(pc,viewer);
161: return(0);
162: }
166: static PetscErrorCode PCSetUp_ILU(PC pc)
167: {
169: PC_ILU *ilu = (PC_ILU*)pc->data;
170: MatInfo info;
171: PetscBool flg;
174: /* ugly hack to change default, since it is not support by some matrix types */
175: if (((PC_Factor*)ilu)->info.shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
176: PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJ,&flg);
177: if (!flg) {
178: PetscObjectTypeCompare((PetscObject)pc->pmat,MATMPIAIJ,&flg);
179: if (!flg) {
180: ((PC_Factor*)ilu)->info.shifttype = (PetscReal)MAT_SHIFT_INBLOCKS;
181: PetscInfo(pc,"Changing shift type from NONZERO to INBLOCKS because block matrices do not support NONZERO");
182: }
183: }
184: }
186: if (ilu->inplace) {
187: CHKMEMQ;
188: if (!pc->setupcalled) {
190: /* In-place factorization only makes sense with the natural ordering,
191: so we only need to get the ordering once, even if nonzero structure changes */
192: MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);
193: if (ilu->row) {PetscLogObjectParent(pc,ilu->row);}
194: if (ilu->col) {PetscLogObjectParent(pc,ilu->col);}
195: }
197: /* In place ILU only makes sense with fill factor of 1.0 because
198: cannot have levels of fill */
199: ((PC_Factor*)ilu)->info.fill = 1.0;
200: ((PC_Factor*)ilu)->info.diagonal_fill = 0.0;
201: MatILUFactor(pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);CHKMEMQ;
202: ((PC_Factor*)ilu)->fact = pc->pmat;
203: } else {
204: if (!pc->setupcalled) {
205: /* first time in so compute reordering and symbolic factorization */
206: MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);
207: if (ilu->row) {PetscLogObjectParent(pc,ilu->row);}
208: if (ilu->col) {PetscLogObjectParent(pc,ilu->col);}
209: /* Remove zeros along diagonal? */
210: if (ilu->nonzerosalongdiagonal) {
211: MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);
212: }
213: if (!((PC_Factor*)ilu)->fact){
214: MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);
215: }
216: MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);
217: MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);
218: ilu->actualfill = info.fill_ratio_needed;
219: PetscLogObjectParent(pc,((PC_Factor*)ilu)->fact);
220: } else if (pc->flag != SAME_NONZERO_PATTERN) {
221: if (!ilu->reuseordering) {
222: /* compute a new ordering for the ILU */
223: ISDestroy(&ilu->row);
224: ISDestroy(&ilu->col);
225: MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);
226: if (ilu->row) {PetscLogObjectParent(pc,ilu->row);}
227: if (ilu->col) {PetscLogObjectParent(pc,ilu->col);}
228: /* Remove zeros along diagonal? */
229: if (ilu->nonzerosalongdiagonal) {
230: MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);
231: }
232: }
233: MatDestroy(&((PC_Factor*)ilu)->fact);
234: MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);
235: MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);
236: MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);
237: ilu->actualfill = info.fill_ratio_needed;
238: PetscLogObjectParent(pc,((PC_Factor*)ilu)->fact);
239: }
240: CHKMEMQ;
241: MatLUFactorNumeric(((PC_Factor*)ilu)->fact,pc->pmat,&((PC_Factor*)ilu)->info);
242: CHKMEMQ;
243: }
244: return(0);
245: }
249: static PetscErrorCode PCDestroy_ILU(PC pc)
250: {
251: PC_ILU *ilu = (PC_ILU*)pc->data;
255: PCReset_ILU(pc);
256: PetscFree(((PC_Factor*)ilu)->solvertype);
257: PetscFree(((PC_Factor*)ilu)->ordering);
258: PetscFree(pc->data);
259: return(0);
260: }
264: static PetscErrorCode PCApply_ILU(PC pc,Vec x,Vec y)
265: {
266: PC_ILU *ilu = (PC_ILU*)pc->data;
270: MatSolve(((PC_Factor*)ilu)->fact,x,y);
271: return(0);
272: }
276: static PetscErrorCode PCApplyTranspose_ILU(PC pc,Vec x,Vec y)
277: {
278: PC_ILU *ilu = (PC_ILU*)pc->data;
282: MatSolveTranspose(((PC_Factor*)ilu)->fact,x,y);
283: return(0);
284: }
288: static PetscErrorCode PCApplySymmetricLeft_ILU(PC pc,Vec x,Vec y)
289: {
291: PC_ILU *icc = (PC_ILU*)pc->data;
294: MatForwardSolve(((PC_Factor*)icc)->fact,x,y);
295: return(0);
296: }
300: static PetscErrorCode PCApplySymmetricRight_ILU(PC pc,Vec x,Vec y)
301: {
303: PC_ILU *icc = (PC_ILU*)pc->data;
306: MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);
307: return(0);
308: }
310: /*MC
311: PCILU - Incomplete factorization preconditioners.
313: Options Database Keys:
314: + -pc_factor_levels <k> - number of levels of fill for ILU(k)
315: . -pc_factor_in_place - only for ILU(0) with natural ordering, reuses the space of the matrix for
316: its factorization (overwrites original matrix)
317: . -pc_factor_diagonal_fill - fill in a zero diagonal even if levels of fill indicate it wouldn't be fill
318: . -pc_factor_reuse_ordering - reuse ordering of factorized matrix from previous factorization
319: . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
320: . -pc_factor_nonzeros_along_diagonal - reorder the matrix before factorization to remove zeros from the diagonal,
321: this decreases the chance of getting a zero pivot
322: . -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
323: . -pc_factor_pivot_in_blocks - for block ILU(k) factorization, i.e. with BAIJ matrices with block size larger
324: than 1 the diagonal blocks are factored with partial pivoting (this increases the
325: stability of the ILU factorization
327: Level: beginner
329: Concepts: incomplete factorization
331: Notes: Only implemented for some matrix formats. (for parallel see PCHYPRE for hypre's ILU)
333: For BAIJ matrices this implements a point block ILU
335: The "symmetric" application of this preconditioner is not actually symmetric since L is not transpose(U)
336: even when the matrix is not symmetric since the U stores the diagonals of the factorization.
338: If you are using MATSEQAIJCUSP matrices (or MATMPIAIJCUSP matrices with block Jacobi) you must have ./configured
339: PETSc with also --download-txpetscgpu to have the triangular solves performed on the GPU (factorization is never
340: done on the GPU).
342: References:
343: T. Dupont, R. Kendall, and H. Rachford. An approximate factorization procedure for solving
344: self-adjoint elliptic difference equations. SIAM J. Numer. Anal., 5:559--573, 1968.
346: T.A. Oliphant. An implicit numerical method for solving two-dimensional time-dependent dif-
347: fusion problems. Quart. Appl. Math., 19:221--229, 1961.
349: Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS, TONY F. CHAN AND HENK A. VAN DER VORST
350: http://igitur-archive.library.uu.nl/math/2001-0621-115821/proc.pdf chapter in Parallel Numerical
351: Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
352: Science and Engineering, Kluwer, pp. 167--202.
355: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
356: PCFactorSetZeroPivot(), PCFactorSetShiftSetType(), PCFactorSetAmount(),
357: PCFactorSetDropTolerance(),PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(),
358: PCFactorSetLevels(), PCFactorSetUseInPlace(), PCFactorSetAllowDiagonalFill(), PCFactorSetPivotInBlocks()
360: M*/
362: EXTERN_C_BEGIN
365: PetscErrorCode PCCreate_ILU(PC pc)
366: {
368: PC_ILU *ilu;
371: PetscNewLog(pc,PC_ILU,&ilu);
373: ((PC_Factor*)ilu)->fact = 0;
374: MatFactorInfoInitialize(&((PC_Factor*)ilu)->info);
375: ((PC_Factor*)ilu)->factortype = MAT_FACTOR_ILU;
376: ((PC_Factor*)ilu)->info.levels = 0.;
377: ((PC_Factor*)ilu)->info.fill = 1.0;
378: ilu->col = 0;
379: ilu->row = 0;
380: ilu->inplace = PETSC_FALSE;
381: PetscStrallocpy(MATSOLVERPETSC,&((PC_Factor*)ilu)->solvertype);
382: PetscStrallocpy(MATORDERINGNATURAL,&((PC_Factor*)ilu)->ordering);
383: ilu->reuseordering = PETSC_FALSE;
384: ((PC_Factor*)ilu)->info.dt = PETSC_DEFAULT;
385: ((PC_Factor*)ilu)->info.dtcount = PETSC_DEFAULT;
386: ((PC_Factor*)ilu)->info.dtcol = PETSC_DEFAULT;
387: ((PC_Factor*)ilu)->info.shifttype = (PetscReal)MAT_SHIFT_NONZERO;
388: ((PC_Factor*)ilu)->info.shiftamount = 100.0*PETSC_MACHINE_EPSILON;
389: ((PC_Factor*)ilu)->info.zeropivot = 100.0*PETSC_MACHINE_EPSILON;
390: ((PC_Factor*)ilu)->info.pivotinblocks = 1.0;
391: ilu->reusefill = PETSC_FALSE;
392: ((PC_Factor*)ilu)->info.diagonal_fill = 0.0;
393: pc->data = (void*)ilu;
395: pc->ops->reset = PCReset_ILU;
396: pc->ops->destroy = PCDestroy_ILU;
397: pc->ops->apply = PCApply_ILU;
398: pc->ops->applytranspose = PCApplyTranspose_ILU;
399: pc->ops->setup = PCSetUp_ILU;
400: pc->ops->setfromoptions = PCSetFromOptions_ILU;
401: pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor;
402: pc->ops->view = PCView_ILU;
403: pc->ops->applysymmetricleft = PCApplySymmetricLeft_ILU;
404: pc->ops->applysymmetricright = PCApplySymmetricRight_ILU;
405: pc->ops->applyrichardson = 0;
407: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Factor",
408: PCFactorSetZeroPivot_Factor);
409: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftType_C","PCFactorSetShiftType_Factor",
410: PCFactorSetShiftType_Factor);
411: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftAmount_C","PCFactorSetShiftAmount_Factor",
412: PCFactorSetShiftAmount_Factor);
413: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorGetMatSolverPackage_C","PCFactorGetMatSolverPackage_Factor",
414: PCFactorGetMatSolverPackage_Factor);
415: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatSolverPackage_C","PCFactorSetMatSolverPackage_Factor",
416: PCFactorSetMatSolverPackage_Factor);
417: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C","PCFactorSetUpMatSolverPackage_Factor",
418: PCFactorSetUpMatSolverPackage_Factor);
419: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetDropTolerance_C","PCFactorSetDropTolerance_ILU",
420: PCFactorSetDropTolerance_ILU);
421: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Factor",
422: PCFactorSetFill_Factor);
423: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Factor",
424: PCFactorSetMatOrderingType_Factor);
425: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseOrdering_C","PCFactorSetReuseOrdering_ILU",
426: PCFactorSetReuseOrdering_ILU);
427: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseFill_C","PCFactorSetReuseFill_ILU",
428: PCFactorSetReuseFill_ILU);
429: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetLevels_C","PCFactorSetLevels_Factor",
430: PCFactorSetLevels_Factor);
431: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUseInPlace_C","PCFactorSetUseInPlace_ILU",
432: PCFactorSetUseInPlace_ILU);
433: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetAllowDiagonalFill_C","PCFactorSetAllowDiagonalFill_Factor",
434: PCFactorSetAllowDiagonalFill_Factor);
435: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetPivotInBlocks_C","PCFactorSetPivotInBlocks_Factor",
436: PCFactorSetPivotInBlocks_Factor);
437: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C","PCFactorReorderForNonzeroDiagonal_ILU",
438: PCFactorReorderForNonzeroDiagonal_ILU);
439: return(0);
440: }
441: EXTERN_C_END