Actual source code: ilu.c

petsc-master 2017-03-28
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>

  7: PetscErrorCode  PCFactorReorderForNonzeroDiagonal_ILU(PC pc,PetscReal z)
  8: {
  9:   PC_ILU *ilu = (PC_ILU*)pc->data;

 12:   ilu->nonzerosalongdiagonal = PETSC_TRUE;
 13:   if (z == PETSC_DECIDE) ilu->nonzerosalongdiagonaltol = 1.e-10;
 14:   else ilu->nonzerosalongdiagonaltol = z;
 15:   return(0);
 16: }

 18: PetscErrorCode PCReset_ILU(PC pc)
 19: {
 20:   PC_ILU         *ilu = (PC_ILU*)pc->data;

 24:   if (!ilu->hdr.inplace) {MatDestroy(&((PC_Factor*)ilu)->fact);}
 25:   if (ilu->row && ilu->col && ilu->row != ilu->col) {ISDestroy(&ilu->row);}
 26:   ISDestroy(&ilu->col);
 27:   return(0);
 28: }

 30: PetscErrorCode  PCFactorSetDropTolerance_ILU(PC pc,PetscReal dt,PetscReal dtcol,PetscInt dtcount)
 31: {
 32:   PC_ILU *ilu = (PC_ILU*)pc->data;

 35:   if (pc->setupcalled && (((PC_Factor*)ilu)->info.dt != dt || ((PC_Factor*)ilu)->info.dtcol != dtcol || ((PC_Factor*)ilu)->info.dtcount != dtcount)) {
 36:     SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot change drop tolerance after using PC");
 37:   }
 38:   ((PC_Factor*)ilu)->info.dt      = dt;
 39:   ((PC_Factor*)ilu)->info.dtcol   = dtcol;
 40:   ((PC_Factor*)ilu)->info.dtcount = dtcount;
 41:   ((PC_Factor*)ilu)->info.usedt   = 1.0;
 42:   return(0);
 43: }

 45: static PetscErrorCode PCSetFromOptions_ILU(PetscOptionItems *PetscOptionsObject,PC pc)
 46: {
 48:   PetscInt       itmp;
 49:   PetscBool      flg,set;
 50:   PC_ILU         *ilu = (PC_ILU*)pc->data;
 51:   PetscReal      tol;

 54:   PetscOptionsHead(PetscOptionsObject,"ILU Options");
 55:   PCSetFromOptions_Factor(PetscOptionsObject,pc);

 57:   PetscOptionsInt("-pc_factor_levels","levels of fill","PCFactorSetLevels",(PetscInt)((PC_Factor*)ilu)->info.levels,&itmp,&flg);
 58:   if (flg) ((PC_Factor*)ilu)->info.levels = itmp;

 60:   PetscOptionsBool("-pc_factor_diagonal_fill","Allow fill into empty diagonal entry","PCFactorSetAllowDiagonalFill",((PC_Factor*)ilu)->info.diagonal_fill ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
 61:   if (set) ((PC_Factor*)ilu)->info.diagonal_fill = (PetscReal) flg;
 62:   PetscOptionsName("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",&flg);
 63:   if (flg) {
 64:     tol  = PETSC_DECIDE;
 65:     PetscOptionsReal("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",ilu->nonzerosalongdiagonaltol,&tol,0);
 66:     PCFactorReorderForNonzeroDiagonal(pc,tol);
 67:   }

 69:   PetscOptionsTail();
 70:   return(0);
 71: }

 73: static PetscErrorCode PCView_ILU(PC pc,PetscViewer viewer)
 74: {

 78:   PCView_Factor(pc,viewer);
 79:   return(0);
 80: }

 82: static PetscErrorCode PCSetUp_ILU(PC pc)
 83: {
 84:   PetscErrorCode         ierr;
 85:   PC_ILU                 *ilu = (PC_ILU*)pc->data;
 86:   MatInfo                info;
 87:   PetscBool              flg;
 88:   const MatSolverPackage stype;
 89:   MatFactorError         err;

 92:   pc->failedreason = PC_NOERROR;
 93:   /* ugly hack to change default, since it is not support by some matrix types */
 94:   if (((PC_Factor*)ilu)->info.shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
 95:     PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJ,&flg);
 96:     if (!flg) {
 97:       PetscObjectTypeCompare((PetscObject)pc->pmat,MATMPIAIJ,&flg);
 98:       if (!flg) {
 99:         ((PC_Factor*)ilu)->info.shifttype = (PetscReal)MAT_SHIFT_INBLOCKS;
100:         PetscInfo(pc,"Changing shift type from NONZERO to INBLOCKS because block matrices do not support NONZERO\n");
101:       }
102:     }
103:   }

105:   MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
106:   if (ilu->hdr.inplace) {
107:     if (!pc->setupcalled) {

109:       /* In-place factorization only makes sense with the natural ordering,
110:          so we only need to get the ordering once, even if nonzero structure changes */
111:       MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);
112:       if (ilu->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->row);}
113:       if (ilu->col) {PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->col);}
114:     }

