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