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: PCFactorSetDefaultOrdering_Factor(pc);
45: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
46: if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */
47: ISDestroy(&dir->col);
48: }
49: if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
50: MatCholeskyFactor(pc->pmat,dir->row,&((PC_Factor*)dir)->info);
51: MatFactorGetError(pc->pmat,&err);
52: if (err) { /* Factor() fails */
53: pc->failedreason = (PCFailedReason)err;
54: return(0);
55: }
57: ((PC_Factor*)dir)->fact = pc->pmat;
58: } else {
59: MatInfo info;
61: if (!pc->setupcalled) {
62: PetscBool canuseordering;
63: if (!((PC_Factor*)dir)->fact) {
64: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
65: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
66: }
67: MatFactorGetCanUseOrdering(((PC_Factor*)dir)->fact,&canuseordering);
68: if (canuseordering) {
69: PCFactorSetDefaultOrdering_Factor(pc);
70: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
71: /* check if dir->row == dir->col */
72: if (dir->row) {
73: ISEqual(dir->row,dir->col,&flg);
74: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"row and column permutations must be equal");
75: }
76: ISDestroy(&dir->col); /* only pass one ordering into CholeskyFactor */
78: flg = PETSC_FALSE;
79: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
80: if (flg) {
81: PetscReal tol = 1.e-10;
82: PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
83: MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
84: }
85: PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);
86: }
87: MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
88: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
89: dir->hdr.actualfill = info.fill_ratio_needed;
90: } else if (pc->flag != SAME_NONZERO_PATTERN) {
91: if (!dir->hdr.reuseordering) {
92: PetscBool canuseordering;
93: MatDestroy(&((PC_Factor*)dir)->fact);
94: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
95: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
96: MatFactorGetCanUseOrdering(((PC_Factor*)dir)->fact,&canuseordering);
97: if (canuseordering) {
98: ISDestroy(&dir->row);
99: PCFactorSetDefaultOrdering_Factor(pc);
100: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
101: ISDestroy(&dir->col); /* only use dir->row ordering in CholeskyFactor */
103: flg = PETSC_FALSE;
104: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
105: if (flg) {
106: PetscReal tol = 1.e-10;
107: PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
108: MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
109: }
110: PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);
111: }
112: }
113: MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
114: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
115: dir->hdr.actualfill = info.fill_ratio_needed;
116: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
117: } else {
118: MatFactorGetError(((PC_Factor*)dir)->fact,&err);
119: if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) {
120: MatFactorClearError(((PC_Factor*)dir)->fact);
121: pc->failedreason = PC_NOERROR;
122: }
123: }
124: MatFactorGetError(((PC_Factor*)dir)->fact,&err);
125: if (err) { /* FactorSymbolic() fails */
126: pc->failedreason = (PCFailedReason)err;
127: return(0);
128: }
130: MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
131: MatFactorGetError(((PC_Factor*)dir)->fact,&err);
132: if (err) { /* FactorNumeric() fails */
133: pc->failedreason = (PCFailedReason)err;
134: }
135: }
137: PCFactorGetMatSolverType(pc,&stype);
138: if (!stype) {
139: MatSolverType solverpackage;
140: MatFactorGetSolverType(((PC_Factor*)dir)->fact,&solverpackage);
141: PCFactorSetMatSolverType(pc,solverpackage);
142: }
143: return(0);
144: }
146: static PetscErrorCode PCReset_Cholesky(PC pc)
147: {
148: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
152: if (!dir->hdr.inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
153: ISDestroy(&dir->row);
154: ISDestroy(&dir->col);
155: return(0);
156: }
158: static PetscErrorCode PCDestroy_Cholesky(PC pc)
159: {
160: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
164: PCReset_Cholesky(pc);
165: PetscFree(((PC_Factor*)dir)->ordering);
166: PetscFree(((PC_Factor*)dir)->solvertype);
167: PetscFree(pc->data);
168: return(0);
169: }
171: static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y)
172: {
173: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
177: if (dir->hdr.inplace) {
178: MatSolve(pc->pmat,x,y);
179: } else {
180: MatSolve(((PC_Factor*)dir)->fact,x,y);
181: }
182: return(0);
183: }
185: static PetscErrorCode PCMatApply_Cholesky(PC pc,Mat X,Mat Y)
186: {
187: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
191: if (dir->hdr.inplace) {
192: MatMatSolve(pc->pmat,X,Y);
193: } else {
194: MatMatSolve(((PC_Factor*)dir)->fact,X,Y);
195: }
196: return(0);
197: }
199: static PetscErrorCode PCApplySymmetricLeft_Cholesky(PC pc,Vec x,Vec y)
200: {
201: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
205: if (dir->hdr.inplace) {
206: MatForwardSolve(pc->pmat,x,y);
207: } else {
208: MatForwardSolve(((PC_Factor*)dir)->fact,x,y);
209: }
210: return(0);
211: }
213: static PetscErrorCode PCApplySymmetricRight_Cholesky(PC pc,Vec x,Vec y)
214: {
215: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
219: if (dir->hdr.inplace) {
220: MatBackwardSolve(pc->pmat,x,y);
221: } else {
222: MatBackwardSolve(((PC_Factor*)dir)->fact,x,y);
223: }
224: return(0);
225: }
227: static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y)
228: {
229: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
233: if (dir->hdr.inplace) {
234: MatSolveTranspose(pc->pmat,x,y);
235: } else {
236: MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);
237: }
238: return(0);
239: }
241: /* -----------------------------------------------------------------------------------*/
243: /* -----------------------------------------------------------------------------------*/
245: /*@
246: PCFactorSetReuseOrdering - When similar matrices are factored, this
247: causes the ordering computed in the first factor to be used for all
248: following factors.
250: Logically Collective on PC
252: Input Parameters:
253: + pc - the preconditioner context
254: - flag - PETSC_TRUE to reuse else PETSC_FALSE
256: Options Database Key:
257: . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
259: Level: intermediate
261: .seealso: PCFactorSetReuseFill()
262: @*/
263: PetscErrorCode PCFactorSetReuseOrdering(PC pc,PetscBool flag)
264: {
270: PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));
271: return(0);
272: }
274: /*MC
275: PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner
277: Options Database Keys:
278: + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
279: . -pc_factor_mat_solver_type - Actives PCFactorSetMatSolverType() to choose the direct solver, like superlu
280: . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
281: . -pc_factor_fill <fill> - Sets fill amount
282: . -pc_factor_in_place - Activates in-place factorization
283: - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
285: Notes:
286: Not all options work for all matrix formats
288: Level: beginner
290: Notes:
291: Usually this will compute an "exact" solution in one iteration and does
292: not need a Krylov method (i.e. you can use -ksp_type preonly, or
293: KSPSetType(ksp,KSPPREONLY) for the Krylov method
295: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
296: PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
297: PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount()
298: PCFactorSetUseInPlace(), PCFactorGetUseInPlace(), PCFactorSetMatOrderingType()
300: M*/
302: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
303: {
305: PC_Cholesky *dir;
308: PetscNewLog(pc,&dir);
309: pc->data = (void*)dir;
310: PCFactorInitialize(pc,MAT_FACTOR_CHOLESKY);
312: ((PC_Factor*)dir)->info.fill = 5.0;
314: pc->ops->destroy = PCDestroy_Cholesky;
315: pc->ops->reset = PCReset_Cholesky;
316: pc->ops->apply = PCApply_Cholesky;
317: pc->ops->matapply = PCMatApply_Cholesky;
318: pc->ops->applysymmetricleft = PCApplySymmetricLeft_Cholesky;
319: pc->ops->applysymmetricright = PCApplySymmetricRight_Cholesky;
320: pc->ops->applytranspose = PCApplyTranspose_Cholesky;
321: pc->ops->setup = PCSetUp_Cholesky;
322: pc->ops->setfromoptions = PCSetFromOptions_Cholesky;
323: pc->ops->view = PCView_Factor;
324: pc->ops->applyrichardson = NULL;
325: return(0);
326: }