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: 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;
92: PetscCall(MatDestroy(&((PC_Factor *)dir)->fact));
93: PetscCall(PCFactorSetUpMatSolverType(pc));
94: PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering));
95: if (canuseordering) {
96: PetscCall(ISDestroy(&dir->row));
97: PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
98: PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col));
99: PetscCall(ISDestroy(&dir->col)); /* only use dir->row ordering in CholeskyFactor */
101: flg = PETSC_FALSE;
102: PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &flg, NULL));
103: if (flg) {
104: PetscReal tol = 1.e-10;
105: PetscCall(PetscOptionsGetReal(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &tol, NULL));
106: PetscCall(MatReorderForNonzeroDiagonal(pc->pmat, tol, dir->row, dir->row));
107: }
108: }
109: }
110: PetscCall(MatCholeskyFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->row, &((PC_Factor *)dir)->info));
111: PetscCall(MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info));
112: dir->hdr.actualfill = info.fill_ratio_needed;
113: } else {
114: PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
115: if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) {
116: PetscCall(MatFactorClearError(((PC_Factor *)dir)->fact));
117: pc->failedreason = PC_NOERROR;
118: }
119: }
120: PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
121: if (err) { /* FactorSymbolic() fails */
122: pc->failedreason = (PCFailedReason)err;
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: PetscCall(MatCholeskyFactorNumeric(((PC_Factor *)dir)->fact, pc->pmat, &((PC_Factor *)dir)->info));
127: PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err));
128: if (err) { /* FactorNumeric() fails */
129: pc->failedreason = (PCFailedReason)err;
130: }
131: }
133: PetscCall(PCFactorGetMatSolverType(pc, &stype));
134: if (!stype) {
135: MatSolverType solverpackage;
136: PetscCall(MatFactorGetSolverType(((PC_Factor *)dir)->fact, &solverpackage));
137: PetscCall(PCFactorSetMatSolverType(pc, solverpackage));
138: }
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
142: static PetscErrorCode PCReset_Cholesky(PC pc)
143: {
144: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
146: PetscFunctionBegin;
147: if (!dir->hdr.inplace && ((PC_Factor *)dir)->fact) PetscCall(MatDestroy(&((PC_Factor *)dir)->fact));
148: PetscCall(ISDestroy(&dir->row));
149: PetscCall(ISDestroy(&dir->col));
150: PetscFunctionReturn(PETSC_SUCCESS);
151: }
153: static PetscErrorCode PCDestroy_Cholesky(PC pc)
154: {
155: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
157: PetscFunctionBegin;
158: PetscCall(PCReset_Cholesky(pc));
159: PetscCall(PetscFree(((PC_Factor *)dir)->ordering));
160: PetscCall(PetscFree(((PC_Factor *)dir)->solvertype));
161: PetscCall(PCFactorClearComposedFunctions(pc));
162: PetscCall(PetscFree(pc->data));
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: static PetscErrorCode PCApply_Cholesky(PC pc, Vec x, Vec y)
167: {
168: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
170: PetscFunctionBegin;
171: if (dir->hdr.inplace) {
172: PetscCall(MatSolve(pc->pmat, x, y));
173: } else {
174: PetscCall(MatSolve(((PC_Factor *)dir)->fact, x, y));
175: }
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: static PetscErrorCode PCMatApply_Cholesky(PC pc, Mat X, Mat Y)
180: {
181: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
183: PetscFunctionBegin;
184: if (dir->hdr.inplace) {
185: PetscCall(MatMatSolve(pc->pmat, X, Y));
186: } else {
187: PetscCall(MatMatSolve(((PC_Factor *)dir)->fact, X, Y));
188: }
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: static PetscErrorCode PCApplySymmetricLeft_Cholesky(PC pc, Vec x, Vec y)
193: {
194: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
196: PetscFunctionBegin;
197: if (dir->hdr.inplace) {
198: PetscCall(MatForwardSolve(pc->pmat, x, y));
199: } else {
200: PetscCall(MatForwardSolve(((PC_Factor *)dir)->fact, x, y));
201: }
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: static PetscErrorCode PCApplySymmetricRight_Cholesky(PC pc, Vec x, Vec y)
206: {
207: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
209: PetscFunctionBegin;
210: if (dir->hdr.inplace) {
211: PetscCall(MatBackwardSolve(pc->pmat, x, y));
212: } else {
213: PetscCall(MatBackwardSolve(((PC_Factor *)dir)->fact, x, y));
214: }
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: static PetscErrorCode PCApplyTranspose_Cholesky(PC pc, Vec x, Vec y)
219: {
220: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
222: PetscFunctionBegin;
223: if (dir->hdr.inplace) {
224: PetscCall(MatSolveTranspose(pc->pmat, x, y));
225: } else {
226: PetscCall(MatSolveTranspose(((PC_Factor *)dir)->fact, x, y));
227: }
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: /*@
232: PCFactorSetReuseOrdering - When similar matrices are factored, this
233: causes the ordering computed in the first factor to be used for all
234: following factors.
236: Logically Collective
238: Input Parameters:
239: + pc - the preconditioner context
240: - flag - `PETSC_TRUE` to reuse else `PETSC_FALSE`
242: Options Database Key:
243: . -pc_factor_reuse_ordering - Activate `PCFactorSetReuseOrdering()`
245: Level: intermediate
247: .seealso: [](ch_ksp), `PCLU`, `PCCHOLESKY`, `PCFactorSetReuseFill()`
248: @*/
249: PetscErrorCode PCFactorSetReuseOrdering(PC pc, PetscBool flag)
250: {
251: PetscFunctionBegin;
254: PetscTryMethod(pc, "PCFactorSetReuseOrdering_C", (PC, PetscBool), (pc, flag));
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: /*MC
259: PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner
261: Options Database Keys:
262: + -pc_factor_reuse_ordering - Activate `PCFactorSetReuseOrdering()`
263: . -pc_factor_mat_solver_type - Actives `PCFactorSetMatSolverType()` to choose the direct solver, like superlu
264: . -pc_factor_reuse_fill - Activates `PCFactorSetReuseFill()`
265: . -pc_factor_fill <fill> - Sets fill amount
266: . -pc_factor_in_place - Activates in-place factorization
267: - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
269: Level: beginner
271: Notes:
272: Not all options work for all matrix formats
274: Usually this will compute an "exact" solution in one iteration and does
275: not need a Krylov method (i.e. you can use -ksp_type preonly, or
276: `KSPSetType`(ksp,`KSPPREONLY`) for the Krylov method
278: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
279: `PCILU`, `PCLU`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`,
280: `PCFactorSetFill()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`
281: `PCFactorSetUseInPlace()`, `PCFactorGetUseInPlace()`, `PCFactorSetMatOrderingType()`, `PCFactorSetReuseOrdering()`
282: M*/
284: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
285: {
286: PC_Cholesky *dir;
288: PetscFunctionBegin;
289: PetscCall(PetscNew(&dir));
290: pc->data = (void *)dir;
291: PetscCall(PCFactorInitialize(pc, MAT_FACTOR_CHOLESKY));
293: ((PC_Factor *)dir)->info.fill = 5.0;
295: pc->ops->destroy = PCDestroy_Cholesky;
296: pc->ops->reset = PCReset_Cholesky;
297: pc->ops->apply = PCApply_Cholesky;
298: pc->ops->matapply = PCMatApply_Cholesky;
299: pc->ops->applysymmetricleft = PCApplySymmetricLeft_Cholesky;
300: pc->ops->applysymmetricright = PCApplySymmetricRight_Cholesky;
301: pc->ops->applytranspose = PCApplyTranspose_Cholesky;
302: pc->ops->setup = PCSetUp_Cholesky;
303: pc->ops->setfromoptions = PCSetFromOptions_Cholesky;
304: pc->ops->view = PCView_Factor;
305: pc->ops->applyrichardson = NULL;
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }