Actual source code: cholesky.c
petsc-3.9.4 2018-09-11
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 PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y)
182: {
183: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
187: if (dir->hdr.inplace) {
188: MatSolveTranspose(pc->pmat,x,y);
189: } else {
190: MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);
191: }
192: return(0);
193: }
195: /* -----------------------------------------------------------------------------------*/
197: /* -----------------------------------------------------------------------------------*/
199: /*@
200: PCFactorSetReuseOrdering - When similar matrices are factored, this
201: causes the ordering computed in the first factor to be used for all
202: following factors.
204: Logically Collective on PC
206: Input Parameters:
207: + pc - the preconditioner context
208: - flag - PETSC_TRUE to reuse else PETSC_FALSE
210: Options Database Key:
211: . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
213: Level: intermediate
215: .keywords: PC, levels, reordering, factorization, incomplete, LU
217: .seealso: PCFactorSetReuseFill()
218: @*/
219: PetscErrorCode PCFactorSetReuseOrdering(PC pc,PetscBool flag)
220: {
226: PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));
227: return(0);
228: }
230: /*MC
231: PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner
233: Options Database Keys:
234: + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
235: . -pc_factor_mat_solver_type - Actives PCFactorSetMatSolverType() to choose the direct solver, like superlu
236: . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
237: . -pc_factor_fill <fill> - Sets fill amount
238: . -pc_factor_in_place - Activates in-place factorization
239: - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
241: Notes: Not all options work for all matrix formats
243: Level: beginner
245: Concepts: Cholesky factorization, direct solver
247: Notes: Usually this will compute an "exact" solution in one iteration and does
248: not need a Krylov method (i.e. you can use -ksp_type preonly, or
249: KSPSetType(ksp,KSPPREONLY) for the Krylov method
251: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
252: PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
253: PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount()
254: PCFactorSetUseInPlace(), PCFactorGetUseInPlace(), PCFactorSetMatOrderingType()
256: M*/
258: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
259: {
261: PC_Cholesky *dir;
264: PetscNewLog(pc,&dir);
265: pc->data = (void*)dir;
266: PCFactorInitialize(pc);
268: ((PC_Factor*)dir)->factortype = MAT_FACTOR_CHOLESKY;
269: ((PC_Factor*)dir)->info.fill = 5.0;
271: dir->col = 0;
272: dir->row = 0;
274: PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)dir)->ordering);
276: pc->ops->destroy = PCDestroy_Cholesky;
277: pc->ops->reset = PCReset_Cholesky;
278: pc->ops->apply = PCApply_Cholesky;
279: pc->ops->applytranspose = PCApplyTranspose_Cholesky;
280: pc->ops->setup = PCSetUp_Cholesky;
281: pc->ops->setfromoptions = PCSetFromOptions_Cholesky;
282: pc->ops->view = PCView_Cholesky;
283: pc->ops->applyrichardson = 0;
284: return(0);
285: }