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: }