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