Actual source code: cholesky.c
petsc-3.4.5 2014-06-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> /*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;
25: lu = (PC_Cholesky*)pc->data;
26: lu->reuseordering = flag;
27: return(0);
28: }
32: static PetscErrorCode PCFactorSetReuseFill_Cholesky(PC pc,PetscBool flag)
33: {
34: PC_Cholesky *lu;
37: lu = (PC_Cholesky*)pc->data;
38: lu->reusefill = flag;
39: return(0);
40: }
44: static PetscErrorCode PCSetFromOptions_Cholesky(PC pc)
45: {
49: PetscOptionsHead("Cholesky options");
50: PCSetFromOptions_Factor(pc);
51: PetscOptionsTail();
52: return(0);
53: }
57: static PetscErrorCode PCView_Cholesky(PC pc,PetscViewer viewer)
58: {
59: PC_Cholesky *chol = (PC_Cholesky*)pc->data;
61: PetscBool iascii;
64: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
65: if (iascii) {
66: if (chol->inplace) {
67: PetscViewerASCIIPrintf(viewer," Cholesky: in-place factorization\n");
68: } else {
69: PetscViewerASCIIPrintf(viewer," Cholesky: out-of-place factorization\n");
70: }
72: if (chol->reusefill) {PetscViewerASCIIPrintf(viewer," Reusing fill from past factorization\n");}
73: if (chol->reuseordering) {PetscViewerASCIIPrintf(viewer," Reusing reordering from past factorization\n");}
74: }
75: PCView_Factor(pc,viewer);
76: return(0);
77: }
82: static PetscErrorCode PCSetUp_Cholesky(PC pc)
83: {
85: PetscBool flg;
86: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
89: if (dir->reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->actualfill;
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(pc,dir->row);}
101: MatCholeskyFactor(pc->pmat,dir->row,&((PC_Factor*)dir)->info);
103: ((PC_Factor*)dir)->fact = pc->pmat;
104: } else {
105: MatInfo info;
106: if (!pc->setupcalled) {
107: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
108: /* check if dir->row == dir->col */
109: ISEqual(dir->row,dir->col,&flg);
110: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"row and column permutations must equal");
111: ISDestroy(&dir->col); /* only pass one ordering into CholeskyFactor */
113: flg = PETSC_FALSE;
114: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
115: if (flg) {
116: PetscReal tol = 1.e-10;
117: PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
118: MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
119: }
120: if (dir->row) {PetscLogObjectParent(pc,dir->row);}
121: if (!((PC_Factor*)dir)->fact) {
122: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
123: }
124: MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
125: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
126: dir->actualfill = info.fill_ratio_needed;
127: PetscLogObjectParent(pc,((PC_Factor*)dir)->fact);
128: } else if (pc->flag != SAME_NONZERO_PATTERN) {
129: if (!dir->reuseordering) {
130: ISDestroy(&dir->row);
131: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
132: ISDestroy(&dir->col); /* only use dir->row ordering in CholeskyFactor */
134: flg = PETSC_FALSE;
135: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
136: if (flg) {
137: PetscReal tol = 1.e-10;
138: PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
139: MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
140: }
141: if (dir->row) {PetscLogObjectParent(pc,dir->row);}
142: }
143: MatDestroy(&((PC_Factor*)dir)->fact);
144: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
145: MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
146: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
147: dir->actualfill = info.fill_ratio_needed;
148: PetscLogObjectParent(pc,((PC_Factor*)dir)->fact);
149: }
150: MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
151: }
152: return(0);
153: }
157: static PetscErrorCode PCReset_Cholesky(PC pc)
158: {
159: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
163: if (!dir->inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
164: ISDestroy(&dir->row);
165: ISDestroy(&dir->col);
166: return(0);
167: }
171: static PetscErrorCode PCDestroy_Cholesky(PC pc)
172: {
173: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
177: PCReset_Cholesky(pc);
178: PetscFree(((PC_Factor*)dir)->ordering);
179: PetscFree(((PC_Factor*)dir)->solvertype);
180: PetscFree(pc->data);
181: return(0);
182: }
186: static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y)
187: {
188: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
192: if (dir->inplace) {
193: MatSolve(pc->pmat,x,y);
194: } else {
195: MatSolve(((PC_Factor*)dir)->fact,x,y);
196: }
197: return(0);
198: }
202: static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y)
203: {
204: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
208: if (dir->inplace) {
209: MatSolveTranspose(pc->pmat,x,y);
210: } else {
211: MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);
212: }
213: return(0);
214: }
216: /* -----------------------------------------------------------------------------------*/
220: static PetscErrorCode PCFactorSetUseInPlace_Cholesky(PC pc)
221: {
222: PC_Cholesky *dir;
225: dir = (PC_Cholesky*)pc->data;
226: dir->inplace = PETSC_TRUE;
227: return(0);
228: }
230: /* -----------------------------------------------------------------------------------*/
234: /*@
235: PCFactorSetReuseOrdering - When similar matrices are factored, this
236: causes the ordering computed in the first factor to be used for all
237: following factors.
239: Logically Collective on PC
241: Input Parameters:
242: + pc - the preconditioner context
243: - flag - PETSC_TRUE to reuse else PETSC_FALSE
245: Options Database Key:
246: . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
248: Level: intermediate
250: .keywords: PC, levels, reordering, factorization, incomplete, LU
252: .seealso: PCFactorSetReuseFill()
253: @*/
254: PetscErrorCode PCFactorSetReuseOrdering(PC pc,PetscBool flag)
255: {
261: PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));
262: return(0);
263: }
266: /*MC
267: PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner
269: Options Database Keys:
270: + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
271: . -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like superlu
272: . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
273: . -pc_factor_fill <fill> - Sets fill amount
274: . -pc_factor_in_place - Activates in-place factorization
275: - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
277: Notes: Not all options work for all matrix formats
279: Level: beginner
281: Concepts: Cholesky factorization, direct solver
283: Notes: Usually this will compute an "exact" solution in one iteration and does
284: not need a Krylov method (i.e. you can use -ksp_type preonly, or
285: KSPSetType(ksp,KSPPREONLY) for the Krylov method
287: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
288: PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
289: PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount()
290: PCFactorSetUseInPlace(), PCFactorSetMatOrderingType()
292: M*/
296: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
297: {
299: PC_Cholesky *dir;
302: PetscNewLog(pc,PC_Cholesky,&dir);
304: ((PC_Factor*)dir)->fact = 0;
305: dir->inplace = PETSC_FALSE;
307: MatFactorInfoInitialize(&((PC_Factor*)dir)->info);
309: ((PC_Factor*)dir)->factortype = MAT_FACTOR_CHOLESKY;
310: ((PC_Factor*)dir)->info.fill = 5.0;
311: ((PC_Factor*)dir)->info.shifttype = (PetscReal) MAT_SHIFT_NONE;
312: ((PC_Factor*)dir)->info.shiftamount = 0.0;
313: ((PC_Factor*)dir)->info.zeropivot = 100.0*PETSC_MACHINE_EPSILON;
314: ((PC_Factor*)dir)->info.pivotinblocks = 1.0;
316: dir->col = 0;
317: dir->row = 0;
319: PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)dir)->ordering);
320: PetscStrallocpy(MATSOLVERPETSC,&((PC_Factor*)dir)->solvertype);
322: dir->reusefill = PETSC_FALSE;
323: dir->reuseordering = PETSC_FALSE;
324: pc->data = (void*)dir;
326: pc->ops->destroy = PCDestroy_Cholesky;
327: pc->ops->reset = PCReset_Cholesky;
328: pc->ops->apply = PCApply_Cholesky;
329: pc->ops->applytranspose = PCApplyTranspose_Cholesky;
330: pc->ops->setup = PCSetUp_Cholesky;
331: pc->ops->setfromoptions = PCSetFromOptions_Cholesky;
332: pc->ops->view = PCView_Cholesky;
333: pc->ops->applyrichardson = 0;
334: pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor;
336: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C",PCFactorSetUpMatSolverPackage_Factor);
337: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatSolverPackage_C",PCFactorSetMatSolverPackage_Factor);
338: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetMatSolverPackage_C",PCFactorGetMatSolverPackage_Factor);
339: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetZeroPivot_C",PCFactorSetZeroPivot_Factor);
340: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Factor);
341: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftAmount_C",PCFactorSetShiftAmount_Factor);
342: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetFill_C",PCFactorSetFill_Factor);
343: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUseInPlace_C",PCFactorSetUseInPlace_Cholesky);
344: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatOrderingType_C",PCFactorSetMatOrderingType_Factor);
345: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C",PCFactorSetReuseOrdering_Cholesky);
346: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseFill_C",PCFactorSetReuseFill_Cholesky);
347: return(0);
348: }