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;

  9:   if (!pc->setupcalled && !((PC_Factor *)icc)->fact) MatGetFactor(pc->pmat, ((PC_Factor *)icc)->solvertype, ((PC_Factor *)icc)->factortype, &((PC_Factor *)icc)->fact);
 10:   return 0;
 11: }

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

 17:   ilu->info.zeropivot = z;
 18:   return 0;
 19: }

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

 25:   if (shifttype == (MatFactorShiftType)PETSC_DECIDE) dir->info.shifttype = (PetscReal)MAT_SHIFT_NONE;
 26:   else {
 27:     dir->info.shifttype = (PetscReal)shifttype;
 28:     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 */ }
 29:   }
 30:   return 0;
 31: }

 33: PetscErrorCode PCFactorSetShiftAmount_Factor(PC pc, PetscReal shiftamount)
 34: {
 35:   PC_Factor *dir = (PC_Factor *)pc->data;

 37:   if (shiftamount == (PetscReal)PETSC_DECIDE) dir->info.shiftamount = 100.0 * PETSC_MACHINE_EPSILON;
 38:   else dir->info.shiftamount = shiftamount;
 39:   return 0;
 40: }

 42: PetscErrorCode PCFactorSetDropTolerance_Factor(PC pc, PetscReal dt, PetscReal dtcol, PetscInt dtcount)
 43: {
 44:   PC_Factor *ilu = (PC_Factor *)pc->data;

 46:   if (pc->setupcalled && (!ilu->info.usedt || ((PC_Factor *)ilu)->info.dt != dt || ((PC_Factor *)ilu)->info.dtcol != dtcol || ((PC_Factor *)ilu)->info.dtcount != dtcount)) {
 47:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot change tolerance after use");
 48:   }
 49:   ilu->info.usedt   = PETSC_TRUE;
 50:   ilu->info.dt      = dt;
 51:   ilu->info.dtcol   = dtcol;
 52:   ilu->info.dtcount = dtcount;
 53:   ilu->info.fill    = PETSC_DEFAULT;
 54:   return 0;
 55: }

 57: PetscErrorCode PCFactorSetFill_Factor(PC pc, PetscReal fill)
 58: {
 59:   PC_Factor *dir = (PC_Factor *)pc->data;

 61:   dir->info.fill = fill;
 62:   return 0;
 63: }

 65: PetscErrorCode PCFactorSetMatOrderingType_Factor(PC pc, MatOrderingType ordering)
 66: {
 67:   PC_Factor *dir = (PC_Factor *)pc->data;
 68:   PetscBool  flg;

 70:   if (!pc->setupcalled) {
 71:     PetscFree(dir->ordering);
 72:     PetscStrallocpy(ordering, (char **)&dir->ordering);
 73:   } else {
 74:     PetscStrcmp(dir->ordering, ordering, &flg);
 76:   }
 77:   return 0;
 78: }

 80: PetscErrorCode PCFactorGetLevels_Factor(PC pc, PetscInt *levels)
 81: {
 82:   PC_Factor *ilu = (PC_Factor *)pc->data;

 84:   *levels = ilu->info.levels;
 85:   return 0;
 86: }

 88: PetscErrorCode PCFactorGetZeroPivot_Factor(PC pc, PetscReal *pivot)
 89: {
 90:   PC_Factor *ilu = (PC_Factor *)pc->data;

 92:   *pivot = ilu->info.zeropivot;
 93:   return 0;
 94: }

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

100:   *shift = ilu->info.shiftamount;
101:   return 0;
102: }

104: PetscErrorCode PCFactorGetShiftType_Factor(PC pc, MatFactorShiftType *type)
105: {
106:   PC_Factor *ilu = (PC_Factor *)pc->data;

108:   *type = (MatFactorShiftType)(int)ilu->info.shifttype;
109:   return 0;
110: }

