Actual source code: cholesky.c
petsc-3.3-p7 2013-05-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> /*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;
18: EXTERN_C_BEGIN
21: PetscErrorCode PCFactorSetReuseOrdering_Cholesky(PC pc,PetscBool flag)
22: {
23: PC_Cholesky *lu;
24:
26: lu = (PC_Cholesky*)pc->data;
27: lu->reuseordering = flag;
28: return(0);
29: }
30: EXTERN_C_END
32: EXTERN_C_BEGIN
35: PetscErrorCode PCFactorSetReuseFill_Cholesky(PC pc,PetscBool flag)
36: {
37: PC_Cholesky *lu;
38:
40: lu = (PC_Cholesky*)pc->data;
41: lu->reusefill = flag;
42: return(0);
43: }
44: EXTERN_C_END
48: static PetscErrorCode PCSetFromOptions_Cholesky(PC pc)
49: {
51:
53: PetscOptionsHead("Cholesky options");
54: PCSetFromOptions_Factor(pc);
55: PetscOptionsTail();
56: return(0);
57: }
61: static PetscErrorCode PCView_Cholesky(PC pc,PetscViewer viewer)
62: {
63: PC_Cholesky *chol = (PC_Cholesky*)pc->data;
65: PetscBool iascii;
66:
68: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
69: if (iascii) {
70: if (chol->inplace) {
71: PetscViewerASCIIPrintf(viewer," Cholesky: in-place factorization\n");
72: } else {
73: PetscViewerASCIIPrintf(viewer," Cholesky: out-of-place factorization\n");
74: }
75:
76: if (chol->reusefill) {PetscViewerASCIIPrintf(viewer," Reusing fill from past factorization\n");}
77: if (chol->reuseordering) {PetscViewerASCIIPrintf(viewer," Reusing reordering from past factorization\n");}
78: }
79: PCView_Factor(pc,viewer);
80: return(0);
81: }
86: static PetscErrorCode PCSetUp_Cholesky(PC pc)
87: {
89: PetscBool flg;
90: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
93: if (dir->reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->actualfill;
94:
95: if (dir->inplace) {
96: if (dir->row && dir->col && (dir->row != dir->col)) {
97: ISDestroy(&dir->row);
98: }
99: ISDestroy(&dir->col);
100: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
101: if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */
102: ISDestroy(&dir->col);
103: }
104: if (dir->row) {PetscLogObjectParent(pc,dir->row);}
105: MatCholeskyFactor(pc->pmat,dir->row,&((PC_Factor*)dir)->info);
106: ((PC_Factor*)dir)->fact = pc->pmat;
107: } else {
108: MatInfo info;
109: if (!pc->setupcalled) {
110: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
111: /* check if dir->row == dir->col */
112: ISEqual(dir->row,dir->col,&flg);
113: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"row and column permutations must equal");
114: ISDestroy(&dir->col); /* only pass one ordering into CholeskyFactor */
116: flg = PETSC_FALSE;
117: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,PETSC_NULL);
118: if (flg) {
119: PetscReal tol = 1.e-10;
120: PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,PETSC_NULL);
121: MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
122: }
123: if (dir->row) {PetscLogObjectParent(pc,dir->row);}
124: if (!((PC_Factor*)dir)->fact){
125: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
126: }
127: MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
128: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
129: dir->actualfill = info.fill_ratio_needed;
130: PetscLogObjectParent(pc,((PC_Factor*)dir)->fact);
131: } else if (pc->flag != SAME_NONZERO_PATTERN) {
132: if (!dir->reuseordering) {
133: ISDestroy(&dir->row);
134: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
135: ISDestroy(&dir->col); /* only use dir->row ordering in CholeskyFactor */
137: flg = PETSC_FALSE;
138: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,PETSC_NULL);
139: if (flg) {
140: PetscReal tol = 1.e-10;
141: PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,PETSC_NULL);
142: MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
143: }
144: if (dir->row) {PetscLogObjectParent(pc,dir->row);}
145: }
146: MatDestroy(&((PC_Factor*)dir)->fact);
147: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
148: MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
149: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
150: dir->actualfill = info.fill_ratio_needed;
151: PetscLogObjectParent(pc,((PC_Factor*)dir)->fact);
152: }
153: MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
154: }
155: return(0);
156: }
160: static PetscErrorCode PCReset_Cholesky(PC pc)
161: {
162: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
166: if (!dir->inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
167: ISDestroy(&dir->row);
168: ISDestroy(&dir->col);
169: return(0);
170: }
174: static PetscErrorCode PCDestroy_Cholesky(PC pc)
175: {
176: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
180: PCReset_Cholesky(pc);
181: PetscFree(((PC_Factor*)dir)->ordering);
182: PetscFree(((PC_Factor*)dir)->solvertype);
183: PetscFree(pc->data);
184: return(0);
185: }
189: static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y)
190: {
191: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
193:
195: if (dir->inplace) {MatSolve(pc->pmat,x,y);}
196: else {MatSolve(((PC_Factor*)dir)->fact,x,y);}
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) {MatSolveTranspose(pc->pmat,x,y);}
209: else {MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);}
210: return(0);
211: }
213: /* -----------------------------------------------------------------------------------*/
215: EXTERN_C_BEGIN
218: PetscErrorCode PCFactorSetUseInPlace_Cholesky(PC pc)
219: {
220: PC_Cholesky *dir;
223: dir = (PC_Cholesky*)pc->data;
224: dir->inplace = PETSC_TRUE;
225: return(0);
226: }
227: EXTERN_C_END
229: /* -----------------------------------------------------------------------------------*/
233: /*@
234: PCFactorSetReuseOrdering - When similar matrices are factored, this
235: causes the ordering computed in the first factor to be used for all
236: following factors.
238: Logically Collective on PC
240: Input Parameters:
241: + pc - the preconditioner context
242: - flag - PETSC_TRUE to reuse else PETSC_FALSE
244: Options Database Key:
245: . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
247: Level: intermediate
249: .keywords: PC, levels, reordering, factorization, incomplete, LU
251: .seealso: PCFactorSetReuseFill()
252: @*/
253: PetscErrorCode PCFactorSetReuseOrdering(PC pc,PetscBool flag)
254: {
260: PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));
261: return(0);
262: }
265: /*MC
266: PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner
268: Options Database Keys:
269: + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
270: . -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like spooles
271: . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
272: . -pc_factor_fill <fill> - Sets fill amount
273: . -pc_factor_in_place - Activates in-place factorization
274: - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
276: Notes: Not all options work for all matrix formats
278: Level: beginner
280: Concepts: Cholesky factorization, direct solver
282: Notes: Usually this will compute an "exact" solution in one iteration and does
283: not need a Krylov method (i.e. you can use -ksp_type preonly, or
284: KSPSetType(ksp,KSPPREONLY) for the Krylov method
286: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
287: PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
288: PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount()
289: PCFactorSetUseInPlace(), PCFactorSetMatOrderingType()
291: M*/
293: EXTERN_C_BEGIN
296: 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;
306: MatFactorInfoInitialize(&((PC_Factor*)dir)->info);
307: ((PC_Factor*)dir)->factortype = MAT_FACTOR_CHOLESKY;
308: ((PC_Factor*)dir)->info.fill = 5.0;
309: ((PC_Factor*)dir)->info.shifttype = (PetscReal) MAT_SHIFT_NONE;
310: ((PC_Factor*)dir)->info.shiftamount = 0.0;
311: ((PC_Factor*)dir)->info.zeropivot = 100.0*PETSC_MACHINE_EPSILON;
312: ((PC_Factor*)dir)->info.pivotinblocks = 1.0;
313: dir->col = 0;
314: dir->row = 0;
315: PetscStrallocpy(MATORDERINGNATURAL,&((PC_Factor*)dir)->ordering);
316: PetscStrallocpy(MATSOLVERPETSC,&((PC_Factor*)dir)->solvertype);
317: dir->reusefill = PETSC_FALSE;
318: dir->reuseordering = PETSC_FALSE;
319: pc->data = (void*)dir;
321: pc->ops->destroy = PCDestroy_Cholesky;
322: pc->ops->reset = PCReset_Cholesky;
323: pc->ops->apply = PCApply_Cholesky;
324: pc->ops->applytranspose = PCApplyTranspose_Cholesky;
325: pc->ops->setup = PCSetUp_Cholesky;
326: pc->ops->setfromoptions = PCSetFromOptions_Cholesky;
327: pc->ops->view = PCView_Cholesky;
328: pc->ops->applyrichardson = 0;
329: pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor;
331: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C","PCFactorSetUpMatSolverPackage_Factor",
332: PCFactorSetUpMatSolverPackage_Factor);
333: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatSolverPackage_C","PCFactorSetMatSolverPackage_Factor",
334: PCFactorSetMatSolverPackage_Factor);
335: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorGetMatSolverPackage_C","PCFactorGetMatSolverPackage_Factor",
336: PCFactorGetMatSolverPackage_Factor);
337: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Factor",
338: PCFactorSetZeroPivot_Factor);
339: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftType_C","PCFactorSetShiftType_Factor",
340: PCFactorSetShiftType_Factor);
341: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftAmount_C","PCFactorSetShiftAmount_Factor",
342: PCFactorSetShiftAmount_Factor);
343: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Factor",
344: PCFactorSetFill_Factor);
345: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUseInPlace_C","PCFactorSetUseInPlace_Cholesky",
346: PCFactorSetUseInPlace_Cholesky);
347: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Factor",
348: PCFactorSetMatOrderingType_Factor);
349: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseOrdering_C","PCFactorSetReuseOrdering_Cholesky",
350: PCFactorSetReuseOrdering_Cholesky);
351: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseFill_C","PCFactorSetReuseFill_Cholesky",
352: PCFactorSetReuseFill_Cholesky);
353: return(0);
354: }
355: EXTERN_C_END