Actual source code: factimpl.c


  2: #include <../src/ksp/pc/impls/factor/factor.h>

  4: PetscErrorCode PCFactorSetUpMatSolverType_Factor(PC pc)
  5: {
  6:   PC_Factor *icc = (PC_Factor *)pc->data;

  8:   PetscFunctionBegin;
  9:   PetscCheck(pc->pmat, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "You can only call this routine after the matrix object has been provided to the solver, for example with KSPSetOperators() or SNESSetJacobian()");
 10:   if (!pc->setupcalled && !((PC_Factor *)icc)->fact) PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)icc)->solvertype, ((PC_Factor *)icc)->factortype, &((PC_Factor *)icc)->fact));
 11:   PetscFunctionReturn(PETSC_SUCCESS);
 12: }

 14: PetscErrorCode PCFactorSetZeroPivot_Factor(PC pc, PetscReal z)
 15: {
 16:   PC_Factor *ilu = (PC_Factor *)pc->data;

 18:   PetscFunctionBegin;
 19:   ilu->info.zeropivot = z;
 20:   PetscFunctionReturn(PETSC_SUCCESS);
 21: }

 23: PetscErrorCode PCFactorSetShiftType_Factor(PC pc, MatFactorShiftType shifttype)
 24: {
 25:   PC_Factor *dir = (PC_Factor *)pc->data;

 27:   PetscFunctionBegin;
 28:   if (shifttype == (MatFactorShiftType)PETSC_DECIDE) dir->info.shifttype = (PetscReal)MAT_SHIFT_NONE;
 29:   else {
 30:     dir->info.shifttype = (PetscReal)shifttype;
 31:     if ((shifttype == MAT_SHIFT_NONZERO || shifttype == MAT_SHIFT_INBLOCKS) && dir->info.shiftamount == 0.0) { dir->info.shiftamount = 100.0 * PETSC_MACHINE_EPSILON; /* set default amount if user has not called PCFactorSetShiftAmount() yet */ }
 32:   }
 33:   PetscFunctionReturn(PETSC_SUCCESS);
 34: }

 36: PetscErrorCode PCFactorSetShiftAmount_Factor(PC pc, PetscReal shiftamount)
 37: {
 38:   PC_Factor *dir = (PC_Factor *)pc->data;

 40:   PetscFunctionBegin;
 41:   if (shiftamount == (PetscReal)PETSC_DECIDE) dir->info.shiftamount = 100.0 * PETSC_MACHINE_EPSILON;
 42:   else dir->info.shiftamount = shiftamount;
 43:   PetscFunctionReturn(PETSC_SUCCESS);
 44: }

 46: PetscErrorCode PCFactorSetDropTolerance_Factor(PC pc, PetscReal dt, PetscReal dtcol, PetscInt dtcount)
 47: {
 48:   PC_Factor *ilu = (PC_Factor *)pc->data;

 50:   PetscFunctionBegin;
 51:   if (pc->setupcalled && (!ilu->info.usedt || ((PC_Factor *)ilu)->info.dt != dt || ((PC_Factor *)ilu)->info.dtcol != dtcol || ((PC_Factor *)ilu)->info.dtcount != dtcount)) {
 52:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot change tolerance after use");
 53:   }
 54:   ilu->info.usedt   = PETSC_TRUE;
 55:   ilu->info.dt      = dt;
 56:   ilu->info.dtcol   = dtcol;
 57:   ilu->info.dtcount = dtcount;
 58:   ilu->info.fill    = PETSC_DEFAULT;
 59:   PetscFunctionReturn(PETSC_SUCCESS);
 60: }

 62: PetscErrorCode PCFactorSetFill_Factor(PC pc, PetscReal fill)
 63: {
 64:   PC_Factor *dir = (PC_Factor *)pc->data;

 66:   PetscFunctionBegin;
 67:   dir->info.fill = fill;
 68:   PetscFunctionReturn(PETSC_SUCCESS);
 69: }

 71: PetscErrorCode PCFactorSetMatOrderingType_Factor(PC pc, MatOrderingType ordering)
 72: {
 73:   PC_Factor *dir = (PC_Factor *)pc->data;
 74:   PetscBool  flg;

 76:   PetscFunctionBegin;
 77:   if (!pc->setupcalled) {
 78:     PetscCall(PetscFree(dir->ordering));
 79:     PetscCall(PetscStrallocpy(ordering, (char **)&dir->ordering));
 80:   } else {
 81:     PetscCall(PetscStrcmp(dir->ordering, ordering, &flg));
 82:     PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot change ordering after use");
 83:   }
 84:   PetscFunctionReturn(PETSC_SUCCESS);
 85: }

 87: PetscErrorCode PCFactorGetLevels_Factor(PC pc, PetscInt *levels)
 88: {
 89:   PC_Factor *ilu = (PC_Factor *)pc->data;

 91:   PetscFunctionBegin;
 92:   *levels = ilu->info.levels;
 93:   PetscFunctionReturn(PETSC_SUCCESS);
 94: }

 96: PetscErrorCode PCFactorGetZeroPivot_Factor(PC pc, PetscReal *pivot)
 97: {
 98:   PC_Factor *ilu = (PC_Factor *)pc->data;

100:   PetscFunctionBegin;
101:   *pivot = ilu->info.zeropivot;
102:   PetscFunctionReturn(PETSC_SUCCESS);
103: }

