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: PetscOptionsHeadBegin(PetscOptionsObject, "Cholesky options");
17: PCSetFromOptions_Factor(pc, PetscOptionsObject);
18: PetscOptionsHeadEnd();
19: return 0;
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: pc->failedreason = PC_NOERROR;
31: if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor *)dir)->info.fill = dir->hdr.actualfill;
33: PCGetOptionsPrefix(pc, &prefix);
34: MatSetOptionsPrefixFactor(pc->pmat, prefix);
36: MatSetErrorIfFailure(pc->pmat, pc->erroriffailure);
37: if (dir->hdr.inplace) {
38: if (dir->row && dir->col && (dir->row != dir->col)) ISDestroy(&dir->row);
39: ISDestroy(&dir->col);
40: /* should only get reordering if the factor matrix uses it but cannot determine because MatGetFactor() not called */
41: PCFactorSetDefaultOrdering_Factor(pc);
42: MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col);
43: if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */
44: ISDestroy(&dir->col);
45: }
46: MatCholeskyFactor(pc->pmat, dir->row, &((PC_Factor *)dir)->info);
47: MatFactorGetError(pc->pmat, &err);
48: if (err) { /* Factor() fails */
49: pc->failedreason = (PCFailedReason)err;
50: return 0;
51: }
53: ((PC_Factor *)dir)->fact = pc->pmat;
54: } else {
55: MatInfo info;
57: if (!pc->setupcalled) {
58: PetscBool canuseordering;
59: if (!((PC_Factor *)dir)->fact) { MatGetFactor(pc->pmat, ((PC_Factor *)dir)->solvertype, MAT_FACTOR_CHOLESKY, &((PC_Factor *)dir)->fact); }
60: MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering);
61: if (canuseordering) {
62: PCFactorSetDefaultOrdering_Factor(pc);
63: MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col);
64: /* check if dir->row == dir->col */
65: if (dir->row) {
66: ISEqual(dir->row, dir->col, &flg);
68: }
69: ISDestroy(&dir->col); /* only pass one ordering into CholeskyFactor */
71: flg = PETSC_FALSE;
72: PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &flg, NULL);
73: if (flg) {
74: PetscReal tol = 1.e-10;
75: PetscOptionsGetReal(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &tol, NULL);
76: MatReorderForNonzeroDiagonal(pc->pmat, tol, dir->row, dir->row);
77: }
78: }
79: MatCholeskyFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->row, &((PC_Factor *)dir)->info);
80: MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info);
81: dir->hdr.actualfill = info.fill_ratio_needed;
82: } else if (pc->flag != SAME_NONZERO_PATTERN) {
83: if (!dir->hdr.reuseordering) {
84: PetscBool canuseordering;
85: MatDestroy(&((PC_Factor *)dir)->fact);
86: MatGetFactor(pc->pmat, ((PC_Factor *)dir)->solvertype, MAT_FACTOR_CHOLESKY, &((PC_Factor *)dir)->fact);
87: MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering);
88: if (canuseordering) {
89: ISDestroy(&dir->row);
90: PCFactorSetDefaultOrdering_Factor(pc);
91: MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col);
92: ISDestroy(&dir->col); /* only use dir->row ordering in CholeskyFactor */
94: flg = PETSC_FALSE;
95: PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &flg, NULL);
96: if (flg) {
97: PetscReal tol = 1.e-10;
98: PetscOptionsGetReal(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_factor_nonzeros_along_diagonal", &tol, NULL);
99: MatReorderForNonzeroDiagonal(pc->pmat, tol, dir->row, dir->row);
100: }
101: }
102: }
103: MatCholeskyFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->row, &((PC_Factor *)dir)->info);
104: MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info);
105: dir->hdr.actualfill = info.fill_ratio_needed;
106: } else {
107: MatFactorGetError(((PC_Factor *)dir)->fact, &err);
108: if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) {
109: MatFactorClearError(((PC_Factor *)dir)->fact);
110: pc->failedreason = PC_NOERROR;
111: }
112: }
113: MatFactorGetError(((PC_Factor *)dir)->fact, &err);
114: if (err) { /* FactorSymbolic() fails */
115: pc->failedreason = (PCFailedReason)err;
116: return 0;
117: }
119: MatCholeskyFactorNumeric(((PC_Factor *)dir)->fact, pc->pmat, &((PC_Factor *)dir)->info);
120: MatFactorGetError(((PC_Factor *)dir)->fact, &err);
121: if (err) { /* FactorNumeric() fails */
122: pc->failedreason = (PCFailedReason)err;
123: }
124: }
126: PCFactorGetMatSolverType(pc, &stype);
127: if (!stype) {
128: MatSolverType solverpackage;
129: MatFactorGetSolverType(((PC_Factor *)dir)->fact, &solverpackage);
130: PCFactorSetMatSolverType(pc, solverpackage);
131: }
132: return 0;
133: }
135: static PetscErrorCode PCReset_Cholesky(PC pc)
136: {
137: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
139: if (!dir->hdr.inplace && ((PC_Factor *)dir)->fact) MatDestroy(&((PC_Factor *)dir)->fact);
140: ISDestroy(&dir->row);
141: ISDestroy(&dir->col);
142: return 0;
143: }
145: static PetscErrorCode PCDestroy_Cholesky(PC pc)
146: {
147: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
149: PCReset_Cholesky(pc);
150: PetscFree(((PC_Factor *)dir)->ordering);
151: PetscFree(((PC_Factor *)dir)->solvertype);
152: PCFactorClearComposedFunctions(pc);
153: PetscFree(pc->data);
154: return 0;
155: }
157: static PetscErrorCode PCApply_Cholesky(PC pc, Vec x, Vec y)
158: {
159: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
161: if (dir->hdr.inplace) {
162: MatSolve(pc->pmat, x, y);
163: } else {
164: MatSolve(((PC_Factor *)dir)->fact, x, y);
165: }
166: return 0;
167: }
169: static PetscErrorCode PCMatApply_Cholesky(PC pc, Mat X, Mat Y)
170: {
171: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
173: if (dir->hdr.inplace) {
174: MatMatSolve(pc->pmat, X, Y);
175: } else {
176: MatMatSolve(((PC_Factor *)dir)->fact, X, Y);
177: }
178: return 0;
179: }
181: static PetscErrorCode PCApplySymmetricLeft_Cholesky(PC pc, Vec x, Vec y)
182: {
183: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
185: if (dir->hdr.inplace) {
186: MatForwardSolve(pc->pmat, x, y);
187: } else {
188: MatForwardSolve(((PC_Factor *)dir)->fact, x, y);
189: }
190: return 0;
191: }
193: static PetscErrorCode PCApplySymmetricRight_Cholesky(PC pc, Vec x, Vec y)
194: {
195: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
197: if (dir->hdr.inplace) {
198: MatBackwardSolve(pc->pmat, x, y);
199: } else {
200: MatBackwardSolve(((PC_Factor *)dir)->fact, x, y);
201: }
202: return 0;
203: }
205: static PetscErrorCode PCApplyTranspose_Cholesky(PC pc, Vec x, Vec y)
206: {
207: PC_Cholesky *dir = (PC_Cholesky *)pc->data;
209: if (dir->hdr.inplace) {
210: MatSolveTranspose(pc->pmat, x, y);
211: } else {
212: MatSolveTranspose(((PC_Factor *)dir)->fact, x, y);
213: }
214: return 0;
215: }
217: /*@
218: PCFactorSetReuseOrdering - When similar matrices are factored, this
219: causes the ordering computed in the first factor to be used for all
220: following factors.
222: Logically Collective
224: Input Parameters:
225: + pc - the preconditioner context
226: - flag - `PETSC_TRUE` to reuse else `PETSC_FALSE`
228: Options Database Key:
229: . -pc_factor_reuse_ordering - Activate `PCFactorSetReuseOrdering()`
231: Level: intermediate
233: .seealso: `PCLU`, `PCCHOLESKY`, `PCFactorSetReuseFill()`
234: @*/
235: PetscErrorCode PCFactorSetReuseOrdering(PC pc, PetscBool flag)
236: {
239: PetscTryMethod(pc, "PCFactorSetReuseOrdering_C", (PC, PetscBool), (pc, flag));
240: return 0;
241: }
243: /*MC
244: PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner
246: Options Database Keys:
247: + -pc_factor_reuse_ordering - Activate `PCFactorSetReuseOrdering()`
248: . -pc_factor_mat_solver_type - Actives `PCFactorSetMatSolverType()` to choose the direct solver, like superlu
249: . -pc_factor_reuse_fill - Activates `PCFactorSetReuseFill()`
250: . -pc_factor_fill <fill> - Sets fill amount
251: . -pc_factor_in_place - Activates in-place factorization
252: - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
254: Level: beginner
256: Notes:
257: Not all options work for all matrix formats
259: Usually this will compute an "exact" solution in one iteration and does
260: not need a Krylov method (i.e. you can use -ksp_type preonly, or
261: `KSPSetType`(ksp,`KSPPREONLY`) for the Krylov method
263: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
264: `PCILU`, `PCLU`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`,
265: `PCFactorSetFill()`, `PCFactorSetShiftNonzero()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`
266: `PCFactorSetUseInPlace()`, `PCFactorGetUseInPlace()`, `PCFactorSetMatOrderingType()`, `PCFactorSetReuseOrdering()`
267: M*/
269: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
270: {
271: PC_Cholesky *dir;
273: PetscNew(&dir);
274: pc->data = (void *)dir;
275: PCFactorInitialize(pc, MAT_FACTOR_CHOLESKY);
277: ((PC_Factor *)dir)->info.fill = 5.0;
279: pc->ops->destroy = PCDestroy_Cholesky;
280: pc->ops->reset = PCReset_Cholesky;
281: pc->ops->apply = PCApply_Cholesky;
282: pc->ops->matapply = PCMatApply_Cholesky;
283: pc->ops->applysymmetricleft = PCApplySymmetricLeft_Cholesky;
284: pc->ops->applysymmetricright = PCApplySymmetricRight_Cholesky;
285: pc->ops->applytranspose = PCApplyTranspose_Cholesky;
286: pc->ops->setup = PCSetUp_Cholesky;
287: pc->ops->setfromoptions = PCSetFromOptions_Cholesky;
288: pc->ops->view = PCView_Factor;
289: pc->ops->applyrichardson = NULL;
290: return 0;
291: }