Actual source code: factimpl.c
1: #include <../src/ksp/pc/impls/factor/factor.h>
3: PetscErrorCode PCFactorSetUpMatSolverType_Factor(PC pc)
4: {
5: PC_Factor *icc = (PC_Factor *)pc->data;
7: PetscFunctionBegin;
8: 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()");
9: if (!icc->fact) PetscCall(MatGetFactor(pc->pmat, icc->solvertype, icc->factortype, &icc->fact));
10: PetscCheck(icc->fact, PetscObjectComm((PetscObject)pc->pmat), PETSC_ERR_SUP, "MatFactor type %s not supported by this matrix instance of type %s and solver type %s. Rerun with \"-info :mat | grep MatGetFactor_\" for additional information.",
11: MatFactorTypes[icc->factortype], ((PetscObject)pc->pmat)->type_name, icc->solvertype);
12: PetscFunctionReturn(PETSC_SUCCESS);
13: }
15: PetscErrorCode PCFactorSetZeroPivot_Factor(PC pc, PetscReal z)
16: {
17: PC_Factor *ilu = (PC_Factor *)pc->data;
19: PetscFunctionBegin;
20: ilu->info.zeropivot = z;
21: PetscFunctionReturn(PETSC_SUCCESS);
22: }
24: PetscErrorCode PCFactorSetShiftType_Factor(PC pc, MatFactorShiftType shifttype)
25: {
26: PC_Factor *dir = (PC_Factor *)pc->data;
28: PetscFunctionBegin;
29: if (shifttype == (MatFactorShiftType)PETSC_DECIDE) dir->info.shifttype = (PetscReal)MAT_SHIFT_NONE;
30: else {
31: dir->info.shifttype = (PetscReal)shifttype;
32: 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 */ }
33: }
34: PetscFunctionReturn(PETSC_SUCCESS);
35: }
37: PetscErrorCode PCFactorSetShiftAmount_Factor(PC pc, PetscReal shiftamount)
38: {
39: PC_Factor *dir = (PC_Factor *)pc->data;
41: PetscFunctionBegin;
42: if (shiftamount == (PetscReal)PETSC_DECIDE) dir->info.shiftamount = 100.0 * PETSC_MACHINE_EPSILON;
43: else dir->info.shiftamount = shiftamount;
44: PetscFunctionReturn(PETSC_SUCCESS);
45: }
47: PetscErrorCode PCFactorSetDropTolerance_Factor(PC pc, PetscReal dt, PetscReal dtcol, PetscInt dtcount)
48: {
49: PC_Factor *ilu = (PC_Factor *)pc->data;
51: PetscFunctionBegin;
52: PetscCheck(pc->setupcalled && (!ilu->info.usedt || ilu->info.dt != dt || ilu->info.dtcol != dtcol || ilu->info.dtcount != dtcount), PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot change tolerance after use");
53: ilu->info.usedt = PETSC_TRUE;
54: ilu->info.dt = dt;
55: ilu->info.dtcol = dtcol;
56: ilu->info.dtcount = dtcount;
57: ilu->info.fill = PETSC_DEFAULT;
58: PetscFunctionReturn(PETSC_SUCCESS);
59: }
61: PetscErrorCode PCFactorSetFill_Factor(PC pc, PetscReal fill)
62: {
63: PC_Factor *dir = (PC_Factor *)pc->data;
65: PetscFunctionBegin;
66: dir->info.fill = fill;
67: PetscFunctionReturn(PETSC_SUCCESS);
68: }
70: PetscErrorCode PCFactorSetMatOrderingType_Factor(PC pc, MatOrderingType ordering)
71: {
72: PC_Factor *dir = (PC_Factor *)pc->data;
73: PetscBool flg;
75: PetscFunctionBegin;
76: if (!pc->setupcalled) {
77: PetscCall(PetscFree(dir->ordering));
78: PetscCall(PetscStrallocpy(ordering, (char **)&dir->ordering));
79: } else {
80: PetscCall(PetscStrcmp(dir->ordering, ordering, &flg));
81: PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot change ordering after use");
82: }
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: PetscErrorCode PCFactorGetLevels_Factor(PC pc, PetscInt *levels)
87: {
88: PC_Factor *ilu = (PC_Factor *)pc->data;
90: PetscFunctionBegin;
91: *levels = ilu->info.levels;
92: PetscFunctionReturn(PETSC_SUCCESS);
93: }
95: PetscErrorCode PCFactorGetZeroPivot_Factor(PC pc, PetscReal *pivot)
96: {
97: PC_Factor *ilu = (PC_Factor *)pc->data;
99: PetscFunctionBegin;
100: *pivot = ilu->info.zeropivot;
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: PetscErrorCode PCFactorGetShiftAmount_Factor(PC pc, PetscReal *shift)
105: {
106: PC_Factor *ilu = (PC_Factor *)pc->data;
108: PetscFunctionBegin;
109: *shift = ilu->info.shiftamount;
110: PetscFunctionReturn(PETSC_SUCCESS);
111: }
113: PetscErrorCode PCFactorGetShiftType_Factor(PC pc, MatFactorShiftType *type)
114: {
115: PC_Factor *ilu = (PC_Factor *)pc->data;
117: PetscFunctionBegin;
118: *type = (MatFactorShiftType)(int)ilu->info.shifttype;
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: PetscErrorCode PCFactorSetLevels_Factor(PC pc, PetscInt levels)
123: {
124: PC_Factor *ilu = (PC_Factor *)pc->data;
126: PetscFunctionBegin;
127: if (!pc->setupcalled) ilu->info.levels = levels;
128: else if (ilu->info.levels != levels) {
129: PetscUseTypeMethod(pc, reset); /* remove previous factored matrices */
130: pc->setupcalled = 0; /* force a complete rebuild of preconditioner factored matrices */
131: ilu->info.levels = levels;
132: } else PetscCheck(!ilu->info.usedt, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot change levels after use with ILUdt");
133: PetscFunctionReturn(PETSC_SUCCESS);
134: }
136: PetscErrorCode PCFactorSetAllowDiagonalFill_Factor(PC pc, PetscBool flg)
137: {
138: PC_Factor *dir = (PC_Factor *)pc->data;
140: PetscFunctionBegin;
141: dir->info.diagonal_fill = (PetscReal)flg;
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
145: PetscErrorCode PCFactorGetAllowDiagonalFill_Factor(PC pc, PetscBool *flg)
146: {
147: PC_Factor *dir = (PC_Factor *)pc->data;
149: PetscFunctionBegin;
150: *flg = dir->info.diagonal_fill ? PETSC_TRUE : PETSC_FALSE;
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: PetscErrorCode PCFactorSetPivotInBlocks_Factor(PC pc, PetscBool pivot)
155: {
156: PC_Factor *dir = (PC_Factor *)pc->data;
158: PetscFunctionBegin;
159: dir->info.pivotinblocks = pivot ? 1.0 : 0.0;
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: PetscErrorCode PCFactorGetMatrix_Factor(PC pc, Mat *mat)
164: {
165: PC_Factor *ilu = (PC_Factor *)pc->data;
167: PetscFunctionBegin;
168: PetscCheck(ilu->fact, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Matrix not yet factored; call after KSPSetUp() or PCSetUp()");
169: *mat = ilu->fact;
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: /* allow access to preallocation information */
174: #include <petsc/private/matimpl.h>
176: PetscErrorCode PCFactorSetMatSolverType_Factor(PC pc, MatSolverType stype)
177: {
178: PC_Factor *lu = (PC_Factor *)pc->data;
180: PetscFunctionBegin;
181: if (lu->fact && lu->fact->assembled) {
182: MatSolverType ltype;
183: PetscBool flg;
184: PetscCall(MatFactorGetSolverType(lu->fact, <ype));
185: PetscCall(PetscStrcmp(stype, ltype, &flg));
186: 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);
187: }
189: PetscCall(PetscFree(lu->solvertype));
190: PetscCall(PetscStrallocpy(stype, &lu->solvertype));
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
194: PetscErrorCode PCFactorGetMatSolverType_Factor(PC pc, MatSolverType *stype)
195: {
196: PC_Factor *lu = (PC_Factor *)pc->data;
198: PetscFunctionBegin;
199: *stype = lu->solvertype;
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: PetscErrorCode PCFactorSetColumnPivot_Factor(PC pc, PetscReal dtcol)
204: {
205: PC_Factor *dir = (PC_Factor *)pc->data;
207: PetscFunctionBegin;
208: 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);
209: dir->info.dtcol = dtcol;
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: PetscErrorCode PCSetFromOptions_Factor(PC pc, PetscOptionItems *PetscOptionsObject)
214: {
215: PC_Factor *factor = (PC_Factor *)pc->data;
216: PetscBool flg, set;
217: char tname[256], solvertype[64];
218: PetscFunctionList ordlist;
219: PetscEnum etmp;
220: PetscBool inplace;
222: PetscFunctionBegin;
223: PetscCall(PCFactorGetUseInPlace(pc, &inplace));
224: PetscCall(PetscOptionsBool("-pc_factor_in_place", "Form factored matrix in the same memory as the matrix", "PCFactorSetUseInPlace", inplace, &flg, &set));
225: if (set) PetscCall(PCFactorSetUseInPlace(pc, flg));
226: PetscCall(PetscOptionsReal("-pc_factor_fill", "Expected non-zeros in factored matrix", "PCFactorSetFill", factor->info.fill, &factor->info.fill, NULL));
228: PetscCall(PetscOptionsEnum("-pc_factor_shift_type", "Type of shift to add to diagonal", "PCFactorSetShiftType", MatFactorShiftTypes, (PetscEnum)(int)factor->info.shifttype, &etmp, &flg));
229: if (flg) PetscCall(PCFactorSetShiftType(pc, (MatFactorShiftType)etmp));
230: PetscCall(PetscOptionsReal("-pc_factor_shift_amount", "Shift added to diagonal", "PCFactorSetShiftAmount", factor->info.shiftamount, &factor->info.shiftamount, NULL));
232: PetscCall(PetscOptionsReal("-pc_factor_zeropivot", "Pivot is considered zero if less than", "PCFactorSetZeroPivot", factor->info.zeropivot, &factor->info.zeropivot, NULL));
233: PetscCall(PetscOptionsReal("-pc_factor_column_pivot", "Column pivot tolerance (used only for some factorization)", "PCFactorSetColumnPivot", factor->info.dtcol, &factor->info.dtcol, &flg));
235: PetscCall(PetscOptionsBool("-pc_factor_pivot_in_blocks", "Pivot inside matrix dense blocks for BAIJ and SBAIJ", "PCFactorSetPivotInBlocks", factor->info.pivotinblocks ? PETSC_TRUE : PETSC_FALSE, &flg, &set));
236: if (set) PetscCall(PCFactorSetPivotInBlocks(pc, flg));
238: PetscCall(PetscOptionsBool("-pc_factor_reuse_fill", "Use fill from previous factorization", "PCFactorSetReuseFill", PETSC_FALSE, &flg, &set));
239: if (set) PetscCall(PCFactorSetReuseFill(pc, flg));
240: PetscCall(PetscOptionsBool("-pc_factor_reuse_ordering", "Reuse ordering from previous factorization", "PCFactorSetReuseOrdering", PETSC_FALSE, &flg, &set));
241: if (set) PetscCall(PCFactorSetReuseOrdering(pc, flg));
243: PetscCall(PetscOptionsDeprecated("-pc_factor_mat_solver_package", "-pc_factor_mat_solver_type", "3.9", NULL));
244: PetscCall(PetscOptionsString("-pc_factor_mat_solver_type", "Specific direct solver to use", "MatGetFactor", factor->solvertype, solvertype, sizeof(solvertype), &flg));
245: if (flg) PetscCall(PCFactorSetMatSolverType(pc, solvertype));
246: PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
247: PetscCall(MatGetOrderingList(&ordlist));
248: PetscCall(PetscOptionsFList("-pc_factor_mat_ordering_type", "Reordering to reduce nonzeros in factored matrix", "PCFactorSetMatOrderingType", ordlist, factor->ordering, tname, sizeof(tname), &flg));
249: if (flg) PetscCall(PCFactorSetMatOrderingType(pc, tname));
250: PetscFunctionReturn(PETSC_SUCCESS);
251: }
253: PetscErrorCode PCView_Factor(PC pc, PetscViewer viewer)
254: {
255: PC_Factor *factor = (PC_Factor *)pc->data;
256: PetscBool isstring, iascii, canuseordering;
257: MatInfo info;
258: MatOrderingType ordering;
259: PetscViewerFormat format;
261: PetscFunctionBegin;
262: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
263: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
264: if (iascii) {
265: if (factor->inplace) {
266: PetscCall(PetscViewerASCIIPrintf(viewer, " in-place factorization\n"));
267: } else {
268: PetscCall(PetscViewerASCIIPrintf(viewer, " out-of-place factorization\n"));
269: }
271: if (factor->reusefill) PetscCall(PetscViewerASCIIPrintf(viewer, " Reusing fill from past factorization\n"));
272: if (factor->reuseordering) PetscCall(PetscViewerASCIIPrintf(viewer, " Reusing reordering from past factorization\n"));
273: if (factor->factortype == MAT_FACTOR_ILU || factor->factortype == MAT_FACTOR_ICC) {
274: if (factor->info.dt > 0) {
275: PetscCall(PetscViewerASCIIPrintf(viewer, " drop tolerance %g\n", (double)factor->info.dt));
276: PetscCall(PetscViewerASCIIPrintf(viewer, " max nonzeros per row %" PetscInt_FMT "\n", (PetscInt)factor->info.dtcount));
277: PetscCall(PetscViewerASCIIPrintf(viewer, " column permutation tolerance %g\n", (double)factor->info.dtcol));
278: } else if (factor->info.levels == 1) {
279: PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " level of fill\n", (PetscInt)factor->info.levels));
280: } else {
281: PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " levels of fill\n", (PetscInt)factor->info.levels));
282: }
283: }
285: PetscCall(PetscViewerASCIIPrintf(viewer, " tolerance for zero pivot %g\n", (double)factor->info.zeropivot));
286: if (MatFactorShiftTypesDetail[(int)factor->info.shifttype]) { /* Only print when using a nontrivial shift */
287: PetscCall(PetscViewerASCIIPrintf(viewer, " using %s [%s]\n", MatFactorShiftTypesDetail[(int)factor->info.shifttype], MatFactorShiftTypes[(int)factor->info.shifttype]));
288: }
290: if (factor->fact) {
291: PetscCall(MatFactorGetCanUseOrdering(factor->fact, &canuseordering));
292: if (!canuseordering) ordering = MATORDERINGEXTERNAL;
293: else ordering = factor->ordering;
294: PetscCall(PetscViewerASCIIPrintf(viewer, " matrix ordering: %s\n", ordering));
295: if (!factor->fact->assembled) {
296: PetscCall(PetscViewerASCIIPrintf(viewer, " matrix solver type: %s\n", factor->fact->solvertype));
297: PetscCall(PetscViewerASCIIPrintf(viewer, " matrix not yet factored; no additional information available\n"));
298: } else {
299: PetscCall(MatGetInfo(factor->fact, MAT_LOCAL, &info));
300: PetscCall(PetscViewerASCIIPrintf(viewer, " factor fill ratio given %g, needed %g\n", info.fill_ratio_given, info.fill_ratio_needed));
301: PetscCall(PetscViewerASCIIPrintf(viewer, " Factored matrix follows:\n"));
302: PetscCall(PetscViewerASCIIPushTab(viewer));
303: PetscCall(PetscViewerASCIIPushTab(viewer));
304: PetscCall(PetscViewerASCIIPushTab(viewer));
305: PetscCall(PetscViewerGetFormat(viewer, &format));
306: PetscCall(PetscViewerPushFormat(viewer, format != PETSC_VIEWER_ASCII_INFO_DETAIL ? PETSC_VIEWER_ASCII_INFO : PETSC_VIEWER_ASCII_INFO_DETAIL));
307: PetscCall(MatView(factor->fact, viewer));
308: PetscCall(PetscViewerPopFormat(viewer));
309: PetscCall(PetscViewerASCIIPopTab(viewer));
310: PetscCall(PetscViewerASCIIPopTab(viewer));
311: PetscCall(PetscViewerASCIIPopTab(viewer));
312: }
313: }
315: } else if (isstring) {
316: MatFactorType t;
317: PetscCall(MatGetFactorType(factor->fact, &t));
318: if (t == MAT_FACTOR_ILU || t == MAT_FACTOR_ICC) PetscCall(PetscViewerStringSPrintf(viewer, " lvls=%" PetscInt_FMT ",order=%s", (PetscInt)factor->info.levels, factor->ordering));
319: }
320: PetscFunctionReturn(PETSC_SUCCESS);
321: }