112: PetscErrorCode PCFactorSetLevels_Factor(PC pc, PetscInt levels)
113: {
114:   PC_Factor *ilu = (PC_Factor *)pc->data;

116:   if (!pc->setupcalled) ilu->info.levels = levels;
117:   else if (ilu->info.levels != levels) {
118:     PetscUseTypeMethod(pc, reset); /* remove previous factored matrices */
119:     pc->setupcalled  = 0;          /* force a complete rebuild of preconditioner factored matrices */
120:     ilu->info.levels = levels;
122:   return 0;
123: }

125: PetscErrorCode PCFactorSetAllowDiagonalFill_Factor(PC pc, PetscBool flg)
126: {
127:   PC_Factor *dir = (PC_Factor *)pc->data;

129:   dir->info.diagonal_fill = (PetscReal)flg;
130:   return 0;
131: }

133: PetscErrorCode PCFactorGetAllowDiagonalFill_Factor(PC pc, PetscBool *flg)
134: {
135:   PC_Factor *dir = (PC_Factor *)pc->data;

137:   *flg = dir->info.diagonal_fill ? PETSC_TRUE : PETSC_FALSE;
138:   return 0;
139: }

141: PetscErrorCode PCFactorSetPivotInBlocks_Factor(PC pc, PetscBool pivot)
142: {
143:   PC_Factor *dir = (PC_Factor *)pc->data;

145:   dir->info.pivotinblocks = pivot ? 1.0 : 0.0;
146:   return 0;
147: }

149: PetscErrorCode PCFactorGetMatrix_Factor(PC pc, Mat *mat)
150: {
151:   PC_Factor *ilu = (PC_Factor *)pc->data;

154:   *mat = ilu->fact;
155:   return 0;
156: }

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

161: PetscErrorCode PCFactorSetMatSolverType_Factor(PC pc, MatSolverType stype)
162: {
163:   PC_Factor *lu = (PC_Factor *)pc->data;

165:   if (lu->fact && lu->fact->assembled) {
166:     MatSolverType ltype;
167:     PetscBool     flg;
168:     MatFactorGetSolverType(lu->fact, &ltype);
169:     PetscStrcmp(stype, ltype, &flg);
171:   }

173:   PetscFree(lu->solvertype);
174:   PetscStrallocpy(stype, &lu->solvertype);
175:   return 0;
176: }

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

182:   *stype = lu->solvertype;
183:   return 0;
184: }

186: PetscErrorCode PCFactorSetColumnPivot_Factor(PC pc, PetscReal dtcol)
187: {
188:   PC_Factor *dir = (PC_Factor *)pc->data;

191:   dir->info.dtcol = dtcol;
192:   return 0;
193: }

195: PetscErrorCode PCSetFromOptions_Factor(PC pc, PetscOptionItems *PetscOptionsObject)
196: {
197:   PC_Factor        *factor = (PC_Factor *)pc->data;
198:   PetscBool         flg, set;
199:   char              tname[256], solvertype[64];
200:   PetscFunctionList ordlist;
201:   PetscEnum         etmp;
202:   PetscBool         inplace;

204:   PCFactorGetUseInPlace(pc, &inplace);
205:   PetscOptionsBool("-pc_factor_in_place", "Form factored matrix in the same memory as the matrix", "PCFactorSetUseInPlace", inplace, &flg, &set);
206:   if (set) PCFactorSetUseInPlace(pc, flg);
207:   PetscOptionsReal("-pc_factor_fill", "Expected non-zeros in factored matrix", "PCFactorSetFill", ((PC_Factor *)factor)->info.fill, &((PC_Factor *)factor)->info.fill, NULL);

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

213:   PetscOptionsReal("-pc_factor_zeropivot", "Pivot is considered zero if less than", "PCFactorSetZeroPivot", ((PC_Factor *)factor)->info.zeropivot, &((PC_Factor *)factor)->info.zeropivot, NULL);
214:   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);

216:   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);
217:   if (set) PCFactorSetPivotInBlocks(pc, flg);

219:   PetscOptionsBool("-pc_factor_reuse_fill", "Use fill from previous factorization", "PCFactorSetReuseFill", PETSC_FALSE, &flg, &set);
220:   if (set) PCFactorSetReuseFill(pc, flg);
221:   PetscOptionsBool("-pc_factor_reuse_ordering", "Reuse ordering from previous factorization", "PCFactorSetReuseOrdering", PETSC_FALSE, &flg, &set);
222:   if (set) PCFactorSetReuseOrdering(pc, flg);

224:   PetscOptionsDeprecated("-pc_factor_mat_solver_package", "-pc_factor_mat_solver_type", "3.9", NULL);
225:   PetscOptionsString("-pc_factor_mat_solver_type", "Specific direct solver to use", "MatGetFactor", ((PC_Factor *)factor)->solvertype, solvertype, sizeof(solvertype), &flg);
226:   if (flg) PCFactorSetMatSolverType(pc, solvertype);
227:   PCFactorSetDefaultOrdering_Factor(pc);
228:   MatGetOrderingList(&ordlist);
229:   PetscOptionsFList("-pc_factor_mat_ordering_type", "Reordering to reduce nonzeros in factored matrix", "PCFactorSetMatOrderingType", ordlist, ((PC_Factor *)factor)->ordering, tname, sizeof(tname), &flg);
230:   if (flg) PCFactorSetMatOrderingType(pc, tname);
231:   return 0;
232: }

