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, <ype));
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: }