116:     /* In place ILU only makes sense with fill factor of 1.0 because
117:        cannot have levels of fill */
118:     ((PC_Factor*)ilu)->info.fill          = 1.0;
119:     ((PC_Factor*)ilu)->info.diagonal_fill = 0.0;

121:     MatILUFactor(pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);
122:     MatFactorGetError(pc->pmat,&err);
123:     if (err) { /* Factor() fails */
124:       pc->failedreason = (PCFailedReason)err;
125:       return(0);
126:     }

128:     ((PC_Factor*)ilu)->fact = pc->pmat;
129:     /* must update the pc record of the matrix state or the PC will attempt to run PCSetUp() yet again */
130:     PetscObjectStateGet((PetscObject)pc->pmat,&pc->matstate);
131:   } else {
132:     if (!pc->setupcalled) {
133:       /* first time in so compute reordering and symbolic factorization */
134:       MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);
135:       if (ilu->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->row);}
136:       if (ilu->col) {PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->col);}
137:       /*  Remove zeros along diagonal?     */
138:       if (ilu->nonzerosalongdiagonal) {
139:         MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);
140:       }
141:       if (!((PC_Factor*)ilu)->fact) {
142:         MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);
143:       }
144:       MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);
145:       MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);
146:       ilu->hdr.actualfill = info.fill_ratio_needed;

148:       PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)ilu)->fact);
149:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
150:       if (!ilu->hdr.reuseordering) {
151:         /* compute a new ordering for the ILU */
152:         ISDestroy(&ilu->row);
153:         ISDestroy(&ilu->col);
154:         MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);
155:         if (ilu->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->row);}
156:         if (ilu->col) {PetscLogObjectParent((PetscObject)pc,(PetscObject)ilu->col);}
157:         /*  Remove zeros along diagonal?     */
158:         if (ilu->nonzerosalongdiagonal) {
159:           MatReorderForNonzeroDiagonal(pc->pmat,ilu->nonzerosalongdiagonaltol,ilu->row,ilu->col);
160:         }
161:       }
162:       MatDestroy(&((PC_Factor*)ilu)->fact);
163:       MatGetFactor(pc->pmat,((PC_Factor*)ilu)->solvertype,MAT_FACTOR_ILU,&((PC_Factor*)ilu)->fact);
164:       MatILUFactorSymbolic(((PC_Factor*)ilu)->fact,pc->pmat,ilu->row,ilu->col,&((PC_Factor*)ilu)->info);
165:       MatGetInfo(((PC_Factor*)ilu)->fact,MAT_LOCAL,&info);
166:       ilu->hdr.actualfill = info.fill_ratio_needed;

168:       PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)ilu)->fact);
169:     }
170:     MatFactorGetError(((PC_Factor*)ilu)->fact,&err);
171:     if (err) { /* FactorSymbolic() fails */
172:       pc->failedreason = (PCFailedReason)err;
173:       return(0);
174:     }

176:     MatLUFactorNumeric(((PC_Factor*)ilu)->fact,pc->pmat,&((PC_Factor*)ilu)->info);
177:     MatFactorGetError(((PC_Factor*)ilu)->fact,&err);
178:     if (err) { /* FactorNumeric() fails */
179:       pc->failedreason = (PCFailedReason)err;
180:     }
181:   }

183:   PCFactorGetMatSolverPackage(pc,&stype);
184:   if (!stype) {
185:     const MatSolverPackage solverpackage;
186:     MatFactorGetSolverPackage(((PC_Factor*)ilu)->fact,&solverpackage);
187:     PCFactorSetMatSolverPackage(pc,solverpackage);
188:   }
189:   return(0);
190: }

192: static PetscErrorCode PCDestroy_ILU(PC pc)
193: {
194:   PC_ILU         *ilu = (PC_ILU*)pc->data;

198:   PCReset_ILU(pc);
199:   PetscFree(((PC_Factor*)ilu)->solvertype);
200:   PetscFree(((PC_Factor*)ilu)->ordering);
201:   PetscFree(pc->data);
202:   return(0);
203: }

205: static PetscErrorCode PCApply_ILU(PC pc,Vec x,Vec y)
206: {
207:   PC_ILU         *ilu = (PC_ILU*)pc->data;

211:   MatSolve(((PC_Factor*)ilu)->fact,x,y);
212:   return(0);
213: }

215: static PetscErrorCode PCApplyTranspose_ILU(PC pc,Vec x,Vec y)
216: {
217:   PC_ILU         *ilu = (PC_ILU*)pc->data;

221:   MatSolveTranspose(((PC_Factor*)ilu)->fact,x,y);
222:   return(0);
223: }

