Actual source code: cholesky.c
petsc-3.13.6 2020-09-29
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: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
44: if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */
45: ISDestroy(&dir->col);
46: }
47: if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
48: MatCholeskyFactor(pc->pmat,dir->row,&((PC_Factor*)dir)->info);
49: MatFactorGetError(pc->pmat,&err);
50: if (err) { /* Factor() fails */
51: pc->failedreason = (PCFailedReason)err;
52: return(0);
53: }
55: ((PC_Factor*)dir)->fact = pc->pmat;
56: } else {
57: MatInfo info;
59: if (!pc->setupcalled) {
60: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
61: /* check if dir->row == dir->col */
62: ISEqual(dir->row,dir->col,&flg);
63: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"row and column permutations must equal");
64: ISDestroy(&dir->col); /* only pass one ordering into CholeskyFactor */
66: flg = PETSC_FALSE;
67: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
68: if (flg) {
69: PetscReal tol = 1.e-10;
70: PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
71: MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
72: }
73: if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
74: if (!((PC_Factor*)dir)->fact) {
75: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
76: }
77: MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
78: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
79: dir->hdr.actualfill = info.fill_ratio_needed;
80: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
81: } else if (pc->flag != SAME_NONZERO_PATTERN) {
82: if (!dir->hdr.reuseordering) {
83: ISDestroy(&dir->row);
84: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
85: ISDestroy(&dir->col); /* only use dir->row ordering in CholeskyFactor */
87: flg = PETSC_FALSE;
88: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
89: if (flg) {
90: PetscReal tol = 1.e-10;
91: PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
92: MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
93: }
94: if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
95: }
96: MatDestroy(&((PC_Factor*)dir)->fact);
97: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
98: MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
99: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
100: dir->hdr.actualfill = info.fill_ratio_needed;
101: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
102: } else {
103: MatFactorGetError(((PC_Factor*)dir)->fact,&err);
104: if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) {
105: MatFactorClearError(((PC_Factor*)dir)->fact);
106: pc->failedreason = PC_NOERROR;
107: }
108: }
109: MatFactorGetError(((PC_Factor*)dir)->fact,&err);
110: if (err) { /* FactorSymbolic() fails */
111: pc->failedreason = (PCFailedReason)err;
112: return(0);
113: }
115: MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
116: MatFactorGetError(((PC_Factor*)dir)->fact,&err);
117: if (err) { /* FactorNumeric() fails */
118: pc->failedreason = (PCFailedReason)err;
119: }
120: }
122: PCFactorGetMatSolverType(pc,&stype);
123: if (!stype) {
124: MatSolverType solverpackage;
125: MatFactorGetSolverType(((PC_Factor*)dir)->fact,&solverpackage);
126: PCFactorSetMatSolverType(pc,solverpackage);
127: }
128: return(0);
129: }
131: static PetscErrorCode PCReset_Cholesky(PC pc)
132: {
133: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
137: if (!dir->hdr.inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
138: ISDestroy(&dir->row);
139: ISDestroy(&dir->col);
140: return(0);
141: }
143: static PetscErrorCode PCDestroy_Cholesky(PC pc)
144: {
145: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
149: PCReset_Cholesky(pc);
150: PetscFree(((PC_Factor*)dir)->ordering);
151: PetscFree(((PC_Factor*)dir)->solvertype);
152: PetscFree(pc->data);
153: return(0);
154: }
156: static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y)
157: {
158: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
162: if (dir->hdr.inplace) {
163: MatSolve(pc->pmat,x,y);
164: } else {
165: MatSolve(((PC_Factor*)dir)->fact,x,y);
166: }
167: return(0);
168: }
170: static PetscErrorCode PCApplySymmetricLeft_Cholesky(PC pc,Vec x,Vec y)
171: {
172: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
176: if (dir->hdr.inplace) {
177: MatForwardSolve(pc->pmat,x,y);
178: } else {
179: MatForwardSolve(((PC_Factor*)dir)->fact,x,y);
180: }
181: return(0);
182: }
184: static PetscErrorCode PCApplySymmetricRight_Cholesky(PC pc,Vec x,Vec y)
185: {
186: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
190: if (dir->hdr.inplace) {
191: MatBackwardSolve(pc->pmat,x,y);
192: } else {
193: MatBackwardSolve(((PC_Factor*)dir)->fact,x,y);
194: }
195: return(0);
196: }
198: static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y)
199: {
200: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
204: if (dir->hdr.inplace) {
205: MatSolveTranspose(pc->pmat,x,y);
206: } else {
207: MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);
208: }
209: return(0);
210: }
212: /* -----------------------------------------------------------------------------------*/
214: /* -----------------------------------------------------------------------------------*/
216: /*@
217: PCFactorSetReuseOrdering - When similar matrices are factored, this
218: causes the ordering computed in the first factor to be used for all
219: following factors.
221: Logically Collective on PC
223: Input Parameters:
224: + pc - the preconditioner context
225: - flag - PETSC_TRUE to reuse else PETSC_FALSE
227: Options Database Key:
228: . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
230: Level: intermediate
232: .seealso: PCFactorSetReuseFill()
233: @*/
234: PetscErrorCode PCFactorSetReuseOrdering(PC pc,PetscBool flag)
235: {
241: PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));
242: return(0);
243: }
245: /*MC
246: PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner
248: Options Database Keys:
249: + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
250: . -pc_factor_mat_solver_type - Actives PCFactorSetMatSolverType() to choose the direct solver, like superlu
251: . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
252: . -pc_factor_fill <fill> - Sets fill amount
253: . -pc_factor_in_place - Activates in-place factorization
254: - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
256: Notes:
257: Not all options work for all matrix formats
259: Level: beginner
261: Notes:
262: Usually this will compute an "exact" solution in one iteration and does
263: not need a Krylov method (i.e. you can use -ksp_type preonly, or
264: KSPSetType(ksp,KSPPREONLY) for the Krylov method
266: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
267: PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
268: PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount()
269: PCFactorSetUseInPlace(), PCFactorGetUseInPlace(), PCFactorSetMatOrderingType()
271: M*/
273: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
274: {
276: PC_Cholesky *dir;
279: PetscNewLog(pc,&dir);
280: pc->data = (void*)dir;
281: PCFactorInitialize(pc);
283: ((PC_Factor*)dir)->factortype = MAT_FACTOR_CHOLESKY;
284: ((PC_Factor*)dir)->info.fill = 5.0;
286: dir->col = 0;
287: dir->row = 0;
289: /* MATORDERINGNATURAL_OR_ND allows selecting type based on matrix type sbaij or aij */
290: PetscStrallocpy(MATORDERINGNATURAL_OR_ND,(char**)&((PC_Factor*)dir)->ordering);
292: pc->ops->destroy = PCDestroy_Cholesky;
293: pc->ops->reset = PCReset_Cholesky;
294: pc->ops->apply = PCApply_Cholesky;
295: pc->ops->applysymmetricleft = PCApplySymmetricLeft_Cholesky;
296: pc->ops->applysymmetricright = PCApplySymmetricRight_Cholesky;
297: pc->ops->applytranspose = PCApplyTranspose_Cholesky;
298: pc->ops->setup = PCSetUp_Cholesky;
299: pc->ops->setfromoptions = PCSetFromOptions_Cholesky;
300: pc->ops->view = PCView_Factor;
301: pc->ops->applyrichardson = 0;
302: return(0);
303: }