Actual source code: cholesky.c
petsc-3.10.5 2019-03-28
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 PCView_Cholesky(PC pc,PetscViewer viewer)
26: {
28: PetscBool iascii;
31: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
32: PCView_Factor(pc,viewer);
33: return(0);
34: }
36: static PetscErrorCode PCSetUp_Cholesky(PC pc)
37: {
38: PetscErrorCode ierr;
39: PetscBool flg;
40: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
41: MatSolverType stype;
42: MatFactorError err;
45: pc->failedreason = PC_NOERROR;
46: if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->hdr.actualfill;
48: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
49: if (dir->hdr.inplace) {
50: if (dir->row && dir->col && (dir->row != dir->col)) {
51: ISDestroy(&dir->row);
52: }
53: ISDestroy(&dir->col);
54: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
55: if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */
56: ISDestroy(&dir->col);
57: }
58: if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
59: MatCholeskyFactor(pc->pmat,dir->row,&((PC_Factor*)dir)->info);
60: MatFactorGetError(pc->pmat,&err);
61: if (err) { /* Factor() fails */
62: pc->failedreason = (PCFailedReason)err;
63: return(0);
64: }
66: ((PC_Factor*)dir)->fact = pc->pmat;
67: } else {
68: MatInfo info;
70: if (!pc->setupcalled) {
71: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
72: /* check if dir->row == dir->col */
73: ISEqual(dir->row,dir->col,&flg);
74: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"row and column permutations must equal");
75: ISDestroy(&dir->col); /* only pass one ordering into CholeskyFactor */
77: flg = PETSC_FALSE;
78: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
79: if (flg) {
80: PetscReal tol = 1.e-10;
81: PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
82: MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
83: }
84: if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
85: if (!((PC_Factor*)dir)->fact) {
86: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
87: }
88: MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
89: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
90: dir->hdr.actualfill = info.fill_ratio_needed;
91: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
92: } else if (pc->flag != SAME_NONZERO_PATTERN) {
93: if (!dir->hdr.reuseordering) {
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: if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
106: }
107: MatDestroy(&((PC_Factor*)dir)->fact);
108: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
109: MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
110: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
111: dir->hdr.actualfill = info.fill_ratio_needed;
112: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
113: } else {
114: MatFactorGetError(((PC_Factor*)dir)->fact,&err);
115: if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) {
116: MatFactorClearError(((PC_Factor*)dir)->fact);
117: pc->failedreason = PC_NOERROR;
118: }
119: }
120: MatFactorGetError(((PC_Factor*)dir)->fact,&err);
121: if (err) { /* FactorSymbolic() fails */
122: pc->failedreason = (PCFailedReason)err;
123: return(0);
124: }
126: MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
127: MatFactorGetError(((PC_Factor*)dir)->fact,&err);
128: if (err) { /* FactorNumeric() fails */
129: pc->failedreason = (PCFailedReason)err;
130: }
131: }
133: PCFactorGetMatSolverType(pc,&stype);
134: if (!stype) {
135: MatSolverType solverpackage;
136: MatFactorGetSolverType(((PC_Factor*)dir)->fact,&solverpackage);
137: PCFactorSetMatSolverType(pc,solverpackage);
138: }
139: return(0);
140: }
142: static PetscErrorCode PCReset_Cholesky(PC pc)
143: {
144: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
148: if (!dir->hdr.inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
149: ISDestroy(&dir->row);
150: ISDestroy(&dir->col);
151: return(0);
152: }
154: static PetscErrorCode PCDestroy_Cholesky(PC pc)
155: {
156: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
160: PCReset_Cholesky(pc);
161: PetscFree(((PC_Factor*)dir)->ordering);
162: PetscFree(((PC_Factor*)dir)->solvertype);
163: PetscFree(pc->data);
164: return(0);
165: }
167: static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y)
168: {
169: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
173: if (dir->hdr.inplace) {
174: MatSolve(pc->pmat,x,y);
175: } else {
176: MatSolve(((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;
187: if (dir->hdr.inplace) {
188: MatForwardSolve(pc->pmat,x,y);
189: } else {
190: MatForwardSolve(((PC_Factor*)dir)->fact,x,y);
191: }
192: return(0);
193: }
195: static PetscErrorCode PCApplySymmetricRight_Cholesky(PC pc,Vec x,Vec y)
196: {
197: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
201: if (dir->hdr.inplace) {
202: MatBackwardSolve(pc->pmat,x,y);
203: } else {
204: MatBackwardSolve(((PC_Factor*)dir)->fact,x,y);
205: }
206: return(0);
207: }
209: static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y)
210: {
211: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
215: if (dir->hdr.inplace) {
216: MatSolveTranspose(pc->pmat,x,y);
217: } else {
218: MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);
219: }
220: return(0);
221: }
223: /* -----------------------------------------------------------------------------------*/
225: /* -----------------------------------------------------------------------------------*/
227: /*@
228: PCFactorSetReuseOrdering - When similar matrices are factored, this
229: causes the ordering computed in the first factor to be used for all
230: following factors.
232: Logically Collective on PC
234: Input Parameters:
235: + pc - the preconditioner context
236: - flag - PETSC_TRUE to reuse else PETSC_FALSE
238: Options Database Key:
239: . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
241: Level: intermediate
243: .keywords: PC, levels, reordering, factorization, incomplete, LU
245: .seealso: PCFactorSetReuseFill()
246: @*/
247: PetscErrorCode PCFactorSetReuseOrdering(PC pc,PetscBool flag)
248: {
254: PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));
255: return(0);
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: Notes:
270: Not all options work for all matrix formats
272: Level: beginner
274: Concepts: Cholesky factorization, direct solver
276: Notes:
277: Usually this will compute an "exact" solution in one iteration and does
278: not need a Krylov method (i.e. you can use -ksp_type preonly, or
279: KSPSetType(ksp,KSPPREONLY) for the Krylov method
281: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
282: PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
283: PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount()
284: PCFactorSetUseInPlace(), PCFactorGetUseInPlace(), PCFactorSetMatOrderingType()
286: M*/
288: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
289: {
291: PC_Cholesky *dir;
294: PetscNewLog(pc,&dir);
295: pc->data = (void*)dir;
296: PCFactorInitialize(pc);
298: ((PC_Factor*)dir)->factortype = MAT_FACTOR_CHOLESKY;
299: ((PC_Factor*)dir)->info.fill = 5.0;
301: dir->col = 0;
302: dir->row = 0;
304: PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)dir)->ordering);
306: pc->ops->destroy = PCDestroy_Cholesky;
307: pc->ops->reset = PCReset_Cholesky;
308: pc->ops->apply = PCApply_Cholesky;
309: pc->ops->applysymmetricleft = PCApplySymmetricLeft_Cholesky;
310: pc->ops->applysymmetricright = PCApplySymmetricRight_Cholesky;
311: pc->ops->applytranspose = PCApplyTranspose_Cholesky;
312: pc->ops->setup = PCSetUp_Cholesky;
313: pc->ops->setfromoptions = PCSetFromOptions_Cholesky;
314: pc->ops->view = PCView_Cholesky;
315: pc->ops->applyrichardson = 0;
316: return(0);
317: }