225: static PetscErrorCode PCApplySymmetricLeft_ILU(PC pc,Vec x,Vec y)
226: {
228:   PC_ILU         *icc = (PC_ILU*)pc->data;

231:   MatForwardSolve(((PC_Factor*)icc)->fact,x,y);
232:   return(0);
233: }

235: static PetscErrorCode PCApplySymmetricRight_ILU(PC pc,Vec x,Vec y)
236: {
238:   PC_ILU         *icc = (PC_ILU*)pc->data;

241:   MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);
242:   return(0);
243: }

245: /*MC
246:      PCILU - Incomplete factorization preconditioners.

248:    Options Database Keys:
249: +  -pc_factor_levels <k> - number of levels of fill for ILU(k)
250: .  -pc_factor_in_place - only for ILU(0) with natural ordering, reuses the space of the matrix for
251:                       its factorization (overwrites original matrix)
252: .  -pc_factor_diagonal_fill - fill in a zero diagonal even if levels of fill indicate it wouldn't be fill
253: .  -pc_factor_reuse_ordering - reuse ordering of factorized matrix from previous factorization
254: .  -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
255: .  -pc_factor_nonzeros_along_diagonal - reorder the matrix before factorization to remove zeros from the diagonal,
256:                                    this decreases the chance of getting a zero pivot
257: .  -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
258: -  -pc_factor_pivot_in_blocks - for block ILU(k) factorization, i.e. with BAIJ matrices with block size larger
259:                              than 1 the diagonal blocks are factored with partial pivoting (this increases the
260:                              stability of the ILU factorization

262:    Level: beginner

264:   Concepts: incomplete factorization

266:    Notes: Only implemented for some matrix formats. (for parallel see PCHYPRE for hypre's ILU)

268:           For BAIJ matrices this implements a point block ILU

270:           The "symmetric" application of this preconditioner is not actually symmetric since L is not transpose(U)
271:           even when the matrix is not symmetric since the U stores the diagonals of the factorization.

273:           If you are using MATSEQAIJCUSPARSE matrices (or MATMPIAIJCUSPARESE matrices with block Jacobi), factorization 
274:           is never done on the GPU).

276:    References:
277: +  1. - T. Dupont, R. Kendall, and H. Rachford. An approximate factorization procedure for solving
278:    self adjoint elliptic difference equations. SIAM J. Numer. Anal., 5, 1968.
279: .  2. -  T.A. Oliphant. An implicit numerical method for solving two dimensional timedependent diffusion problems. Quart. Appl. Math., 19, 1961.
280: -  3. -  TONY F. CHAN AND HENK A. VAN DER VORST, APPROXIMATE AND INCOMPLETE FACTORIZATIONS, 
281:       Chapter in Parallel Numerical
282:       Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
283:       Science and Engineering, Kluwer.


286: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
287:            PCFactorSetZeroPivot(), PCFactorSetShiftSetType(), PCFactorSetAmount(),
288:            PCFactorSetDropTolerance(),PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(),
289:            PCFactorSetLevels(), PCFactorSetUseInPlace(), PCFactorSetAllowDiagonalFill(), PCFactorSetPivotInBlocks(),
290:            PCFactorGetAllowDiagonalFill(), PCFactorGetUseInPlace()

292: M*/

294: PETSC_EXTERN PetscErrorCode PCCreate_ILU(PC pc)
295: {
297:   PC_ILU         *ilu;

300:   PetscNewLog(pc,&ilu);
301:   pc->data = (void*)ilu;
302:   PCFactorInitialize(pc);

304:   ((PC_Factor*)ilu)->factortype         = MAT_FACTOR_ILU;
305:   ((PC_Factor*)ilu)->info.levels        = 0.;
306:   ((PC_Factor*)ilu)->info.fill          = 1.0;
307:   ilu->col                              = 0;
308:   ilu->row                              = 0;
309:   PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)ilu)->ordering);
310:   ((PC_Factor*)ilu)->info.dt            = PETSC_DEFAULT;
311:   ((PC_Factor*)ilu)->info.dtcount       = PETSC_DEFAULT;
312:   ((PC_Factor*)ilu)->info.dtcol         = PETSC_DEFAULT;

314:   pc->ops->reset               = PCReset_ILU;
315:   pc->ops->destroy             = PCDestroy_ILU;
316:   pc->ops->apply               = PCApply_ILU;
317:   pc->ops->applytranspose      = PCApplyTranspose_ILU;
318:   pc->ops->setup               = PCSetUp_ILU;
319:   pc->ops->setfromoptions      = PCSetFromOptions_ILU;
320:   pc->ops->view                = PCView_ILU;
321:   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_ILU;
322:   pc->ops->applysymmetricright = PCApplySymmetricRight_ILU;
323:   pc->ops->applyrichardson     = 0;
324:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetDropTolerance_C",PCFactorSetDropTolerance_ILU);
325:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C",PCFactorReorderForNonzeroDiagonal_ILU);
326:   return(0);
327: }