Actual source code: ilu.c

petsc-3.7.4 2016-10-02
Report Typos and Errors
  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: }