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(PC pc, PetscOptionItems *PetscOptionsObject)
 42: {
 43:   PetscInt  itmp;
 44:   PetscBool flg, set;
 45:   PC_ILU   *ilu = (PC_ILU *)pc->data;
 46:   PetscReal tol;

 48:   PetscOptionsHeadBegin(PetscOptionsObject, "ILU Options");
 49:   PCSetFromOptions_Factor(pc, PetscOptionsObject);

 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:   PetscOptionsHeadEnd();
 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;
 74:   const char    *prefix;

 76:   pc->failedreason = PC_NOERROR;
 77:   /* ugly hack to change default, since it is not support by some matrix types */
 78:   if (((PC_Factor *)ilu)->info.shifttype == (PetscReal)MAT_SHIFT_NONZERO) {
 79:     PetscObjectTypeCompare((PetscObject)pc->pmat, MATSEQAIJ, &flg);
 80:     if (!flg) {
 81:       PetscObjectTypeCompare((PetscObject)pc->pmat, MATMPIAIJ, &flg);
 82:       if (!flg) {
 83:         ((PC_Factor *)ilu)->info.shifttype = (PetscReal)MAT_SHIFT_INBLOCKS;
 84:         PetscInfo(pc, "Changing shift type from NONZERO to INBLOCKS because block matrices do not support NONZERO\n");
 85:       }
 86:     }
 87:   }

 89:   PCGetOptionsPrefix(pc, &prefix);
 90:   MatSetOptionsPrefixFactor(pc->pmat, prefix);

 92:   MatSetErrorIfFailure(pc->pmat, pc->erroriffailure);
 93:   if (ilu->hdr.inplace) {
 94:     if (!pc->setupcalled) {
 95:       /* In-place factorization only makes sense with the natural ordering,
 96:          so we only need to get the ordering once, even if nonzero structure changes */
 97:       /* Should not get the ordering if the factorization routine does not use it, but do not yet have access to the factor matrix */
 98:       PCFactorSetDefaultOrdering_Factor(pc);
 99:       MatDestroy(&((PC_Factor *)ilu)->fact);
100:       MatGetOrdering(pc->pmat, ((PC_Factor *)ilu)->ordering, &ilu->row, &ilu->col);
101:     }

103:     /* In place ILU only makes sense with fill factor of 1.0 because
104:        cannot have levels of fill */
105:     ((PC_Factor *)ilu)->info.fill          = 1.0;
106:     ((PC_Factor *)ilu)->info.diagonal_fill = 0.0;

108:     MatILUFactor(pc->pmat, ilu->row, ilu->col, &((PC_Factor *)ilu)->info);
109:     MatFactorGetError(pc->pmat, &err);
110:     if (err) { /* Factor() fails */
111:       pc->failedreason = (PCFailedReason)err;
112:       return 0;
113:     }

115:     ((PC_Factor *)ilu)->fact = pc->pmat;
116:     /* must update the pc record of the matrix state or the PC will attempt to run PCSetUp() yet again */
117:     PetscObjectStateGet((PetscObject)pc->pmat, &pc->matstate);
118:   } else {
119:     if (!pc->setupcalled) {
120:       /* first time in so compute reordering and symbolic factorization */
121:       PetscBool canuseordering;
122:       if (!((PC_Factor *)ilu)->fact) { MatGetFactor(pc->pmat, ((PC_Factor *)ilu)->solvertype, MAT_FACTOR_ILU, &((PC_Factor *)ilu)->fact); }
123:       MatFactorGetCanUseOrdering(((PC_Factor *)ilu)->fact, &canuseordering);
124:       if (canuseordering) {
125:         PCFactorSetDefaultOrdering_Factor(pc);
126:         MatGetOrdering(pc->pmat, ((PC_Factor *)ilu)->ordering, &ilu->row, &ilu->col);
127:         /*  Remove zeros along diagonal?     */
128:         if (ilu->nonzerosalongdiagonal) MatReorderForNonzeroDiagonal(pc->pmat, ilu->nonzerosalongdiagonaltol, ilu->row, ilu->col);
129:       }
130:       MatILUFactorSymbolic(((PC_Factor *)ilu)->fact, pc->pmat, ilu->row, ilu->col, &((PC_Factor *)ilu)->info);
131:       MatGetInfo(((PC_Factor *)ilu)->fact, MAT_LOCAL, &info);
132:       ilu->hdr.actualfill = info.fill_ratio_needed;
133:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
134:       if (!ilu->hdr.reuseordering) {
135:         PetscBool canuseordering;
136:         MatDestroy(&((PC_Factor *)ilu)->fact);
137:         MatGetFactor(pc->pmat, ((PC_Factor *)ilu)->solvertype, MAT_FACTOR_ILU, &((PC_Factor *)ilu)->fact);
138:         MatFactorGetCanUseOrdering(((PC_Factor *)ilu)->fact, &canuseordering);
139:         if (canuseordering) {
140:           /* compute a new ordering for the ILU */
141:           ISDestroy(&ilu->row);
142:           ISDestroy(&ilu->col);
143:           PCFactorSetDefaultOrdering_Factor(pc);
144:           MatGetOrdering(pc->pmat, ((PC_Factor *)ilu)->ordering, &ilu->row, &ilu->col);
145:           /*  Remove zeros along diagonal?     */
146:           if (ilu->nonzerosalongdiagonal) MatReorderForNonzeroDiagonal(pc->pmat, ilu->nonzerosalongdiagonaltol, ilu->row, ilu->col);
147:         }
148:       }
149:       MatILUFactorSymbolic(((PC_Factor *)ilu)->fact, pc->pmat, ilu->row, ilu->col, &((PC_Factor *)ilu)->info);
150:       MatGetInfo(((PC_Factor *)ilu)->fact, MAT_LOCAL, &info);
151:       ilu->hdr.actualfill = info.fill_ratio_needed;
152:     }
153:     MatFactorGetError(((PC_Factor *)ilu)->fact, &err);
154:     if (err) { /* FactorSymbolic() fails */
155:       pc->failedreason = (PCFailedReason)err;
156:       return 0;
157:     }

159:     MatLUFactorNumeric(((PC_Factor *)ilu)->fact, pc->pmat, &((PC_Factor *)ilu)->info);
160:     MatFactorGetError(((PC_Factor *)ilu)->fact, &err);
161:     if (err) { /* FactorNumeric() fails */
162:       pc->failedreason = (PCFailedReason)err;
163:     }
164:   }

166:   PCFactorGetMatSolverType(pc, &stype);
167:   if (!stype) {
168:     MatSolverType solverpackage;
169:     MatFactorGetSolverType(((PC_Factor *)ilu)->fact, &solverpackage);
170:     PCFactorSetMatSolverType(pc, solverpackage);
171:   }
172:   return 0;
173: }

