1 #include <../src/ksp/pc/impls/factor/factor.h> /*I "petscpc.h" I*/
2
PCFactorSetUpMatSolverType_Factor(PC pc)3 PetscErrorCode PCFactorSetUpMatSolverType_Factor(PC pc)
4 {
5 PC_Factor *icc = (PC_Factor *)pc->data;
6
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 }
14
PCFactorSetZeroPivot_Factor(PC pc,PetscReal z)15 PetscErrorCode PCFactorSetZeroPivot_Factor(PC pc, PetscReal z)
16 {
17 PC_Factor *ilu = (PC_Factor *)pc->data;
18
19 PetscFunctionBegin;
20 ilu->info.zeropivot = z;
21 PetscFunctionReturn(PETSC_SUCCESS);
22 }
23
PCFactorSetShiftType_Factor(PC pc,MatFactorShiftType shifttype)24 PetscErrorCode PCFactorSetShiftType_Factor(PC pc, MatFactorShiftType shifttype)
25 {
26 PC_Factor *dir = (PC_Factor *)pc->data;
27
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 }
36
PCFactorSetShiftAmount_Factor(PC pc,PetscReal shiftamount)37 PetscErrorCode PCFactorSetShiftAmount_Factor(PC pc, PetscReal shiftamount)
38 {
39 PC_Factor *dir = (PC_Factor *)pc->data;
40
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 }
46
PCFactorSetDropTolerance_Factor(PC pc,PetscReal dt,PetscReal dtcol,PetscInt dtcount)47 PetscErrorCode PCFactorSetDropTolerance_Factor(PC pc, PetscReal dt, PetscReal dtcol, PetscInt dtcount)
48 {
49 PC_Factor *ilu = (PC_Factor *)pc->data;
50
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 }
60
PCFactorSetFill_Factor(PC pc,PetscReal fill)61 PetscErrorCode PCFactorSetFill_Factor(PC pc, PetscReal fill)
62 {
63 PC_Factor *dir = (PC_Factor *)pc->data;
64
65 PetscFunctionBegin;
66 dir->info.fill = fill;
67 PetscFunctionReturn(PETSC_SUCCESS);
68 }
69
PCFactorSetMatOrderingType_Factor(PC pc,MatOrderingType ordering)70 PetscErrorCode PCFactorSetMatOrderingType_Factor(PC pc, MatOrderingType ordering)
71 {
72 PC_Factor *dir = (PC_Factor *)pc->data;
73 PetscBool flg;
74
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 }
85
PCFactorGetLevels_Factor(PC pc,PetscInt * levels)86 PetscErrorCode PCFactorGetLevels_Factor(PC pc, PetscInt *levels)
87 {
88 PC_Factor *ilu = (PC_Factor *)pc->data;
89
90 PetscFunctionBegin;
91 *levels = ilu->info.levels;
92 PetscFunctionReturn(PETSC_SUCCESS);
93 }
94
PCFactorGetZeroPivot_Factor(PC pc,PetscReal * pivot)95 PetscErrorCode PCFactorGetZeroPivot_Factor(PC pc, PetscReal *pivot)
96 {
97 PC_Factor *ilu = (PC_Factor *)pc->data;
98
99 PetscFunctionBegin;
100 *pivot = ilu->info.zeropivot;
101 PetscFunctionReturn(PETSC_SUCCESS);
102 }
103
PCFactorGetShiftAmount_Factor(PC pc,PetscReal * shift)104 PetscErrorCode PCFactorGetShiftAmount_Factor(PC pc, PetscReal *shift)
105 {
106 PC_Factor *ilu = (PC_Factor *)pc->data;
107
108 PetscFunctionBegin;
109 *shift = ilu->info.shiftamount;
110 PetscFunctionReturn(PETSC_SUCCESS);
111 }
112
PCFactorGetShiftType_Factor(PC pc,MatFactorShiftType * type)113 PetscErrorCode PCFactorGetShiftType_Factor(PC pc, MatFactorShiftType *type)
114 {
115 PC_Factor *ilu = (PC_Factor *)pc->data;
116
117 PetscFunctionBegin;
118 *type = (MatFactorShiftType)(int)ilu->info.shifttype;
119 PetscFunctionReturn(PETSC_SUCCESS);
120 }
121
PCFactorSetLevels_Factor(PC pc,PetscInt levels)122 PetscErrorCode PCFactorSetLevels_Factor(PC pc, PetscInt levels)
123 {
124 PC_Factor *ilu = (PC_Factor *)pc->data;
125
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 = PETSC_FALSE; /* 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 }
135
PCFactorSetAllowDiagonalFill_Factor(PC pc,PetscBool flg)136 PetscErrorCode PCFactorSetAllowDiagonalFill_Factor(PC pc, PetscBool flg)
137 {
138 PC_Factor *dir = (PC_Factor *)pc->data;
139
140 PetscFunctionBegin;
141 dir->info.diagonal_fill = (PetscReal)flg;
142 PetscFunctionReturn(PETSC_SUCCESS);
143 }
144
PCFactorGetAllowDiagonalFill_Factor(PC pc,PetscBool * flg)145 PetscErrorCode PCFactorGetAllowDiagonalFill_Factor(PC pc, PetscBool *flg)
146 {
147 PC_Factor *dir = (PC_Factor *)pc->data;
148
149 PetscFunctionBegin;
150 *flg = dir->info.diagonal_fill ? PETSC_TRUE : PETSC_FALSE;
151 PetscFunctionReturn(PETSC_SUCCESS);
152 }
153
PCFactorSetPivotInBlocks_Factor(PC pc,PetscBool pivot)154 PetscErrorCode PCFactorSetPivotInBlocks_Factor(PC pc, PetscBool pivot)
155 {
156 PC_Factor *dir = (PC_Factor *)pc->data;
157
158 PetscFunctionBegin;
159 dir->info.pivotinblocks = pivot ? 1.0 : 0.0;
160 PetscFunctionReturn(PETSC_SUCCESS);
161 }
162
PCFactorGetMatrix_Factor(PC pc,Mat * mat)163 PetscErrorCode PCFactorGetMatrix_Factor(PC pc, Mat *mat)
164 {
165 PC_Factor *ilu = (PC_Factor *)pc->data;
166
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 }
172
173 /* allow access to preallocation information */
174 #include <petsc/private/matimpl.h>
175
PCFactorSetMatSolverType_Factor(PC pc,MatSolverType stype)176 PetscErrorCode PCFactorSetMatSolverType_Factor(PC pc, MatSolverType stype)
177 {
178 PC_Factor *lu = (PC_Factor *)pc->data;
179
180 PetscFunctionBegin;
181 if (lu->fact && lu->fact->assembled) {
182 MatSolverType ltype;
183 PetscBool flg;
184
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 }
189
190 PetscCall(PetscFree(lu->solvertype));
191 PetscCall(PetscStrallocpy(stype, &lu->solvertype));
192 PetscFunctionReturn(PETSC_SUCCESS);
193 }
194
PCFactorGetMatSolverType_Factor(PC pc,MatSolverType * stype)195 PetscErrorCode PCFactorGetMatSolverType_Factor(PC pc, MatSolverType *stype)
196 {
197 PC_Factor *lu = (PC_Factor *)pc->data;
198
199 PetscFunctionBegin;
200 *stype = lu->solvertype;
201 PetscFunctionReturn(PETSC_SUCCESS);
202 }
203
PCFactorSetColumnPivot_Factor(PC pc,PetscReal dtcol)204 PetscErrorCode PCFactorSetColumnPivot_Factor(PC pc, PetscReal dtcol)
205 {
206 PC_Factor *dir = (PC_Factor *)pc->data;
207
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 }
213
PCSetFromOptions_Factor(PC pc,PetscOptionItems PetscOptionsObject)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;
222
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", factor->info.fill, &factor->info.fill, NULL));
228
229 PetscCall(PetscOptionsEnum("-pc_factor_shift_type", "Type of shift to add to diagonal", "PCFactorSetShiftType", MatFactorShiftTypes, (PetscEnum)(int)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", factor->info.shiftamount, &factor->info.shiftamount, NULL));
232
233 PetscCall(PetscOptionsReal("-pc_factor_zeropivot", "Pivot is considered zero if less than", "PCFactorSetZeroPivot", factor->info.zeropivot, &factor->info.zeropivot, NULL));
234 PetscCall(PetscOptionsReal("-pc_factor_column_pivot", "Column pivot tolerance (used only for some factorization)", "PCFactorSetColumnPivot", factor->info.dtcol, &factor->info.dtcol, &flg));
235
236 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));
237 if (set) PetscCall(PCFactorSetPivotInBlocks(pc, flg));
238
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));
243
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", 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, factor->ordering, tname, sizeof(tname), &flg));
250 if (flg) PetscCall(PCFactorSetMatOrderingType(pc, tname));
251 PetscCall(PetscOptionsBool("-pc_factor_mat_factor_on_host", "Do mat factorization on host (with device matrix types)", "MatGetFactor", factor->info.factoronhost, &factor->info.factoronhost, NULL));
252 PetscCall(PetscOptionsBool("-pc_factor_mat_solve_on_host", "Do mat solve on host with the factor (with device matrix types)", "MatGetFactor", factor->info.solveonhost, &factor->info.solveonhost, NULL));
253 PetscFunctionReturn(PETSC_SUCCESS);
254 }
255
PCView_Factor(PC pc,PetscViewer viewer)256 PetscErrorCode PCView_Factor(PC pc, PetscViewer viewer)
257 {
258 PC_Factor *factor = (PC_Factor *)pc->data;
259 PetscBool isstring, isascii, canuseordering;
260 MatInfo info;
261 MatOrderingType ordering;
262 PetscViewerFormat format;
263
264 PetscFunctionBegin;
265 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
266 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
267 if (isascii) {
268 if (factor->inplace) {
269 PetscCall(PetscViewerASCIIPrintf(viewer, " in-place factorization\n"));
270 } else {
271 PetscCall(PetscViewerASCIIPrintf(viewer, " out-of-place factorization\n"));
272 }
273
274 if (factor->reusefill) PetscCall(PetscViewerASCIIPrintf(viewer, " Reusing fill from past factorization\n"));
275 if (factor->reuseordering) PetscCall(PetscViewerASCIIPrintf(viewer, " Reusing reordering from past factorization\n"));
276 if (factor->factortype == MAT_FACTOR_ILU || factor->factortype == MAT_FACTOR_ICC) {
277 if (factor->info.dt > 0) {
278 PetscCall(PetscViewerASCIIPrintf(viewer, " drop tolerance %g\n", (double)factor->info.dt));
279 PetscCall(PetscViewerASCIIPrintf(viewer, " max nonzeros per row %" PetscInt_FMT "\n", (PetscInt)factor->info.dtcount));
280 PetscCall(PetscViewerASCIIPrintf(viewer, " column permutation tolerance %g\n", (double)factor->info.dtcol));
281 } else if (factor->info.levels == 1) {
282 PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " level of fill\n", (PetscInt)factor->info.levels));
283 } else {
284 PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " levels of fill\n", (PetscInt)factor->info.levels));
285 }
286 }
287
288 PetscCall(PetscViewerASCIIPrintf(viewer, " tolerance for zero pivot %g\n", (double)factor->info.zeropivot));
289 if (MatFactorShiftTypesDetail[(int)factor->info.shifttype]) { /* Only print when using a nontrivial shift */
290 PetscCall(PetscViewerASCIIPrintf(viewer, " using %s [%s]\n", MatFactorShiftTypesDetail[(int)factor->info.shifttype], MatFactorShiftTypes[(int)factor->info.shifttype]));
291 }
292
293 if (factor->fact) {
294 PetscCall(MatFactorGetCanUseOrdering(factor->fact, &canuseordering));
295 if (!canuseordering) ordering = MATORDERINGEXTERNAL;
296 else ordering = factor->ordering;
297 PetscCall(PetscViewerASCIIPrintf(viewer, " matrix ordering: %s\n", ordering));
298 if (!factor->fact->assembled) {
299 PetscCall(PetscViewerASCIIPrintf(viewer, " matrix solver type: %s\n", factor->fact->solvertype));
300 PetscCall(PetscViewerASCIIPrintf(viewer, " matrix not yet factored; no additional information available\n"));
301 } else {
302 PetscCall(MatGetInfo(factor->fact, MAT_LOCAL, &info));
303 PetscCall(PetscViewerASCIIPrintf(viewer, " factor fill ratio given %g, needed %g\n", info.fill_ratio_given, info.fill_ratio_needed));
304 PetscCall(PetscViewerASCIIPrintf(viewer, " Factored matrix:\n"));
305 PetscCall(PetscViewerASCIIPushTab(viewer));
306 PetscCall(PetscViewerASCIIPushTab(viewer));
307 PetscCall(PetscViewerASCIIPushTab(viewer));
308 PetscCall(PetscViewerGetFormat(viewer, &format));
309 PetscCall(PetscViewerPushFormat(viewer, format != PETSC_VIEWER_ASCII_INFO_DETAIL ? PETSC_VIEWER_ASCII_INFO : PETSC_VIEWER_ASCII_INFO_DETAIL));
310 PetscCall(MatView(factor->fact, viewer));
311 PetscCall(PetscViewerPopFormat(viewer));
312 PetscCall(PetscViewerASCIIPopTab(viewer));
313 PetscCall(PetscViewerASCIIPopTab(viewer));
314 PetscCall(PetscViewerASCIIPopTab(viewer));
315 }
316 }
317
318 } else if (isstring) {
319 MatFactorType t;
320 PetscCall(MatGetFactorType(factor->fact, &t));
321 if (t == MAT_FACTOR_ILU || t == MAT_FACTOR_ICC) PetscCall(PetscViewerStringSPrintf(viewer, " lvls=%" PetscInt_FMT ",order=%s", (PetscInt)factor->info.levels, factor->ordering));
322 }
323 PetscFunctionReturn(PETSC_SUCCESS);
324 }
325