234: PetscErrorCode PCView_Factor(PC pc, PetscViewer viewer)
235: {
236:   PC_Factor        *factor = (PC_Factor *)pc->data;
237:   PetscBool         isstring, iascii, canuseordering;
238:   MatInfo           info;
239:   MatOrderingType   ordering;
240:   PetscViewerFormat format;

242:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring);
243:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
244:   if (iascii) {
245:     if (factor->inplace) {
246:       PetscViewerASCIIPrintf(viewer, "  in-place factorization\n");
247:     } else {
248:       PetscViewerASCIIPrintf(viewer, "  out-of-place factorization\n");
249:     }

251:     if (factor->reusefill) PetscViewerASCIIPrintf(viewer, "  Reusing fill from past factorization\n");
252:     if (factor->reuseordering) PetscViewerASCIIPrintf(viewer, "  Reusing reordering from past factorization\n");
253:     if (factor->factortype == MAT_FACTOR_ILU || factor->factortype == MAT_FACTOR_ICC) {
254:       if (factor->info.dt > 0) {
255:         PetscViewerASCIIPrintf(viewer, "  drop tolerance %g\n", (double)factor->info.dt);
256:         PetscViewerASCIIPrintf(viewer, "  max nonzeros per row %" PetscInt_FMT "\n", (PetscInt)factor->info.dtcount);
257:         PetscViewerASCIIPrintf(viewer, "  column permutation tolerance %g\n", (double)factor->info.dtcol);
258:       } else if (factor->info.levels == 1) {
259:         PetscViewerASCIIPrintf(viewer, "  %" PetscInt_FMT " level of fill\n", (PetscInt)factor->info.levels);
260:       } else {
261:         PetscViewerASCIIPrintf(viewer, "  %" PetscInt_FMT " levels of fill\n", (PetscInt)factor->info.levels);
262:       }
263:     }

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

270:     if (factor->fact) {
271:       MatFactorGetCanUseOrdering(factor->fact, &canuseordering);
272:       if (!canuseordering) ordering = MATORDERINGEXTERNAL;
273:       else ordering = factor->ordering;
274:       PetscViewerASCIIPrintf(viewer, "  matrix ordering: %s\n", ordering);
275:       if (!factor->fact->assembled) {
276:         PetscViewerASCIIPrintf(viewer, "  matrix solver type: %s\n", factor->fact->solvertype);
277:         PetscViewerASCIIPrintf(viewer, "  matrix not yet factored; no additional information available\n");
278:       } else {
279:         MatGetInfo(factor->fact, MAT_LOCAL, &info);
280:         PetscViewerASCIIPrintf(viewer, "  factor fill ratio given %g, needed %g\n", (double)info.fill_ratio_given, (double)info.fill_ratio_needed);
281:         PetscViewerASCIIPrintf(viewer, "    Factored matrix follows:\n");
282:         PetscViewerASCIIPushTab(viewer);
283:         PetscViewerASCIIPushTab(viewer);
284:         PetscViewerASCIIPushTab(viewer);
285:         PetscViewerGetFormat(viewer, &format);
286:         PetscViewerPushFormat(viewer, format != PETSC_VIEWER_ASCII_INFO_DETAIL ? PETSC_VIEWER_ASCII_INFO : PETSC_VIEWER_ASCII_INFO_DETAIL);
287:         MatView(factor->fact, viewer);
288:         PetscViewerPopFormat(viewer);
289:         PetscViewerASCIIPopTab(viewer);
290:         PetscViewerASCIIPopTab(viewer);
291:         PetscViewerASCIIPopTab(viewer);
292:       }
293:     }

295:   } else if (isstring) {
296:     MatFactorType t;
297:     MatGetFactorType(factor->fact, &t);
298:     if (t == MAT_FACTOR_ILU || t == MAT_FACTOR_ICC) PetscViewerStringSPrintf(viewer, " lvls=%" PetscInt_FMT ",order=%s", (PetscInt)factor->info.levels, factor->ordering);
299:   }
300:   return 0;
301: }