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 #include <../src/ksp/pc/impls/factor/factor.h> /*I "petscpc.h" I*/
7
8 typedef struct {
9 PC_Factor hdr;
10 IS row, col; /* index sets used for reordering */
11 } PC_Cholesky;
12
PCSetFromOptions_Cholesky(PC pc,PetscOptionItems PetscOptionsObject)13 static PetscErrorCode PCSetFromOptions_Cholesky(PC pc, PetscOptionItems PetscOptionsObject)
14 {
15 PetscFunctionBegin;
16 PetscOptionsHeadBegin(PetscOptionsObject, "Cholesky options");
17 PetscCall(PCSetFromOptions_Factor(pc, PetscOptionsObject));
18 PetscOptionsHeadEnd();
19 PetscFunctionReturn(PETSC_SUCCESS);
20 }
21
PCSetUp_Cholesky(PC pc)22 static PetscErrorCode PCSetUp_Cholesky(PC pc)
23 {
24 PetscBool flg;
25 PC_Cholesky *dir = (PC_Cholesky *)pc->data;
26 MatSolverType stype;
27 MatFactorError err;
28 const char *prefix;
29
30 PetscFunctionBegin;
31 pc->failedreason = PC_NOERROR;
32 if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor *)dir)->info.fill = dir->hdr.actualfill;
33
34 PetscCall(PCGetOptionsPrefix(pc, &prefix));
35 PetscCall(MatSetOptionsPrefixFactor(pc->pmat, prefix));
36
37 PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
38 if (dir->hdr.inplace) {
39 MatFactorType ftype;
40
41 PetscCall(MatGetFactorType(pc->pmat, &ftype));
42 if (ftype == MAT_FACTOR_NONE) {
43 if (dir->row && dir->col && (dir->row != dir->col)) PetscCall(ISDestroy(&dir->row));
44 PetscCall(ISDestroy(&dir->col));
45 /* This should only get the ordering if needed, but since MatGetFactor() is not called we can't know if it is needed */
46 PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
47 PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col));
48 if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */
49 PetscCall(ISDestroy(&dir->col));
50 }
51 PetscCall(MatCholeskyFactor(pc->pmat, dir->row, &((PC_Factor *)dir)->info));
52 PetscCall(MatFactorGetError(pc->pmat, &err));
53 if (err) { /* Factor() fails */
54 pc->failedreason = (PCFailedReason)err;
55 PetscFunctionReturn(PETSC_SUCCESS);
56 }
57 }
58 ((PC_Factor *)dir)->fact = pc->pmat;
59 } else {
60 MatInfo info;
61
62 if (!pc->setupcalled) {
63 PetscBool canuseordering;
64
65 PetscCall(PCFactorSetUpMatSolverType(pc));
66 PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering));
67 if (canuseordering) {
68 PetscBool external;
69
70 PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
71 PetscCall(PetscStrcmp(((PC_Factor *)dir)->ordering, MATORDERINGEXTERNAL, &external));
72 if (!external) {
73 PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col));
74 /* check if dir->row == dir->col */
75 if (dir->row) {
76 PetscCall(ISEqual(dir->row, dir->col, &flg));
77 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "row and column permutations must be equal");
78 }
79 PetscCall(ISDestroy(&dir->col)); /* only pass one ordering into CholeskyFactor */
80
81 flg = PETSC_FALSE;
82 PetscCall(PetscOptionsHasName(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &flg));
83 if (flg) {
84 PetscReal tol = 1.e-10;
85 PetscCall(PetscOptionsGetReal(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &tol, NULL));
86 PetscCall(MatReorderForNonzeroDiagonal(pc->pmat, tol, dir->row, dir->row));
87 }
88 }
89 }
90 PetscCall(MatCholeskyFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->row, &((PC_Factor *)dir)->info));
91 PetscCall(MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info));
92 dir->hdr.actualfill = info.fill_ratio_needed;
93 } else if (pc->flag != SAME_NONZERO_PATTERN) {
94 PetscCall(MatDestroy(&((PC_Factor *)dir)->fact));
95 PetscCall(PCFactorSetUpMatSolverType(pc));
96 if (!dir->hdr.reuseordering) {
97 PetscBool canuseordering;
98
99 PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering));
100 if (canuseordering) {
101 PetscBool external;
102
103 PetscCall(ISDestroy(&dir->row));
104 PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
105 PetscCall(PetscStrcmp(((PC_Factor *)dir)->ordering, MATORDERINGEXTERNAL, &external));
106 if (!external) {
107 PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col));
108 PetscCall(ISDestroy(&dir->col)); /* only use dir->row ordering in CholeskyFactor */
109
110 flg = PETSC_FALSE;
111 PetscCall(PetscOptionsHasName(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &flg));
112 if (flg) {
113 PetscReal tol = 1.e-10;
114 PetscCall(PetscOptionsGetReal(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &tol, NULL));
115 PetscCall(MatReorderForNonzeroDiagonal(pc->pmat, tol, dir->row, dir->row));
116 }
117 }
118 }
119 }
120 PetscCall(MatCholeskyFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->row, &((PC_Factor *)dir)->info));
121 PetscCall(MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info));
122 dir->hdr.actualfill = info.fill_ratio_needed;
123 } else {
124 PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
125 if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) {
126 PetscCall(MatFactorClearError(((PC_Factor *)dir)->fact));
127 pc->failedreason = PC_NOERROR;
128 }
129 }
130 PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
131 if (err) { /* FactorSymbolic() fails */
132 pc->failedreason = (PCFailedReason)err;
133 PetscFunctionReturn(PETSC_SUCCESS);
134 }
135
136 PetscCall(MatCholeskyFactorNumeric(((PC_Factor *)dir)->fact, pc->pmat, &((PC_Factor *)dir)->info));
137 PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
138 if (err) { /* FactorNumeric() fails */
139 pc->failedreason = (PCFailedReason)err;
140 }
141 }
142
143 PetscCall(PCFactorGetMatSolverType(pc, &stype));
144 if (!stype) {
145 MatSolverType solverpackage;
146 PetscCall(MatFactorGetSolverType(((PC_Factor *)dir)->fact, &solverpackage));
147 PetscCall(PCFactorSetMatSolverType(pc, solverpackage));
148 }
149 PetscFunctionReturn(PETSC_SUCCESS);
150 }
151
PCReset_Cholesky(PC pc)152 static PetscErrorCode PCReset_Cholesky(PC pc)
153 {
154 PC_Cholesky *dir = (PC_Cholesky *)pc->data;
155
156 PetscFunctionBegin;
157 if (!dir->hdr.inplace && ((PC_Factor *)dir)->fact) PetscCall(MatDestroy(&((PC_Factor *)dir)->fact));
158 PetscCall(ISDestroy(&dir->row));
159 PetscCall(ISDestroy(&dir->col));
160 PetscFunctionReturn(PETSC_SUCCESS);
161 }
162
PCDestroy_Cholesky(PC pc)163 static PetscErrorCode PCDestroy_Cholesky(PC pc)
164 {
165 PC_Cholesky *dir = (PC_Cholesky *)pc->data;
166
167 PetscFunctionBegin;
168 PetscCall(PCReset_Cholesky(pc));
169 PetscCall(PetscFree(((PC_Factor *)dir)->ordering));
170 PetscCall(PetscFree(((PC_Factor *)dir)->solvertype));
171 PetscCall(PCFactorClearComposedFunctions(pc));
172 PetscCall(PetscFree(pc->data));
173 PetscFunctionReturn(PETSC_SUCCESS);
174 }
175
PCApply_Cholesky(PC pc,Vec x,Vec y)176 static PetscErrorCode PCApply_Cholesky(PC pc, Vec x, Vec y)
177 {
178 PC_Cholesky *dir = (PC_Cholesky *)pc->data;
179
180 PetscFunctionBegin;
181 if (dir->hdr.inplace) {
182 PetscCall(MatSolve(pc->pmat, x, y));
183 } else {
184 PetscCall(MatSolve(((PC_Factor *)dir)->fact, x, y));
185 }
186 PetscFunctionReturn(PETSC_SUCCESS);
187 }
188
PCMatApply_Cholesky(PC pc,Mat X,Mat Y)189 static PetscErrorCode PCMatApply_Cholesky(PC pc, Mat X, Mat Y)
190 {
191 PC_Cholesky *dir = (PC_Cholesky *)pc->data;
192
193 PetscFunctionBegin;
194 if (dir->hdr.inplace) {
195 PetscCall(MatMatSolve(pc->pmat, X, Y));
196 } else {
197 PetscCall(MatMatSolve(((PC_Factor *)dir)->fact, X, Y));
198 }
199 PetscFunctionReturn(PETSC_SUCCESS);
200 }
201
PCApplySymmetricLeft_Cholesky(PC pc,Vec x,Vec y)202 static PetscErrorCode PCApplySymmetricLeft_Cholesky(PC pc, Vec x, Vec y)
203 {
204 PC_Cholesky *dir = (PC_Cholesky *)pc->data;
205
206 PetscFunctionBegin;
207 if (dir->hdr.inplace) {
208 PetscCall(MatForwardSolve(pc->pmat, x, y));
209 } else {
210 PetscCall(MatForwardSolve(((PC_Factor *)dir)->fact, x, y));
211 }
212 PetscFunctionReturn(PETSC_SUCCESS);
213 }
214
PCApplySymmetricRight_Cholesky(PC pc,Vec x,Vec y)215 static PetscErrorCode PCApplySymmetricRight_Cholesky(PC pc, Vec x, Vec y)
216 {
217 PC_Cholesky *dir = (PC_Cholesky *)pc->data;
218
219 PetscFunctionBegin;
220 if (dir->hdr.inplace) {
221 PetscCall(MatBackwardSolve(pc->pmat, x, y));
222 } else {
223 PetscCall(MatBackwardSolve(((PC_Factor *)dir)->fact, x, y));
224 }
225 PetscFunctionReturn(PETSC_SUCCESS);
226 }
227
PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y)228 static PetscErrorCode PCApplyTranspose_Cholesky(PC pc, Vec x, Vec y)
229 {
230 PC_Cholesky *dir = (PC_Cholesky *)pc->data;
231
232 PetscFunctionBegin;
233 if (dir->hdr.inplace) {
234 PetscCall(MatSolveTranspose(pc->pmat, x, y));
235 } else {
236 PetscCall(MatSolveTranspose(((PC_Factor *)dir)->fact, x, y));
237 }
238 PetscFunctionReturn(PETSC_SUCCESS);
239 }
240
PCMatApplyTranspose_Cholesky(PC pc,Mat X,Mat Y)241 static PetscErrorCode PCMatApplyTranspose_Cholesky(PC pc, Mat X, Mat Y)
242 {
243 PC_Cholesky *dir = (PC_Cholesky *)pc->data;
244
245 PetscFunctionBegin;
246 if (dir->hdr.inplace) {
247 PetscCall(MatMatSolveTranspose(pc->pmat, X, Y));
248 } else {
249 PetscCall(MatMatSolveTranspose(((PC_Factor *)dir)->fact, X, Y));
250 }
251 PetscFunctionReturn(PETSC_SUCCESS);
252 }
253
254 /*@
255 PCFactorSetReuseOrdering - When similar matrices are factored, this
256 causes the ordering computed in the first factor to be used for all
257 following factors.
258
259 Logically Collective
260
261 Input Parameters:
262 + pc - the preconditioner context
263 - flag - `PETSC_TRUE` to reuse else `PETSC_FALSE`
264
265 Options Database Key:
266 . -pc_factor_reuse_ordering - Activate `PCFactorSetReuseOrdering()`
267
268 Level: intermediate
269
270 .seealso: [](ch_ksp), `PCLU`, `PCCHOLESKY`, `PCFactorSetReuseFill()`
271 @*/
PCFactorSetReuseOrdering(PC pc,PetscBool flag)272 PetscErrorCode PCFactorSetReuseOrdering(PC pc, PetscBool flag)
273 {
274 PetscFunctionBegin;
275 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
276 PetscValidLogicalCollectiveBool(pc, flag, 2);
277 PetscTryMethod(pc, "PCFactorSetReuseOrdering_C", (PC, PetscBool), (pc, flag));
278 PetscFunctionReturn(PETSC_SUCCESS);
279 }
280
281 /*MC
282 PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner
283
284 Options Database Keys:
285 + -pc_factor_reuse_ordering - Activate `PCFactorSetReuseOrdering()`
286 . -pc_factor_mat_solver_type - Actives `PCFactorSetMatSolverType()` to choose the direct solver, like superlu
287 . -pc_factor_reuse_fill - Activates `PCFactorSetReuseFill()`
288 . -pc_factor_fill <fill> - Sets the explected fill amount
289 . -pc_factor_in_place - Activates in-place factorization
290 - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine used to determine the order the rows are used in the factorization to reduce fill
291 and thus be more effective
292
293 Level: beginner
294
295 Notes:
296 The Cholesky factorization direct solver, `PCCHOLESKY` is only for symmetric positive-definite (SPD) matrices. For such
297 SPD matrices it is more efficient than using the LU factorization direct solver, `PCLU`.
298
299 Not all options work for all matrix formats
300
301 Usually this will compute an "exact" solution in one iteration and does
302 not need a Krylov method (i.e. you can use -ksp_type preonly, or
303 `KSPSetType`(ksp,`KSPPREONLY`) for the Krylov method
304
305 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
306 `PCILU`, `PCLU`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`,
307 `PCFactorSetFill()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`
308 `PCFactorSetUseInPlace()`, `PCFactorGetUseInPlace()`, `PCFactorSetMatOrderingType()`, `PCFactorSetReuseOrdering()`
309 M*/
310
PCCreate_Cholesky(PC pc)311 PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
312 {
313 PC_Cholesky *dir;
314
315 PetscFunctionBegin;
316 PetscCall(PetscNew(&dir));
317 pc->data = (void *)dir;
318 PetscCall(PCFactorInitialize(pc, MAT_FACTOR_CHOLESKY));
319
320 ((PC_Factor *)dir)->info.fill = 5.0;
321
322 pc->ops->destroy = PCDestroy_Cholesky;
323 pc->ops->reset = PCReset_Cholesky;
324 pc->ops->apply = PCApply_Cholesky;
325 pc->ops->matapply = PCMatApply_Cholesky;
326 pc->ops->applysymmetricleft = PCApplySymmetricLeft_Cholesky;
327 pc->ops->applysymmetricright = PCApplySymmetricRight_Cholesky;
328 pc->ops->applytranspose = PCApplyTranspose_Cholesky;
329 pc->ops->matapplytranspose = PCMatApplyTranspose_Cholesky;
330 pc->ops->setup = PCSetUp_Cholesky;
331 pc->ops->setfromoptions = PCSetFromOptions_Cholesky;
332 pc->ops->view = PCView_Factor;
333 pc->ops->applyrichardson = NULL;
334 PetscFunctionReturn(PETSC_SUCCESS);
335 }
336