Actual source code: cholesky.c
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>
8: typedef struct {
9: PC_Factor hdr;
10: IS row, col; /* index sets used for reordering */
11: } PC_Cholesky;
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: }
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;
30: PetscFunctionBegin;
31: pc->failedreason = PC_NOERROR;
32: if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor *)dir)->info.fill = dir->hdr.actualfill;
34: PetscCall(PCGetOptionsPrefix(pc, &prefix));
35: PetscCall(MatSetOptionsPrefixFactor(pc->pmat, prefix));
37: PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
38: if (dir->hdr.inplace) {
39: MatFactorType ftype;
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;
62: if (!pc->setupcalled) {
63: PetscBool canuseordering;
65: PetscCall(PCFactorSetUpMatSolverType(pc));
66: PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering));
67: if (canuseordering) {
68: PetscBool external;
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 */
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: if (!dir->hdr.reuseordering) {
95: PetscBool canuseordering;
97: PetscCall(MatDestroy(&((PC_Factor *)dir)->fact));
98: PetscCall(PCFactorSetUpMatSolverType(pc));
99: PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering));
100: if (canuseordering) {
101: PetscBool external;
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 */
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: }
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: }
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: }
152: static PetscErrorCode PCReset_Cholesky(PC pc)
153: {
154: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
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: }
163: static PetscErrorCode PCDestroy_Cholesky(PC pc)
164: {
165: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
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: }
176: static PetscErrorCode PCApply_Cholesky(PC pc, Vec x, Vec y)
177: {
178: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
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: }
189: static PetscErrorCode PCMatApply_Cholesky(PC pc, Mat X, Mat Y)
190: {
191: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
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: }
202: static PetscErrorCode PCApplySymmetricLeft_Cholesky(PC pc, Vec x, Vec y)
203: {
204: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
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: }
215: static PetscErrorCode PCApplySymmetricRight_Cholesky(PC pc, Vec x, Vec y)
216: {
217: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
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: }
228: static PetscErrorCode PCApplyTranspose_Cholesky(PC pc, Vec x, Vec y)
229: {
230: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
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: }
241: static PetscErrorCode PCMatApplyTranspose_Cholesky(PC pc, Mat X, Mat Y)
242: {
243: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
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: }
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.
259: Logically Collective
261: Input Parameters:
262: + pc - the preconditioner context
263: - flag - `PETSC_TRUE` to reuse else `PETSC_FALSE`
265: Options Database Key:
266: . -pc_factor_reuse_ordering - Activate `PCFactorSetReuseOrdering()`
268: Level: intermediate
270: .seealso: [](ch_ksp), `PCLU`, `PCCHOLESKY`, `PCFactorSetReuseFill()`
271: @*/
272: PetscErrorCode PCFactorSetReuseOrdering(PC pc, PetscBool flag)
273: {
274: PetscFunctionBegin;
277: PetscTryMethod(pc, "PCFactorSetReuseOrdering_C", (PC, PetscBool), (pc, flag));
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: /*MC
282: PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner
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
293: Level: beginner
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`.
299: Not all options work for all matrix formats
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
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*/
311: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
312: {
313: PC_Cholesky *dir;
315: PetscFunctionBegin;
316: PetscCall(PetscNew(&dir));
317: pc->data = (void *)dir;
318: PetscCall(PCFactorInitialize(pc, MAT_FACTOR_CHOLESKY));
320: ((PC_Factor *)dir)->info.fill = 5.0;
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: }