175: static PetscErrorCode PCDestroy_ILU(PC pc)
176: {
177:   PC_ILU *ilu = (PC_ILU *)pc->data;

179:   PCReset_ILU(pc);
180:   PetscFree(((PC_Factor *)ilu)->solvertype);
181:   PetscFree(((PC_Factor *)ilu)->ordering);
182:   PetscFree(pc->data);
183:   PCFactorClearComposedFunctions(pc);
184:   return 0;
185: }

187: static PetscErrorCode PCApply_ILU(PC pc, Vec x, Vec y)
188: {
189:   PC_ILU *ilu = (PC_ILU *)pc->data;

191:   MatSolve(((PC_Factor *)ilu)->fact, x, y);
192:   return 0;
193: }

195: static PetscErrorCode PCMatApply_ILU(PC pc, Mat X, Mat Y)
196: {
197:   PC_ILU *ilu = (PC_ILU *)pc->data;

199:   MatMatSolve(((PC_Factor *)ilu)->fact, X, Y);
200:   return 0;
201: }

203: static PetscErrorCode PCApplyTranspose_ILU(PC pc, Vec x, Vec y)
204: {
205:   PC_ILU *ilu = (PC_ILU *)pc->data;

207:   MatSolveTranspose(((PC_Factor *)ilu)->fact, x, y);
208:   return 0;
209: }

211: static PetscErrorCode PCApplySymmetricLeft_ILU(PC pc, Vec x, Vec y)
212: {
213:   PC_ILU *icc = (PC_ILU *)pc->data;

215:   MatForwardSolve(((PC_Factor *)icc)->fact, x, y);
216:   return 0;
217: }

219: static PetscErrorCode PCApplySymmetricRight_ILU(PC pc, Vec x, Vec y)
220: {
221:   PC_ILU *icc = (PC_ILU *)pc->data;

223:   MatBackwardSolve(((PC_Factor *)icc)->fact, x, y);
224:   return 0;
225: }

