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: {
19: PetscOptionsHead(PetscOptionsObject,"Cholesky options");
20: PCSetFromOptions_Factor(PetscOptionsObject,pc);
21: PetscOptionsTail();
22: return(0);
23: }
25: static PetscErrorCode PCSetUp_Cholesky(PC pc)
26: {
27: PetscErrorCode ierr;
28: PetscBool flg;
29: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
30: MatSolverType stype;
31: MatFactorError err;
34: pc->failedreason = PC_NOERROR;
35: if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->hdr.actualfill;
37: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
38: if (dir->hdr.inplace) {
39: if (dir->row && dir->col && (dir->row != dir->col)) {
40: ISDestroy(&dir->row);
41: }
42: ISDestroy(&dir->col);
43: /* should only get reordering if the factor matrix uses it but cannot determine because MatGetFactor() not called */
44: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
45: if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */
46: ISDestroy(&dir->col);
47: }
48: if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
49: MatCholeskyFactor(pc->pmat,dir->row,&((PC_Factor*)dir)->info);
50: MatFactorGetError(pc->pmat,&err);
51: if (err) { /* Factor() fails */
52: pc->failedreason = (PCFailedReason)err;
53: return(0);
54: }
56: ((PC_Factor*)dir)->fact = pc->pmat;
57: } else {
58: MatInfo info;
60: if (!pc->setupcalled) {
61: PetscBool useordering;
62: if (!((PC_Factor*)dir)->fact) {
63: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
64: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
65: }
66: MatFactorGetUseOrdering(((PC_Factor*)dir)->fact,&useordering);
67: if (useordering) {
68: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
69: /* check if dir->row == dir->col */
70: ISEqual(dir->row,dir->col,&flg);
71: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"row and column permutations must equal");
72: ISDestroy(&dir->col); /* only pass one ordering into CholeskyFactor */
74: flg = PETSC_FALSE;
75: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
76: if (flg) {
77: PetscReal tol = 1.e-10;
78: PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
79: MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
80: }
81: PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);
82: }
83: MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
84: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
85: dir->hdr.actualfill = info.fill_ratio_needed;
86: } else if (pc->flag != SAME_NONZERO_PATTERN) {
87: if (!dir->hdr.reuseordering) {
88: PetscBool useordering;
89: MatDestroy(&((PC_Factor*)dir)->fact);
90: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
91: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
92: MatFactorGetUseOrdering(((PC_Factor*)dir)->fact,&useordering);
93: if (useordering) {
94: ISDestroy(&dir->row);
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;
147: if (!dir->hdr.inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
148: ISDestroy(&dir->row);
149: ISDestroy(&dir->col);
150: return(0);
151: }
153: static PetscErrorCode PCDestroy_Cholesky(PC pc)
154: {
155: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
159: PCReset_Cholesky(pc);
160: PetscFree(((PC_Factor*)dir)->ordering);
161: PetscFree(((PC_Factor*)dir)->solvertype);
162: PetscFree(pc->data);
163: return(0);
164: }
166: static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y)
167: {
168: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
172: if (dir->hdr.inplace) {
173: MatSolve(pc->pmat,x,y);
174: } else {
175: MatSolve(((PC_Factor*)dir)->fact,x,y);
176: }
177: return(0);
178: }
180: static PetscErrorCode PCMatApply_Cholesky(PC pc,Mat X,Mat Y)
181: {
182: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
186: if (dir->hdr.inplace) {
187: MatMatSolve(pc->pmat,X,Y);
188: } else {
189: MatMatSolve(((PC_Factor*)dir)->fact,X,Y);
190: }
191: return(0);
192: }
194: static PetscErrorCode PCApplySymmetricLeft_Cholesky(PC pc,Vec x,Vec y)
195: {
196: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
200: if (dir->hdr.inplace) {
201: MatForwardSolve(pc->pmat,x,y);
202: } else {
203: MatForwardSolve(((PC_Factor*)dir)->fact,x,y);
204: }
205: return(0);
206: }
208: static PetscErrorCode PCApplySymmetricRight_Cholesky(PC pc,Vec x,Vec y)
209: {
210: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
214: if (dir->hdr.inplace) {
215: MatBackwardSolve(pc->pmat,x,y);
216: } else {
217: MatBackwardSolve(((PC_Factor*)dir)->fact,x,y);
218: }
219: return(0);
220: }
222: static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y)
223: {
224: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
228: if (dir->hdr.inplace) {
229: MatSolveTranspose(pc->pmat,x,y);
230: } else {
231: MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);
232: }
233: return(0);
234: }
236: /* -----------------------------------------------------------------------------------*/
238: /* -----------------------------------------------------------------------------------*/
240: /*@
241: PCFactorSetReuseOrdering - When similar matrices are factored, this
242: causes the ordering computed in the first factor to be used for all
243: following factors.
245: Logically Collective on PC
247: Input Parameters:
248: + pc - the preconditioner context
249: - flag - PETSC_TRUE to reuse else PETSC_FALSE
251: Options Database Key:
252: . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
254: Level: intermediate
256: .seealso: PCFactorSetReuseFill()
257: @*/
258: PetscErrorCode PCFactorSetReuseOrdering(PC pc,PetscBool flag)
259: {
265: PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));
266: return(0);
267: }
269: /*MC
270: PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner
272: Options Database Keys:
273: + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
274: . -pc_factor_mat_solver_type - Actives PCFactorSetMatSolverType() to choose the direct solver, like superlu
275: . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
276: . -pc_factor_fill <fill> - Sets fill amount
277: . -pc_factor_in_place - Activates in-place factorization
278: - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
280: Notes:
281: Not all options work for all matrix formats
283: Level: beginner
285: Notes:
286: Usually this will compute an "exact" solution in one iteration and does
287: not need a Krylov method (i.e. you can use -ksp_type preonly, or
288: KSPSetType(ksp,KSPPREONLY) for the Krylov method
290: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
291: PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
292: PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount()
293: PCFactorSetUseInPlace(), PCFactorGetUseInPlace(), PCFactorSetMatOrderingType()
295: M*/
297: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
298: {
300: PC_Cholesky *dir;
303: PetscNewLog(pc,&dir);
304: pc->data = (void*)dir;
305: PCFactorInitialize(pc);
307: ((PC_Factor*)dir)->factortype = MAT_FACTOR_CHOLESKY;
308: ((PC_Factor*)dir)->info.fill = 5.0;
310: dir->col = NULL;
311: dir->row = NULL;
313: /* MATORDERINGNATURAL_OR_ND allows selecting type based on matrix type sbaij or aij */
314: PetscStrallocpy(MATORDERINGNATURAL_OR_ND,(char**)&((PC_Factor*)dir)->ordering);
316: pc->ops->destroy = PCDestroy_Cholesky;
317: pc->ops->reset = PCReset_Cholesky;
318: pc->ops->apply = PCApply_Cholesky;
319: pc->ops->matapply = PCMatApply_Cholesky;
320: pc->ops->applysymmetricleft = PCApplySymmetricLeft_Cholesky;
321: pc->ops->applysymmetricright = PCApplySymmetricRight_Cholesky;
322: pc->ops->applytranspose = PCApplyTranspose_Cholesky;
323: pc->ops->setup = PCSetUp_Cholesky;
324: pc->ops->setfromoptions = PCSetFromOptions_Cholesky;
325: pc->ops->view = PCView_Factor;
326: pc->ops->applyrichardson = NULL;
327: return(0);
328: }