Actual source code: cholesky.c
petsc-3.6.1 2015-08-06
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(PetscOptions *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: }
80: static PetscErrorCode PCSetUp_Cholesky(PC pc)
81: {
83: PetscBool flg;
84: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
87: if (dir->reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->actualfill;
89: if (dir->inplace) {
90: if (dir->row && dir->col && (dir->row != dir->col)) {
91: ISDestroy(&dir->row);
92: }
93: ISDestroy(&dir->col);
94: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
95: if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */
96: ISDestroy(&dir->col);
97: }
98: if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
99: MatCholeskyFactor(pc->pmat,dir->row,&((PC_Factor*)dir)->info);
101: ((PC_Factor*)dir)->fact = pc->pmat;
102: } else {
103: MatInfo info;
104: if (!pc->setupcalled) {
105: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
106: /* check if dir->row == dir->col */
107: ISEqual(dir->row,dir->col,&flg);
108: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"row and column permutations must equal");
109: ISDestroy(&dir->col); /* only pass one ordering into CholeskyFactor */
111: flg = PETSC_FALSE;
112: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
113: if (flg) {
114: PetscReal tol = 1.e-10;
115: PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
116: MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
117: }
118: if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
119: if (!((PC_Factor*)dir)->fact) {
120: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
121: }
122: MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
123: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
124: dir->actualfill = info.fill_ratio_needed;
125: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
126: } else if (pc->flag != SAME_NONZERO_PATTERN) {
127: if (!dir->reuseordering) {
128: ISDestroy(&dir->row);
129: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
130: ISDestroy(&dir->col); /* only use dir->row ordering in CholeskyFactor */
132: flg = PETSC_FALSE;
133: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
134: if (flg) {
135: PetscReal tol = 1.e-10;
136: PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
137: MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
138: }
139: if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
140: }
141: MatDestroy(&((PC_Factor*)dir)->fact);
142: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
143: MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
144: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
145: dir->actualfill = info.fill_ratio_needed;
146: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
147: }
148: MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
149: }
150: return(0);
151: }
155: static PetscErrorCode PCReset_Cholesky(PC pc)
156: {
157: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
161: if (!dir->inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
162: ISDestroy(&dir->row);
163: ISDestroy(&dir->col);
164: return(0);
165: }
169: static PetscErrorCode PCDestroy_Cholesky(PC pc)
170: {
171: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
175: PCReset_Cholesky(pc);
176: PetscFree(((PC_Factor*)dir)->ordering);
177: PetscFree(((PC_Factor*)dir)->solvertype);
178: PetscFree(pc->data);
179: return(0);
180: }
184: static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y)
185: {
186: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
190: if (dir->inplace) {
191: MatSolve(pc->pmat,x,y);
192: } else {
193: MatSolve(((PC_Factor*)dir)->fact,x,y);
194: }
195: return(0);
196: }
200: static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y)
201: {
202: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
206: if (dir->inplace) {
207: MatSolveTranspose(pc->pmat,x,y);
208: } else {
209: MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);
210: }
211: return(0);
212: }
214: /* -----------------------------------------------------------------------------------*/
218: static PetscErrorCode PCFactorSetUseInPlace_Cholesky(PC pc,PetscBool flg)
219: {
220: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
223: dir->inplace = flg;
224: return(0);
225: }
229: static PetscErrorCode PCFactorGetUseInPlace_Cholesky(PC pc,PetscBool *flg)
230: {
231: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
234: *flg = dir->inplace;
235: return(0);
236: }
238: /* -----------------------------------------------------------------------------------*/
242: /*@
243: PCFactorSetReuseOrdering - When similar matrices are factored, this
244: causes the ordering computed in the first factor to be used for all
245: following factors.
247: Logically Collective on PC
249: Input Parameters:
250: + pc - the preconditioner context
251: - flag - PETSC_TRUE to reuse else PETSC_FALSE
253: Options Database Key:
254: . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
256: Level: intermediate
258: .keywords: PC, levels, reordering, factorization, incomplete, LU
260: .seealso: PCFactorSetReuseFill()
261: @*/
262: PetscErrorCode PCFactorSetReuseOrdering(PC pc,PetscBool flag)
263: {
269: PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));
270: return(0);
271: }
273: /*MC
274: PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner
276: Options Database Keys:
277: + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
278: . -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like superlu
279: . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
280: . -pc_factor_fill <fill> - Sets fill amount
281: . -pc_factor_in_place - Activates in-place factorization
282: - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
284: Notes: Not all options work for all matrix formats
286: Level: beginner
288: Concepts: Cholesky factorization, direct solver
290: Notes: Usually this will compute an "exact" solution in one iteration and does
291: not need a Krylov method (i.e. you can use -ksp_type preonly, or
292: KSPSetType(ksp,KSPPREONLY) for the Krylov method
294: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
295: PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
296: PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount()
297: PCFactorSetUseInPlace(), PCFactorGetUseInPlace(), PCFactorSetMatOrderingType()
299: M*/
303: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
304: {
306: PC_Cholesky *dir;
309: PetscNewLog(pc,&dir);
311: ((PC_Factor*)dir)->fact = 0;
312: dir->inplace = PETSC_FALSE;
314: MatFactorInfoInitialize(&((PC_Factor*)dir)->info);
316: ((PC_Factor*)dir)->factortype = MAT_FACTOR_CHOLESKY;
317: ((PC_Factor*)dir)->info.fill = 5.0;
318: ((PC_Factor*)dir)->info.shifttype = (PetscReal) MAT_SHIFT_NONE;
319: ((PC_Factor*)dir)->info.shiftamount = 0.0;
320: ((PC_Factor*)dir)->info.zeropivot = 100.0*PETSC_MACHINE_EPSILON;
321: ((PC_Factor*)dir)->info.pivotinblocks = 1.0;
323: dir->col = 0;
324: dir->row = 0;
326: PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)dir)->ordering);
327: PetscStrallocpy(MATSOLVERPETSC,&((PC_Factor*)dir)->solvertype);
329: dir->reusefill = PETSC_FALSE;
330: dir->reuseordering = PETSC_FALSE;
331: pc->data = (void*)dir;
333: pc->ops->destroy = PCDestroy_Cholesky;
334: pc->ops->reset = PCReset_Cholesky;
335: pc->ops->apply = PCApply_Cholesky;
336: pc->ops->applytranspose = PCApplyTranspose_Cholesky;
337: pc->ops->setup = PCSetUp_Cholesky;
338: pc->ops->setfromoptions = PCSetFromOptions_Cholesky;
339: pc->ops->view = PCView_Cholesky;
340: pc->ops->applyrichardson = 0;
341: pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor;
343: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C",PCFactorSetUpMatSolverPackage_Factor);
344: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatSolverPackage_C",PCFactorSetMatSolverPackage_Factor);
345: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetMatSolverPackage_C",PCFactorGetMatSolverPackage_Factor);
346: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetZeroPivot_C",PCFactorSetZeroPivot_Factor);
347: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Factor);
348: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftAmount_C",PCFactorSetShiftAmount_Factor);
349: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetFill_C",PCFactorSetFill_Factor);
350: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUseInPlace_C",PCFactorSetUseInPlace_Cholesky);
351: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetUseInPlace_C",PCFactorGetUseInPlace_Cholesky);
352: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatOrderingType_C",PCFactorSetMatOrderingType_Factor);
353: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C",PCFactorSetReuseOrdering_Cholesky);
354: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseFill_C",PCFactorSetReuseFill_Cholesky);
355: return(0);
356: }