Actual source code: cholesky.c
petsc-3.7.7 2017-09-25
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> /*I "petscpc.h" I*/
9: typedef struct {
10: PC_Factor hdr;
11: PetscReal actualfill; /* actual fill in factor */
12: PetscBool inplace; /* flag indicating in-place factorization */
13: IS row,col; /* index sets used for reordering */
14: PetscBool reuseordering; /* reuses previous reordering computed */
15: PetscBool reusefill; /* reuse fill from previous Cholesky */
16: } PC_Cholesky;
20: static PetscErrorCode PCFactorSetReuseOrdering_Cholesky(PC pc,PetscBool flag)
21: {
22: PC_Cholesky *lu = (PC_Cholesky*)pc->data;
25: lu->reuseordering = flag;
26: return(0);
27: }
31: static PetscErrorCode PCFactorSetReuseFill_Cholesky(PC pc,PetscBool flag)
32: {
33: PC_Cholesky *lu = (PC_Cholesky*)pc->data;
36: lu->reusefill = flag;
37: return(0);
38: }
42: static PetscErrorCode PCSetFromOptions_Cholesky(PetscOptionItems *PetscOptionsObject,PC pc)
43: {
47: PetscOptionsHead(PetscOptionsObject,"Cholesky options");
48: PCSetFromOptions_Factor(PetscOptionsObject,pc);
49: PetscOptionsTail();
50: return(0);
51: }
55: static PetscErrorCode PCView_Cholesky(PC pc,PetscViewer viewer)
56: {
57: PC_Cholesky *chol = (PC_Cholesky*)pc->data;
59: PetscBool iascii;
62: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
63: if (iascii) {
64: if (chol->inplace) {
65: PetscViewerASCIIPrintf(viewer," Cholesky: in-place factorization\n");
66: } else {
67: PetscViewerASCIIPrintf(viewer," Cholesky: out-of-place factorization\n");
68: }
70: if (chol->reusefill) {PetscViewerASCIIPrintf(viewer," Reusing fill from past factorization\n");}
71: if (chol->reuseordering) {PetscViewerASCIIPrintf(viewer," Reusing reordering from past factorization\n");}
72: }
73: PCView_Factor(pc,viewer);
74: return(0);
75: }
79: static PetscErrorCode PCSetUp_Cholesky(PC pc)
80: {
81: PetscErrorCode ierr;
82: PetscBool flg;
83: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
84: const MatSolverPackage stype;
87: pc->failedreason = PC_NOERROR;
88: if (dir->reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->actualfill;
90: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
91: if (dir->inplace) {
92: if (dir->row && dir->col && (dir->row != dir->col)) {
93: ISDestroy(&dir->row);
94: }
95: ISDestroy(&dir->col);
96: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
97: if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */
98: ISDestroy(&dir->col);
99: }
100: if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
101: MatCholeskyFactor(pc->pmat,dir->row,&((PC_Factor*)dir)->info);
102: if (pc->pmat->errortype) { /* Factor() fails */
103: pc->failedreason = (PCFailedReason)pc->pmat->errortype;
104: return(0);
105: }
107: ((PC_Factor*)dir)->fact = pc->pmat;
108: } else {
109: MatInfo info;
110: Mat F;
111: if (!pc->setupcalled) {
112: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
113: /* check if dir->row == dir->col */
114: ISEqual(dir->row,dir->col,&flg);
115: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"row and column permutations must equal");
116: ISDestroy(&dir->col); /* only pass one ordering into CholeskyFactor */
118: flg = PETSC_FALSE;
119: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
120: if (flg) {
121: PetscReal tol = 1.e-10;
122: PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
123: MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
124: }
125: if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
126: if (!((PC_Factor*)dir)->fact) {
127: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
128: }
129: MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
130: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
131: dir->actualfill = info.fill_ratio_needed;
132: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
133: } else if (pc->flag != SAME_NONZERO_PATTERN) {
134: if (!dir->reuseordering) {
135: ISDestroy(&dir->row);
136: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
137: ISDestroy(&dir->col); /* only use dir->row ordering in CholeskyFactor */
139: flg = PETSC_FALSE;
140: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
141: if (flg) {
142: PetscReal tol = 1.e-10;
143: PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
144: MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
145: }
146: if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
147: }
148: MatDestroy(&((PC_Factor*)dir)->fact);
149: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
150: MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
151: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
152: dir->actualfill = info.fill_ratio_needed;
153: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
154: } else {
155: F = ((PC_Factor*)dir)->fact;
156: if ((PCFailedReason)F->errortype == PC_FACTOR_NUMERIC_ZEROPIVOT) {
157: F->errortype = MAT_FACTOR_NOERROR;
158: pc->failedreason = (PCFailedReason)F->errortype;
159: }
160: }
161: F = ((PC_Factor*)dir)->fact;
162: if (F->errortype) { /* FactorSymbolic() fails */
163: pc->failedreason = (PCFailedReason)F->errortype;
164: return(0);
165: }
167: MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
168: if (F->errortype) { /* FactorNumeric() fails */
169: pc->failedreason = (PCFailedReason)F->errortype;
170: }
171: }
173: PCFactorGetMatSolverPackage(pc,&stype);
174: if (!stype) {
175: PCFactorSetMatSolverPackage(pc,((PC_Factor*)dir)->fact->solvertype);
176: }
177: return(0);
178: }
182: static PetscErrorCode PCReset_Cholesky(PC pc)
183: {
184: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
188: if (!dir->inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
189: ISDestroy(&dir->row);
190: ISDestroy(&dir->col);
191: return(0);
192: }
196: static PetscErrorCode PCDestroy_Cholesky(PC pc)
197: {
198: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
202: PCReset_Cholesky(pc);
203: PetscFree(((PC_Factor*)dir)->ordering);
204: PetscFree(((PC_Factor*)dir)->solvertype);
205: PetscFree(pc->data);
206: return(0);
207: }
211: static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y)
212: {
213: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
217: if (dir->inplace) {
218: MatSolve(pc->pmat,x,y);
219: } else {
220: MatSolve(((PC_Factor*)dir)->fact,x,y);
221: }
222: return(0);
223: }
227: static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y)
228: {
229: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
233: if (dir->inplace) {
234: MatSolveTranspose(pc->pmat,x,y);
235: } else {
236: MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);
237: }
238: return(0);
239: }
241: /* -----------------------------------------------------------------------------------*/
245: static PetscErrorCode PCFactorSetUseInPlace_Cholesky(PC pc,PetscBool flg)
246: {
247: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
250: dir->inplace = flg;
251: return(0);
252: }
256: static PetscErrorCode PCFactorGetUseInPlace_Cholesky(PC pc,PetscBool *flg)
257: {
258: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
261: *flg = dir->inplace;
262: return(0);
263: }
265: /* -----------------------------------------------------------------------------------*/
269: /*@
270: PCFactorSetReuseOrdering - When similar matrices are factored, this
271: causes the ordering computed in the first factor to be used for all
272: following factors.
274: Logically Collective on PC
276: Input Parameters:
277: + pc - the preconditioner context
278: - flag - PETSC_TRUE to reuse else PETSC_FALSE
280: Options Database Key:
281: . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
283: Level: intermediate
285: .keywords: PC, levels, reordering, factorization, incomplete, LU
287: .seealso: PCFactorSetReuseFill()
288: @*/
289: PetscErrorCode PCFactorSetReuseOrdering(PC pc,PetscBool flag)
290: {
296: PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));
297: return(0);
298: }
300: /*MC
301: PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner
303: Options Database Keys:
304: + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
305: . -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like superlu
306: . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
307: . -pc_factor_fill <fill> - Sets fill amount
308: . -pc_factor_in_place - Activates in-place factorization
309: - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
311: Notes: Not all options work for all matrix formats
313: Level: beginner
315: Concepts: Cholesky factorization, direct solver
317: Notes: Usually this will compute an "exact" solution in one iteration and does
318: not need a Krylov method (i.e. you can use -ksp_type preonly, or
319: KSPSetType(ksp,KSPPREONLY) for the Krylov method
321: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
322: PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
323: PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount()
324: PCFactorSetUseInPlace(), PCFactorGetUseInPlace(), PCFactorSetMatOrderingType()
326: M*/
330: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
331: {
333: PC_Cholesky *dir;
336: PetscNewLog(pc,&dir);
338: ((PC_Factor*)dir)->fact = 0;
339: dir->inplace = PETSC_FALSE;
341: MatFactorInfoInitialize(&((PC_Factor*)dir)->info);
343: ((PC_Factor*)dir)->factortype = MAT_FACTOR_CHOLESKY;
344: ((PC_Factor*)dir)->info.fill = 5.0;
345: ((PC_Factor*)dir)->info.shifttype = (PetscReal) MAT_SHIFT_NONE;
346: ((PC_Factor*)dir)->info.shiftamount = 0.0;
347: ((PC_Factor*)dir)->info.zeropivot = 100.0*PETSC_MACHINE_EPSILON;
348: ((PC_Factor*)dir)->info.pivotinblocks = 1.0;
350: dir->col = 0;
351: dir->row = 0;
353: PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)dir)->ordering);
354: dir->reusefill = PETSC_FALSE;
355: dir->reuseordering = PETSC_FALSE;
356: pc->data = (void*)dir;
358: pc->ops->destroy = PCDestroy_Cholesky;
359: pc->ops->reset = PCReset_Cholesky;
360: pc->ops->apply = PCApply_Cholesky;
361: pc->ops->applytranspose = PCApplyTranspose_Cholesky;
362: pc->ops->setup = PCSetUp_Cholesky;
363: pc->ops->setfromoptions = PCSetFromOptions_Cholesky;
364: pc->ops->view = PCView_Cholesky;
365: pc->ops->applyrichardson = 0;
366: pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor;
368: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetZeroPivot_C",PCFactorSetZeroPivot_Factor);
369: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetZeroPivot_C",PCFactorGetZeroPivot_Factor);
370: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Factor);
371: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetShiftType_C",PCFactorGetShiftType_Factor);
372: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftAmount_C",PCFactorSetShiftAmount_Factor);
373: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetShiftAmount_C",PCFactorGetShiftAmount_Factor);
374: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C",PCFactorSetUpMatSolverPackage_Factor);
375: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C",PCFactorSetUpMatSolverPackage_Factor);
376: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatSolverPackage_C",PCFactorSetMatSolverPackage_Factor);
377: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetMatSolverPackage_C",PCFactorGetMatSolverPackage_Factor);
378: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetFill_C",PCFactorSetFill_Factor);
379: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUseInPlace_C",PCFactorSetUseInPlace_Cholesky);
380: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetUseInPlace_C",PCFactorGetUseInPlace_Cholesky);
381: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatOrderingType_C",PCFactorSetMatOrderingType_Factor);
382: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C",PCFactorSetReuseOrdering_Cholesky);
383: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseFill_C",PCFactorSetReuseFill_Cholesky);
384: return(0);
385: }