Actual source code: ilu.c
petsc-3.7.5 2017-01-01
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,PetscBool flg)
77: {
78: PC_ILU *dir = (PC_ILU*)pc->data;
81: dir->inplace = flg;
82: return(0);
83: }
87: PetscErrorCode PCFactorGetUseInPlace_ILU(PC pc,PetscBool *flg)
88: {
89: PC_ILU *dir = (PC_ILU*)pc->data;
92: *flg = dir->inplace;
93: return(0);
94: }
98: static PetscErrorCode PCSetFromOptions_ILU(PetscOptionItems *PetscOptionsObject,PC pc)
99: {
101: PetscInt itmp;
102: PetscBool flg,set;
103: PC_ILU *ilu = (PC_ILU*)pc->data;
104: PetscReal tol;
107: PetscOptionsHead(PetscOptionsObject,"ILU Options");
108: PCSetFromOptions_Factor(PetscOptionsObject,pc);
110: PetscOptionsInt("-pc_factor_levels","levels of fill","PCFactorSetLevels",(PetscInt)((PC_Factor*)ilu)->info.levels,&itmp,&flg);
111: if (flg) ((PC_Factor*)ilu)->info.levels = itmp;
113: PetscOptionsBool("-pc_factor_diagonal_fill","Allow fill into empty diagonal entry","PCFactorSetAllowDiagonalFill",((PC_Factor*)ilu)->info.diagonal_fill ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
114: if (set) ((PC_Factor*)ilu)->info.diagonal_fill = (PetscReal) flg;
115: PetscOptionsName("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",&flg);
116: if (flg) {
117: tol = PETSC_DECIDE;
118: PetscOptionsReal("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",ilu->nonzerosalongdiagonaltol,&tol,0);
119: PCFactorReorderForNonzeroDiagonal(pc,tol);
120: }
122: PetscOptionsTail();
123: return(0);
124: }
128: static PetscErrorCode PCView_ILU(PC pc,PetscViewer viewer)
129: {
130: PC_ILU *ilu = (PC_ILU*)pc->data;
132: PetscBool iascii;
135: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
136: if (iascii) {
137: if (ilu->inplace) {
138: PetscViewerASCIIPrintf(viewer," ILU: in-place factorization\n");
139: } else {
140: PetscViewerASCIIPrintf(viewer," ILU: out-of-place factorization\n");
141: }
143: if (ilu->reusefill) {PetscViewerASCIIPrintf(viewer," ILU: Reusing fill from past factorization\n");}
144: if (ilu->reuseordering) {PetscViewerASCIIPrintf(viewer," ILU: Reusing reordering from past factorization\n");}
145: }
146: PCView_Factor(pc,viewer);
147: return(0);
148: }
152: static PetscErrorCode PCSetUp_ILU(PC pc)
153: {
155: PC_ILU *ilu = (PC_ILU*)pc->data;
156: MatInfo info;
157: PetscBool flg;
158: const MatSolverPackage stype;
161: pc->failedreason = PC_NOERROR;
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\n");
170: }
171: }
172: }
174: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
175: if (ilu->inplace) {
176: if (!pc->setupcalled) {
178: /* In-place factorization only makes sense with the natural ordering,
179: so we only need to get the ordering once, even if nonzero structure changes */
180: MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);
181: if (ilu->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->row);}
182: if (ilu->col) {PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->col);}
183: }
185: /* In place ILU only makes sense with fill factor of 1.0 because
186: cannot have levels of fill */
187: ((PC_Factor*)ilu)->info.fill = 1.0;
188: ((PC_Factor*)ilu)->info.diagonal_fill = 0.0;
190: MatILUFactor(pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);
191: if (pc->pmat->errortype) { /* Factor() fails */
192: pc->failedreason = (PCFailedReason)pc->pmat->errortype;
193: return(0);
194: }
196: ((PC_Factor*)ilu)->fact = pc->pmat;
197: /* must update the pc record of the matrix state or the PC will attempt to run PCSetUp() yet again */
198: PetscObjectStateGet((PetscObject)pc->pmat,&pc->matstate);
199: } else {
200: Mat F;
201: if (!pc->setupcalled) {
202: /* first time in so compute reordering and symbolic factorization */
203: MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);
204: if (ilu->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->row);}
205: if (ilu->col) {PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->col);}
206: /* Remove zeros along diagonal? */
207: if (ilu->nonzerosalongdiagonal) {
208: MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);
209: }
210: if (!((PC_Factor*)ilu)->fact) {
211: MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);
212: }
213: MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);
214: MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);
216: ilu->actualfill = info.fill_ratio_needed;
218: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)ilu)->fact);
219: } else if (pc->flag != SAME_NONZERO_PATTERN) {
220: if (!ilu->reuseordering) {
221: /* compute a new ordering for the ILU */
222: ISDestroy(&ilu->row);
223: ISDestroy(&ilu->col);
224: MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);
225: if (ilu->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->row);}
226: if (ilu->col) {PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->col);}
227: /* Remove zeros along diagonal? */
228: if (ilu->nonzerosalongdiagonal) {
229: MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);
230: }
231: }
232: MatDestroy(&((PC_Factor*)ilu)->fact);
233: MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);
234: MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);
235: MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);
237: ilu->actualfill = info.fill_ratio_needed;
239: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)ilu)->fact);
240: }
241: F = ((PC_Factor*)ilu)->fact;
242: if (F->errortype) { /* FactorSymbolic() fails */
243: pc->failedreason = (PCFailedReason)F->errortype;
244: return(0);
245: }
247: MatLUFactorNumeric(((PC_Factor*)ilu)->fact,pc->pmat,&((PC_Factor*)ilu)->info);
248: if (F->errortype) { /* FactorNumeric() fails */
249: pc->failedreason = (PCFailedReason)F->errortype;
250: }
251: }
253: PCFactorGetMatSolverPackage(pc,&stype);
254: if (!stype) {
255: PCFactorSetMatSolverPackage(pc,((PC_Factor*)ilu)->fact->solvertype);
256: }
257: return(0);
258: }
262: static PetscErrorCode PCDestroy_ILU(PC pc)
263: {
264: PC_ILU *ilu = (PC_ILU*)pc->data;
268: PCReset_ILU(pc);
269: PetscFree(((PC_Factor*)ilu)->solvertype);
270: PetscFree(((PC_Factor*)ilu)->ordering);
271: PetscFree(pc->data);
272: return(0);
273: }
277: static PetscErrorCode PCApply_ILU(PC pc,Vec x,Vec y)
278: {
279: PC_ILU *ilu = (PC_ILU*)pc->data;
283: MatSolve(((PC_Factor*)ilu)->fact,x,y);
284: return(0);
285: }
289: static PetscErrorCode PCApplyTranspose_ILU(PC pc,Vec x,Vec y)
290: {
291: PC_ILU *ilu = (PC_ILU*)pc->data;
295: MatSolveTranspose(((PC_Factor*)ilu)->fact,x,y);
296: return(0);
297: }
301: static PetscErrorCode PCApplySymmetricLeft_ILU(PC pc,Vec x,Vec y)
302: {
304: PC_ILU *icc = (PC_ILU*)pc->data;
307: MatForwardSolve(((PC_Factor*)icc)->fact,x,y);
308: return(0);
309: }
313: static PetscErrorCode PCApplySymmetricRight_ILU(PC pc,Vec x,Vec y)
314: {
316: PC_ILU *icc = (PC_ILU*)pc->data;
319: MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);
320: return(0);
321: }
323: /*MC
324: PCILU - Incomplete factorization preconditioners.
326: Options Database Keys:
327: + -pc_factor_levels <k> - number of levels of fill for ILU(k)
328: . -pc_factor_in_place - only for ILU(0) with natural ordering, reuses the space of the matrix for
329: its factorization (overwrites original matrix)
330: . -pc_factor_diagonal_fill - fill in a zero diagonal even if levels of fill indicate it wouldn't be fill
331: . -pc_factor_reuse_ordering - reuse ordering of factorized matrix from previous factorization
332: . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
333: . -pc_factor_nonzeros_along_diagonal - reorder the matrix before factorization to remove zeros from the diagonal,
334: this decreases the chance of getting a zero pivot
335: . -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
336: - -pc_factor_pivot_in_blocks - for block ILU(k) factorization, i.e. with BAIJ matrices with block size larger
337: than 1 the diagonal blocks are factored with partial pivoting (this increases the
338: stability of the ILU factorization
340: Level: beginner
342: Concepts: incomplete factorization
344: Notes: Only implemented for some matrix formats. (for parallel see PCHYPRE for hypre's ILU)
346: For BAIJ matrices this implements a point block ILU
348: The "symmetric" application of this preconditioner is not actually symmetric since L is not transpose(U)
349: even when the matrix is not symmetric since the U stores the diagonals of the factorization.
351: If you are using MATSEQAIJCUSPARSE matrices (or MATMPIAIJCUSPARESE matrices with block Jacobi), factorization
352: is never done on the GPU).
354: References:
355: + 1. - T. Dupont, R. Kendall, and H. Rachford. An approximate factorization procedure for solving
356: self adjoint elliptic difference equations. SIAM J. Numer. Anal., 5, 1968.
357: . 2. - T.A. Oliphant. An implicit numerical method for solving two dimensional timedependent diffusion problems. Quart. Appl. Math., 19, 1961.
358: - 3. - TONY F. CHAN AND HENK A. VAN DER VORST, APPROXIMATE AND INCOMPLETE FACTORIZATIONS,
359: Chapter in Parallel Numerical
360: Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
361: Science and Engineering, Kluwer.
364: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
365: PCFactorSetZeroPivot(), PCFactorSetShiftSetType(), PCFactorSetAmount(),
366: PCFactorSetDropTolerance(),PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(),
367: PCFactorSetLevels(), PCFactorSetUseInPlace(), PCFactorSetAllowDiagonalFill(), PCFactorSetPivotInBlocks(),
368: PCFactorGetAllowDiagonalFill(), PCFactorGetUseInPlace()
370: M*/
374: PETSC_EXTERN PetscErrorCode PCCreate_ILU(PC pc)
375: {
377: PC_ILU *ilu;
380: PetscNewLog(pc,&ilu);
382: ((PC_Factor*)ilu)->fact = 0;
383: MatFactorInfoInitialize(&((PC_Factor*)ilu)->info);
384: ((PC_Factor*)ilu)->factortype = MAT_FACTOR_ILU;
385: ((PC_Factor*)ilu)->info.levels = 0.;
386: ((PC_Factor*)ilu)->info.fill = 1.0;
387: ilu->col = 0;
388: ilu->row = 0;
389: ilu->inplace = PETSC_FALSE;
390: PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)ilu)->ordering);
391: ilu->reuseordering = PETSC_FALSE;
392: ((PC_Factor*)ilu)->info.dt = PETSC_DEFAULT;
393: ((PC_Factor*)ilu)->info.dtcount = PETSC_DEFAULT;
394: ((PC_Factor*)ilu)->info.dtcol = PETSC_DEFAULT;
395: ((PC_Factor*)ilu)->info.shifttype = (PetscReal)MAT_SHIFT_NONE;
396: ((PC_Factor*)ilu)->info.shiftamount = 100.0*PETSC_MACHINE_EPSILON;
397: ((PC_Factor*)ilu)->info.zeropivot = 100.0*PETSC_MACHINE_EPSILON;
398: ((PC_Factor*)ilu)->info.pivotinblocks = 1.0;
399: ilu->reusefill = PETSC_FALSE;
400: ((PC_Factor*)ilu)->info.diagonal_fill = 0.0;
401: pc->data = (void*)ilu;
403: pc->ops->reset = PCReset_ILU;
404: pc->ops->destroy = PCDestroy_ILU;
405: pc->ops->apply = PCApply_ILU;
406: pc->ops->applytranspose = PCApplyTranspose_ILU;
407: pc->ops->setup = PCSetUp_ILU;
408: pc->ops->setfromoptions = PCSetFromOptions_ILU;
409: pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor;
410: pc->ops->view = PCView_ILU;
411: pc->ops->applysymmetricleft = PCApplySymmetricLeft_ILU;
412: pc->ops->applysymmetricright = PCApplySymmetricRight_ILU;
413: pc->ops->applyrichardson = 0;
415: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetZeroPivot_C",PCFactorSetZeroPivot_Factor);
416: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetZeroPivot_C",PCFactorGetZeroPivot_Factor);
417: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Factor);
418: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetShiftType_C",PCFactorGetShiftType_Factor);
419: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftAmount_C",PCFactorSetShiftAmount_Factor);
420: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetShiftAmount_C",PCFactorGetShiftAmount_Factor);
421: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetMatSolverPackage_C",PCFactorGetMatSolverPackage_Factor);
422: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatSolverPackage_C",PCFactorSetMatSolverPackage_Factor);
423: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C",PCFactorSetUpMatSolverPackage_Factor);
424: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetDropTolerance_C",PCFactorSetDropTolerance_ILU);
425: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetFill_C",PCFactorSetFill_Factor);
426: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatOrderingType_C",PCFactorSetMatOrderingType_Factor);
427: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C",PCFactorSetReuseOrdering_ILU);
428: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseFill_C",PCFactorSetReuseFill_ILU);
429: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetLevels_C",PCFactorSetLevels_Factor);
430: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetLevels_C",PCFactorGetLevels_Factor);
431: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUseInPlace_C",PCFactorSetUseInPlace_ILU);
432: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetUseInPlace_C",PCFactorGetUseInPlace_ILU);
433: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetAllowDiagonalFill_C",PCFactorSetAllowDiagonalFill_Factor);
434: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetAllowDiagonalFill_C",PCFactorGetAllowDiagonalFill_Factor);
435: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetPivotInBlocks_C",PCFactorSetPivotInBlocks_Factor);
436: PetscObjectComposeFunction((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C",PCFactorReorderForNonzeroDiagonal_ILU);
437: return(0);
438: }