1 /*
2 Defines a direct factorization preconditioner for any Mat implementation
3 Note: this need not be considered a preconditioner since it supplies
4 a direct solver.
5 */
6
7 #include <../src/ksp/pc/impls/factor/lu/lu.h> /*I "petscpc.h" I*/
8
PCFactorReorderForNonzeroDiagonal_LU(PC pc,PetscReal z)9 static PetscErrorCode PCFactorReorderForNonzeroDiagonal_LU(PC pc, PetscReal z)
10 {
11 PC_LU *lu = (PC_LU *)pc->data;
12
13 PetscFunctionBegin;
14 lu->nonzerosalongdiagonal = PETSC_TRUE;
15 if (z == (PetscReal)PETSC_DECIDE) lu->nonzerosalongdiagonaltol = 1.e-10;
16 else lu->nonzerosalongdiagonaltol = z;
17 PetscFunctionReturn(PETSC_SUCCESS);
18 }
19
PCSetFromOptions_LU(PC pc,PetscOptionItems PetscOptionsObject)20 static PetscErrorCode PCSetFromOptions_LU(PC pc, PetscOptionItems PetscOptionsObject)
21 {
22 PC_LU *lu = (PC_LU *)pc->data;
23 PetscBool flg = PETSC_FALSE;
24 PetscReal tol;
25
26 PetscFunctionBegin;
27 PetscOptionsHeadBegin(PetscOptionsObject, "LU options");
28 PetscCall(PCSetFromOptions_Factor(pc, PetscOptionsObject));
29
30 PetscCall(PetscOptionsName("-pc_factor_nonzeros_along_diagonal", "Reorder to remove zeros from diagonal", "PCFactorReorderForNonzeroDiagonal", &flg));
31 if (flg) {
32 tol = PETSC_DECIDE;
33 PetscCall(PetscOptionsReal("-pc_factor_nonzeros_along_diagonal", "Reorder to remove zeros from diagonal", "PCFactorReorderForNonzeroDiagonal", lu->nonzerosalongdiagonaltol, &tol, NULL));
34 PetscCall(PCFactorReorderForNonzeroDiagonal(pc, tol));
35 }
36 PetscOptionsHeadEnd();
37 PetscFunctionReturn(PETSC_SUCCESS);
38 }
39
PCSetUp_LU(PC pc)40 static PetscErrorCode PCSetUp_LU(PC pc)
41 {
42 PC_LU *dir = (PC_LU *)pc->data;
43 MatSolverType stype;
44 MatFactorError err;
45 const char *prefix;
46
47 PetscFunctionBegin;
48 pc->failedreason = PC_NOERROR;
49 if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor *)dir)->info.fill = dir->hdr.actualfill;
50
51 PetscCall(PCGetOptionsPrefix(pc, &prefix));
52 PetscCall(MatSetOptionsPrefixFactor(pc->pmat, prefix));
53
54 PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
55 if (dir->hdr.inplace) {
56 MatFactorType ftype;
57
58 PetscCall(MatGetFactorType(pc->pmat, &ftype));
59 if (ftype == MAT_FACTOR_NONE) {
60 if (dir->row && dir->col && dir->row != dir->col) PetscCall(ISDestroy(&dir->row));
61 PetscCall(ISDestroy(&dir->col));
62 /* This should only get the ordering if needed, but since MatGetFactor() is not called we can't know if it is needed */
63 PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
64 PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col));
65 PetscCall(MatLUFactor(pc->pmat, dir->row, dir->col, &((PC_Factor *)dir)->info));
66 PetscCall(MatFactorGetError(pc->pmat, &err));
67 if (err) { /* Factor() fails */
68 pc->failedreason = (PCFailedReason)err;
69 PetscFunctionReturn(PETSC_SUCCESS);
70 }
71 }
72 ((PC_Factor *)dir)->fact = pc->pmat;
73 } else {
74 MatInfo info;
75
76 if (!pc->setupcalled) {
77 PetscBool canuseordering;
78
79 PetscCall(PCFactorSetUpMatSolverType(pc));
80 PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering));
81 if (canuseordering) {
82 PetscBool external;
83
84 PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
85 PetscCall(PetscStrcmp(((PC_Factor *)dir)->ordering, MATORDERINGEXTERNAL, &external));
86 if (!external) {
87 PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col));
88 if (dir->nonzerosalongdiagonal) PetscCall(MatReorderForNonzeroDiagonal(pc->pmat, dir->nonzerosalongdiagonaltol, dir->row, dir->col));
89 }
90 }
91 PetscCall(MatLUFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->row, dir->col, &((PC_Factor *)dir)->info));
92 PetscCall(MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info));
93 dir->hdr.actualfill = info.fill_ratio_needed;
94 } else if (pc->flag != SAME_NONZERO_PATTERN) {
95 PetscBool canuseordering;
96
97 PetscCall(MatDestroy(&((PC_Factor *)dir)->fact));
98 PetscCall(PCFactorSetUpMatSolverType(pc));
99 if (!dir->hdr.reuseordering) {
100 PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering));
101 if (canuseordering) {
102 PetscBool external;
103
104 if (dir->row && dir->col && dir->row != dir->col) PetscCall(ISDestroy(&dir->row));
105 PetscCall(ISDestroy(&dir->col));
106 PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
107 PetscCall(PetscStrcmp(((PC_Factor *)dir)->ordering, MATORDERINGEXTERNAL, &external));
108 if (!external) {
109 PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col));
110 if (dir->nonzerosalongdiagonal) PetscCall(MatReorderForNonzeroDiagonal(pc->pmat, dir->nonzerosalongdiagonaltol, dir->row, dir->col));
111 }
112 }
113 }
114 PetscCall(MatLUFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->row, dir->col, &((PC_Factor *)dir)->info));
115 PetscCall(MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info));
116 dir->hdr.actualfill = info.fill_ratio_needed;
117 } else {
118 PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
119 if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) {
120 PetscCall(MatFactorClearError(((PC_Factor *)dir)->fact));
121 pc->failedreason = PC_NOERROR;
122 }
123 }
124 PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
125 if (err) { /* FactorSymbolic() fails */
126 pc->failedreason = (PCFailedReason)err;
127 PetscFunctionReturn(PETSC_SUCCESS);
128 }
129
130 PetscCall(MatLUFactorNumeric(((PC_Factor *)dir)->fact, pc->pmat, &((PC_Factor *)dir)->info));
131 PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
132 if (err) { /* FactorNumeric() fails */
133 pc->failedreason = (PCFailedReason)err;
134 }
135 }
136
137 PetscCall(PCFactorGetMatSolverType(pc, &stype));
138 if (!stype) {
139 MatSolverType solverpackage;
140 PetscCall(MatFactorGetSolverType(((PC_Factor *)dir)->fact, &solverpackage));
141 PetscCall(PCFactorSetMatSolverType(pc, solverpackage));
142 }
143 PetscFunctionReturn(PETSC_SUCCESS);
144 }
145
PCReset_LU(PC pc)146 static PetscErrorCode PCReset_LU(PC pc)
147 {
148 PC_LU *dir = (PC_LU *)pc->data;
149
150 PetscFunctionBegin;
151 if (!dir->hdr.inplace && ((PC_Factor *)dir)->fact) PetscCall(MatDestroy(&((PC_Factor *)dir)->fact));
152 if (dir->row && dir->col && dir->row != dir->col) PetscCall(ISDestroy(&dir->row));
153 PetscCall(ISDestroy(&dir->col));
154 PetscFunctionReturn(PETSC_SUCCESS);
155 }
156
PCDestroy_LU(PC pc)157 static PetscErrorCode PCDestroy_LU(PC pc)
158 {
159 PC_LU *dir = (PC_LU *)pc->data;
160
161 PetscFunctionBegin;
162 PetscCall(PCReset_LU(pc));
163 PetscCall(PetscFree(((PC_Factor *)dir)->ordering));
164 PetscCall(PetscFree(((PC_Factor *)dir)->solvertype));
165 PetscCall(PCFactorClearComposedFunctions(pc));
166 PetscCall(PetscFree(pc->data));
167 PetscFunctionReturn(PETSC_SUCCESS);
168 }
169
PCApply_LU(PC pc,Vec x,Vec y)170 static PetscErrorCode PCApply_LU(PC pc, Vec x, Vec y)
171 {
172 PC_LU *dir = (PC_LU *)pc->data;
173
174 PetscFunctionBegin;
175 if (dir->hdr.inplace) {
176 PetscCall(MatSolve(pc->pmat, x, y));
177 } else {
178 PetscCall(MatSolve(((PC_Factor *)dir)->fact, x, y));
179 }
180 PetscFunctionReturn(PETSC_SUCCESS);
181 }
182
PCMatApply_LU(PC pc,Mat X,Mat Y)183 static PetscErrorCode PCMatApply_LU(PC pc, Mat X, Mat Y)
184 {
185 PC_LU *dir = (PC_LU *)pc->data;
186
187 PetscFunctionBegin;
188 if (dir->hdr.inplace) {
189 PetscCall(MatMatSolve(pc->pmat, X, Y));
190 } else {
191 PetscCall(MatMatSolve(((PC_Factor *)dir)->fact, X, Y));
192 }
193 PetscFunctionReturn(PETSC_SUCCESS);
194 }
195
PCApplyTranspose_LU(PC pc,Vec x,Vec y)196 static PetscErrorCode PCApplyTranspose_LU(PC pc, Vec x, Vec y)
197 {
198 PC_LU *dir = (PC_LU *)pc->data;
199
200 PetscFunctionBegin;
201 if (dir->hdr.inplace) {
202 PetscCall(MatSolveTranspose(pc->pmat, x, y));
203 } else {
204 PetscCall(MatSolveTranspose(((PC_Factor *)dir)->fact, x, y));
205 }
206 PetscFunctionReturn(PETSC_SUCCESS);
207 }
208
PCMatApplyTranspose_LU(PC pc,Mat X,Mat Y)209 static PetscErrorCode PCMatApplyTranspose_LU(PC pc, Mat X, Mat Y)
210 {
211 PC_LU *dir = (PC_LU *)pc->data;
212
213 PetscFunctionBegin;
214 if (dir->hdr.inplace) {
215 PetscCall(MatMatSolveTranspose(pc->pmat, X, Y));
216 } else {
217 PetscCall(MatMatSolveTranspose(((PC_Factor *)dir)->fact, X, Y));
218 }
219 PetscFunctionReturn(PETSC_SUCCESS);
220 }
221
222 /*MC
223 PCLU - Uses a direct solver, based on LU factorization, as a preconditioner
224
225 Options Database Keys:
226 + -pc_factor_reuse_ordering - Activate `PCFactorSetReuseOrdering()`
227 . -pc_factor_mat_solver_type - Actives `PCFactorSetMatSolverType()` to choose the direct solver, like superlu
228 . -pc_factor_reuse_fill - Activates `PCFactorSetReuseFill()`
229 . -pc_factor_fill <fill> - Sets fill amount
230 . -pc_factor_in_place - Activates in-place factorization
231 . -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
232 . -pc_factor_pivot_in_blocks <true,false> - allow pivoting within the small blocks during factorization (may increase
233 stability of factorization.
234 . -pc_factor_shift_type <shifttype> - Sets shift type or -1 for the default; use '-help' for a list of available types
235 . -pc_factor_shift_amount <shiftamount> - Sets shift amount or -1 for the default
236 . -pc_factor_nonzeros_along_diagonal - permutes the rows and columns to try to put nonzero value along the diagonal.
237 . -pc_factor_mat_solver_type <packagename> - use an external package for the solve, see `MatSolverType` for possibilities
238 - -mat_solvertype_optionname - options for a specific solver package, for example -mat_mumps_cntl_1
239
240 Level: beginner
241
242 Notes:
243 Not all options work for all matrix formats
244
245 Run with `-help` to see additional options for particular matrix formats or factorization algorithms
246
247 The Cholesky factorization direct solver, `PCCHOLESKY` will be more efficient than `PCLU` for symmetric positive-definite (SPD) matrices
248
249 Usually this will compute an "exact" solution in one iteration and does
250 not need a Krylov method (i.e. you can use -ksp_type preonly, or
251 `KSPSetType`(ksp,`KSPPREONLY`) for the Krylov method.
252
253 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `MatSolverType`, `MatGetFactor()`, `PCQR`, `PCSVD`,
254 `PCILU`, `PCCHOLESKY`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`,
255 `PCFactorSetFill()`, `PCFactorSetUseInPlace()`, `PCFactorSetMatOrderingType()`, `PCFactorSetColumnPivot()`,
256 `PCFactorSetPivotInBlocks()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`
257 `PCFactorReorderForNonzeroDiagonal()`
258 M*/
259
PCCreate_LU(PC pc)260 PETSC_EXTERN PetscErrorCode PCCreate_LU(PC pc)
261 {
262 PC_LU *dir;
263
264 PetscFunctionBegin;
265 PetscCall(PetscNew(&dir));
266 pc->data = (void *)dir;
267 PetscCall(PCFactorInitialize(pc, MAT_FACTOR_LU));
268 dir->nonzerosalongdiagonal = PETSC_FALSE;
269
270 ((PC_Factor *)dir)->info.fill = 5.0;
271 ((PC_Factor *)dir)->info.dtcol = 1.e-6; /* default to pivoting; this is only thing PETSc LU supports */
272 ((PC_Factor *)dir)->info.shifttype = (PetscReal)MAT_SHIFT_NONE;
273 dir->col = NULL;
274 dir->row = NULL;
275
276 pc->ops->reset = PCReset_LU;
277 pc->ops->destroy = PCDestroy_LU;
278 pc->ops->apply = PCApply_LU;
279 pc->ops->matapply = PCMatApply_LU;
280 pc->ops->applytranspose = PCApplyTranspose_LU;
281 pc->ops->matapplytranspose = PCMatApplyTranspose_LU;
282 pc->ops->setup = PCSetUp_LU;
283 pc->ops->setfromoptions = PCSetFromOptions_LU;
284 pc->ops->view = PCView_Factor;
285 pc->ops->applyrichardson = NULL;
286 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorReorderForNonzeroDiagonal_C", PCFactorReorderForNonzeroDiagonal_LU));
287 PetscFunctionReturn(PETSC_SUCCESS);
288 }
289