227: /*MC
228:      PCILU - Incomplete factorization preconditioners.

230:    Options Database Keys:
231: +  -pc_factor_levels <k> - number of levels of fill for ILU(k)
232: .  -pc_factor_in_place - only for ILU(0) with natural ordering, reuses the space of the matrix for
233:                       its factorization (overwrites original matrix)
234: .  -pc_factor_diagonal_fill - fill in a zero diagonal even if levels of fill indicate it wouldn't be fill
235: .  -pc_factor_reuse_ordering - reuse ordering of factorized matrix from previous factorization
236: .  -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
237: .  -pc_factor_nonzeros_along_diagonal - reorder the matrix before factorization to remove zeros from the diagonal,
238:                                    this decreases the chance of getting a zero pivot
239: .  -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
240: -  -pc_factor_pivot_in_blocks - for block ILU(k) factorization, i.e. with BAIJ matrices with block size larger
241:                              than 1 the diagonal blocks are factored with partial pivoting (this increases the
242:                              stability of the ILU factorization

244:    Level: beginner

246:    Notes:
247:    Only implemented for some matrix format and sequential. For parallel see `PCHYPRE` for hypre's ILU

249:    For `MATSEQBAIJ` matrices this implements a point block ILU

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

254:    If you are using `MATSEQAIJCUSPARSE` matrices (or `MATMPIAIJCUSPARSE` matrices with block Jacobi), factorization
255:    is never done on the GPU).

257:    References:
258: +  * - T. Dupont, R. Kendall, and H. Rachford. An approximate factorization procedure for solving
259:    self adjoint elliptic difference equations. SIAM J. Numer. Anal., 5, 1968.
260: .  * -  T.A. Oliphant. An implicit numerical method for solving two dimensional timedependent diffusion problems. Quart. Appl. Math., 19, 1961.
261: -  * -  TONY F. CHAN AND HENK A. VAN DER VORST, APPROXIMATE AND INCOMPLETE FACTORIZATIONS,
262:       Chapter in Parallel Numerical
263:       Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
264:       Science and Engineering, Kluwer.

266: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSOR`, `MatOrderingType`, `PCLU`, `PCICC`, `PCCHOLESKY`,
267:           `PCFactorSetZeroPivot()`, `PCFactorSetShiftSetType()`, `PCFactorSetAmount()`,
268:           `PCFactorSetDropTolerance()`, `PCFactorSetFill()`, `PCFactorSetMatOrderingType()`, `PCFactorSetReuseOrdering()`,
269:           `PCFactorSetLevels()`, `PCFactorSetUseInPlace()`, `PCFactorSetAllowDiagonalFill()`, `PCFactorSetPivotInBlocks()`,
270:           `PCFactorGetAllowDiagonalFill()`, `PCFactorGetUseInPlace()`
271: M*/

273: PETSC_EXTERN PetscErrorCode PCCreate_ILU(PC pc)
274: {
275:   PC_ILU *ilu;

277:   PetscNew(&ilu);
278:   pc->data = (void *)ilu;
279:   PCFactorInitialize(pc, MAT_FACTOR_ILU);

281:   ((PC_Factor *)ilu)->info.levels  = 0.;
282:   ((PC_Factor *)ilu)->info.fill    = 1.0;
283:   ilu->col                         = NULL;
284:   ilu->row                         = NULL;
285:   ((PC_Factor *)ilu)->info.dt      = PETSC_DEFAULT;
286:   ((PC_Factor *)ilu)->info.dtcount = PETSC_DEFAULT;
287:   ((PC_Factor *)ilu)->info.dtcol   = PETSC_DEFAULT;

289:   pc->ops->reset               = PCReset_ILU;
290:   pc->ops->destroy             = PCDestroy_ILU;
291:   pc->ops->apply               = PCApply_ILU;
292:   pc->ops->matapply            = PCMatApply_ILU;
293:   pc->ops->applytranspose      = PCApplyTranspose_ILU;
294:   pc->ops->setup               = PCSetUp_ILU;
295:   pc->ops->setfromoptions      = PCSetFromOptions_ILU;
296:   pc->ops->view                = PCView_Factor;
297:   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_ILU;
298:   pc->ops->applysymmetricright = PCApplySymmetricRight_ILU;
299:   pc->ops->applyrichardson     = NULL;
300:   PetscObjectComposeFunction((PetscObject)pc, "PCFactorSetDropTolerance_C", PCFactorSetDropTolerance_ILU);
301:   PetscObjectComposeFunction((PetscObject)pc, "PCFactorReorderForNonzeroDiagonal_C", PCFactorReorderForNonzeroDiagonal_ILU);
302:   return 0;
303: }