105: PetscErrorCode PCFactorGetShiftAmount_Factor(PC pc, PetscReal *shift)
106: {
107:   PC_Factor *ilu = (PC_Factor *)pc->data;

109:   PetscFunctionBegin;
110:   *shift = ilu->info.shiftamount;
111:   PetscFunctionReturn(PETSC_SUCCESS);
112: }

114: PetscErrorCode PCFactorGetShiftType_Factor(PC pc, MatFactorShiftType *type)
115: {
116:   PC_Factor *ilu = (PC_Factor *)pc->data;

118:   PetscFunctionBegin;
119:   *type = (MatFactorShiftType)(int)ilu->info.shifttype;
120:   PetscFunctionReturn(PETSC_SUCCESS);
121: }

123: PetscErrorCode PCFactorSetLevels_Factor(PC pc, PetscInt levels)
124: {
125:   PC_Factor *ilu = (PC_Factor *)pc->data;

127:   PetscFunctionBegin;
128:   if (!pc->setupcalled) ilu->info.levels = levels;
129:   else if (ilu->info.levels != levels) {
130:     PetscUseTypeMethod(pc, reset); /* remove previous factored matrices */
131:     pc->setupcalled  = 0;          /* force a complete rebuild of preconditioner factored matrices */
132:     ilu->info.levels = levels;
133:   } else PetscCheck(!ilu->info.usedt, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot change levels after use with ILUdt");
134:   PetscFunctionReturn(PETSC_SUCCESS);
135: }

137: PetscErrorCode PCFactorSetAllowDiagonalFill_Factor(PC pc, PetscBool flg)
138: {
139:   PC_Factor *dir = (PC_Factor *)pc->data;

141:   PetscFunctionBegin;
142:   dir->info.diagonal_fill = (PetscReal)flg;
143:   PetscFunctionReturn(PETSC_SUCCESS);
144: }

146: PetscErrorCode PCFactorGetAllowDiagonalFill_Factor(PC pc, PetscBool *flg)
147: {
148:   PC_Factor *dir = (PC_Factor *)pc->data;

150:   PetscFunctionBegin;
151:   *flg = dir->info.diagonal_fill ? PETSC_TRUE : PETSC_FALSE;
152:   PetscFunctionReturn(PETSC_SUCCESS);
153: }

155: PetscErrorCode PCFactorSetPivotInBlocks_Factor(PC pc, PetscBool pivot)
156: {
157:   PC_Factor *dir = (PC_Factor *)pc->data;

159:   PetscFunctionBegin;
160:   dir->info.pivotinblocks = pivot ? 1.0 : 0.0;
161:   PetscFunctionReturn(PETSC_SUCCESS);
162: }

164: PetscErrorCode PCFactorGetMatrix_Factor(PC pc, Mat *mat)
165: {
166:   PC_Factor *ilu = (PC_Factor *)pc->data;

168:   PetscFunctionBegin;
169:   PetscCheck(ilu->fact, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Matrix not yet factored; call after KSPSetUp() or PCSetUp()");
170:   *mat = ilu->fact;
171:   PetscFunctionReturn(PETSC_SUCCESS);
172: }

174: /* allow access to preallocation information */
175: #include <petsc/private/matimpl.h>

177: PetscErrorCode PCFactorSetMatSolverType_Factor(PC pc, MatSolverType stype)
178: {
179:   PC_Factor *lu = (PC_Factor *)pc->data;

181:   PetscFunctionBegin;
182:   if (lu->fact && lu->fact->assembled) {
183:     MatSolverType ltype;
184:     PetscBool     flg;
185:     PetscCall(MatFactorGetSolverType(lu->fact, &ltype));
186:     PetscCall(PetscStrcmp(stype, ltype, &flg));
187:     PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot change solver matrix package from %s to %s after PC has been setup or used", ltype, stype);
188:   }

190:   PetscCall(PetscFree(lu->solvertype));
191:   PetscCall(PetscStrallocpy(stype, &lu->solvertype));
192:   PetscFunctionReturn(PETSC_SUCCESS);
193: }

195: PetscErrorCode PCFactorGetMatSolverType_Factor(PC pc, MatSolverType *stype)
196: {
197:   PC_Factor *lu = (PC_Factor *)pc->data;

199:   PetscFunctionBegin;
200:   *stype = lu->solvertype;
201:   PetscFunctionReturn(PETSC_SUCCESS);
202: }

204: PetscErrorCode PCFactorSetColumnPivot_Factor(PC pc, PetscReal dtcol)
205: {
206:   PC_Factor *dir = (PC_Factor *)pc->data;

208:   PetscFunctionBegin;
209:   PetscCheck(dtcol >= 0.0 && dtcol <= 1.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Column pivot tolerance is %g must be between 0 and 1", (double)dtcol);
210:   dir->info.dtcol = dtcol;
211:   PetscFunctionReturn(PETSC_SUCCESS);
212: }

214: PetscErrorCode PCSetFromOptions_Factor(PC pc, PetscOptionItems *PetscOptionsObject)
215: {
216:   PC_Factor        *factor = (PC_Factor *)pc->data;
217:   PetscBool         flg, set;
218:   char              tname[256], solvertype[64];
219:   PetscFunctionList ordlist;
220:   PetscEnum         etmp;
221:   PetscBool         inplace;

223:   PetscFunctionBegin;
224:   PetscCall(PCFactorGetUseInPlace(pc, &inplace));
225:   PetscCall(PetscOptionsBool("-pc_factor_in_place", "Form factored matrix in the same memory as the matrix", "PCFactorSetUseInPlace", inplace, &flg, &set));
226:   if (set) PetscCall(PCFactorSetUseInPlace(pc, flg));
227:   PetscCall(PetscOptionsReal("-pc_factor_fill", "Expected non-zeros in factored matrix", "PCFactorSetFill", ((PC_Factor *)factor)->info.fill, &((PC_Factor *)factor)->info.fill, NULL));

229:   PetscCall(PetscOptionsEnum("-pc_factor_shift_type", "Type of shift to add to diagonal", "PCFactorSetShiftType", MatFactorShiftTypes, (PetscEnum)(int)((PC_Factor *)factor)->info.shifttype, &etmp, &flg));
230:   if (flg) PetscCall(PCFactorSetShiftType(pc, (MatFactorShiftType)etmp));
231:   PetscCall(PetscOptionsReal("-pc_factor_shift_amount", "Shift added to diagonal", "PCFactorSetShiftAmount", ((PC_Factor *)factor)->info.shiftamount, &((PC_Factor *)factor)->info.shiftamount, NULL));

233:   PetscCall(PetscOptionsReal("-pc_factor_zeropivot", "Pivot is considered zero if less than", "PCFactorSetZeroPivot", ((PC_Factor *)factor)->info.zeropivot, &((PC_Factor *)factor)->info.zeropivot, NULL));
234:   PetscCall(PetscOptionsReal("-pc_factor_column_pivot", "Column pivot tolerance (used only for some factorization)", "PCFactorSetColumnPivot", ((PC_Factor *)factor)->info.dtcol, &((PC_Factor *)factor)->info.dtcol, &flg));

236:   PetscCall(PetscOptionsBool("-pc_factor_pivot_in_blocks", "Pivot inside matrix dense blocks for BAIJ and SBAIJ", "PCFactorSetPivotInBlocks", ((PC_Factor *)factor)->info.pivotinblocks ? PETSC_TRUE : PETSC_FALSE, &flg, &set));
237:   if (set) PetscCall(PCFactorSetPivotInBlocks(pc, flg));

239:   PetscCall(PetscOptionsBool("-pc_factor_reuse_fill", "Use fill from previous factorization", "PCFactorSetReuseFill", PETSC_FALSE, &flg, &set));
240:   if (set) PetscCall(PCFactorSetReuseFill(pc, flg));
241:   PetscCall(PetscOptionsBool("-pc_factor_reuse_ordering", "Reuse ordering from previous factorization", "PCFactorSetReuseOrdering", PETSC_FALSE, &flg, &set));
242:   if (set) PetscCall(PCFactorSetReuseOrdering(pc, flg));

244:   PetscCall(PetscOptionsDeprecated("-pc_factor_mat_solver_package", "-pc_factor_mat_solver_type", "3.9", NULL));
245:   PetscCall(PetscOptionsString("-pc_factor_mat_solver_type", "Specific direct solver to use", "MatGetFactor", ((PC_Factor *)factor)->solvertype, solvertype, sizeof(solvertype), &flg));
246:   if (flg) PetscCall(PCFactorSetMatSolverType(pc, solvertype));
247:   PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
248:   PetscCall(MatGetOrderingList(&ordlist));
249:   PetscCall(PetscOptionsFList("-pc_factor_mat_ordering_type", "Reordering to reduce nonzeros in factored matrix", "PCFactorSetMatOrderingType", ordlist, ((PC_Factor *)factor)->ordering, tname, sizeof(tname), &flg));
250:   if (flg) PetscCall(PCFactorSetMatOrderingType(pc, tname));
251:   PetscFunctionReturn(PETSC_SUCCESS);
252: }

254: PetscErrorCode PCView_Factor(PC pc, PetscViewer viewer)
255: {
256:   PC_Factor        *factor = (PC_Factor *)pc->data;
257:   PetscBool         isstring, iascii, canuseordering;
258:   MatInfo           info;
259:   MatOrderingType   ordering;
260:   PetscViewerFormat format;

262:   PetscFunctionBegin;
263:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
264:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
265:   if (iascii) {
266:     if (factor->inplace) {
267:       PetscCall(PetscViewerASCIIPrintf(viewer, "  in-place factorization\n"));
268:     } else {
269:       PetscCall(PetscViewerASCIIPrintf(viewer, "  out-of-place factorization\n"));
270:     }

272:     if (factor->reusefill) PetscCall(PetscViewerASCIIPrintf(viewer, "  Reusing fill from past factorization\n"));
273:     if (factor->reuseordering) PetscCall(PetscViewerASCIIPrintf(viewer, "  Reusing reordering from past factorization\n"));
274:     if (factor->factortype == MAT_FACTOR_ILU || factor->factortype == MAT_FACTOR_ICC) {
275:       if (factor->info.dt > 0) {
276:         PetscCall(PetscViewerASCIIPrintf(viewer, "  drop tolerance %g\n", (double)factor->info.dt));
277:         PetscCall(PetscViewerASCIIPrintf(viewer, "  max nonzeros per row %" PetscInt_FMT "\n", (PetscInt)factor->info.dtcount));
278:         PetscCall(PetscViewerASCIIPrintf(viewer, "  column permutation tolerance %g\n", (double)factor->info.dtcol));
279:       } else if (factor->info.levels == 1) {
280:         PetscCall(PetscViewerASCIIPrintf(viewer, "  %" PetscInt_FMT " level of fill\n", (PetscInt)factor->info.levels));
281:       } else {
282:         PetscCall(PetscViewerASCIIPrintf(viewer, "  %" PetscInt_FMT " levels of fill\n", (PetscInt)factor->info.levels));
283:       }
284:     }

286:     PetscCall(PetscViewerASCIIPrintf(viewer, "  tolerance for zero pivot %g\n", (double)factor->info.zeropivot));
287:     if (MatFactorShiftTypesDetail[(int)factor->info.shifttype]) { /* Only print when using a nontrivial shift */
288:       PetscCall(PetscViewerASCIIPrintf(viewer, "  using %s [%s]\n", MatFactorShiftTypesDetail[(int)factor->info.shifttype], MatFactorShiftTypes[(int)factor->info.shifttype]));
289:     }

291:     if (factor->fact) {
292:       PetscCall(MatFactorGetCanUseOrdering(factor->fact, &canuseordering));
293:       if (!canuseordering) ordering = MATORDERINGEXTERNAL;
294:       else ordering = factor->ordering;
295:       PetscCall(PetscViewerASCIIPrintf(viewer, "  matrix ordering: %s\n", ordering));
296:       if (!factor->fact->assembled) {
297:         PetscCall(PetscViewerASCIIPrintf(viewer, "  matrix solver type: %s\n", factor->fact->solvertype));
298:         PetscCall(PetscViewerASCIIPrintf(viewer, "  matrix not yet factored; no additional information available\n"));
299:       } else {
300:         PetscCall(MatGetInfo(factor->fact, MAT_LOCAL, &info));
301:         PetscCall(PetscViewerASCIIPrintf(viewer, "  factor fill ratio given %g, needed %g\n", (double)info.fill_ratio_given, (double)info.fill_ratio_needed));
302:         PetscCall(PetscViewerASCIIPrintf(viewer, "    Factored matrix follows:\n"));
303:         PetscCall(PetscViewerASCIIPushTab(viewer));
304:         PetscCall(PetscViewerASCIIPushTab(viewer));
305:         PetscCall(PetscViewerASCIIPushTab(viewer));
306:         PetscCall(PetscViewerGetFormat(viewer, &format));
307:         PetscCall(PetscViewerPushFormat(viewer, format != PETSC_VIEWER_ASCII_INFO_DETAIL ? PETSC_VIEWER_ASCII_INFO : PETSC_VIEWER_ASCII_INFO_DETAIL));
308:         PetscCall(MatView(factor->fact, viewer));
309:         PetscCall(PetscViewerPopFormat(viewer));
310:         PetscCall(PetscViewerASCIIPopTab(viewer));
311:         PetscCall(PetscViewerASCIIPopTab(viewer));
312:         PetscCall(PetscViewerASCIIPopTab(viewer));
313:       }
314:     }

316:   } else if (isstring) {
317:     MatFactorType t;
318:     PetscCall(MatGetFactorType(factor->fact, &t));
319:     if (t == MAT_FACTOR_ILU || t == MAT_FACTOR_ICC) PetscCall(PetscViewerStringSPrintf(viewer, " lvls=%" PetscInt_FMT ",order=%s", (PetscInt)factor->info.levels, factor->ordering));
320:   }
321:   PetscFunctionReturn(PETSC_SUCCESS);
322: }