Actual source code: cholesky.c
petsc-3.7.3 2016-08-01
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: {
82: PetscBool flg;
83: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
84: const MatSolverPackage stype;
87: if (dir->reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->actualfill;
89: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
90: if (dir->inplace) {
91: if (dir->row && dir->col && (dir->row != dir->col)) {
92: ISDestroy(&dir->row);
93: }
94: ISDestroy(&dir->col);
95: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
96: if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */
97: ISDestroy(&dir->col);
98: }
99: if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
100: MatCholeskyFactor(pc->pmat,dir->row,&((PC_Factor*)dir)->info);
101: if (pc->pmat->errortype) { /* Factor() fails */
102: pc->failedreason = (PCFailedReason)pc->pmat->errortype;
103: return(0);
104: }
106: ((PC_Factor*)dir)->fact = pc->pmat;
107: } else {
108: MatInfo info;
109: Mat F;
110: if (!pc->setupcalled) {
111: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
112: /* check if dir->row == dir->col */
113: ISEqual(dir->row,dir->col,&flg);
114: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"row and column permutations must equal");
115: ISDestroy(&dir->col); /* only pass one ordering into CholeskyFactor */
117: flg = PETSC_FALSE;
118: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
119: if (flg) {
120: PetscReal tol = 1.e-10;
121: PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
122: MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
123: }
124: if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
125: if (!((PC_Factor*)dir)->fact) {
126: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
127: }
128: MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
129: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
130: dir->actualfill = info.fill_ratio_needed;
131: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
132: } else if (pc->flag != SAME_NONZERO_PATTERN) {
133: if (!dir->reuseordering) {
134: ISDestroy(&dir->row);
135: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
136: ISDestroy(&dir->col); /* only use dir->row ordering in CholeskyFactor */
138: flg = PETSC_FALSE;
139: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
140: if (flg) {
141: PetscReal tol = 1.e-10;
142: PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
143: MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
144: }
145: if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
146: }
147: MatDestroy(&((PC_Factor*)dir)->fact);
148: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
149: MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
150: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
151: dir->actualfill = info.fill_ratio_needed;
152: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
153: }
154: F = ((PC_Factor*)dir)->fact;
155: if (F->errortype) { /* FactorSymbolic() fails */
156: pc->failedreason = (PCFailedReason)F->errortype;
157: return(0);
158: }
160: MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
161: if (F->errortype) { /* FactorNumeric() fails */
162: pc->failedreason = (PCFailedReason)F->errortype;
163: }
164: }
166: PCFactorGetMatSolverPackage(pc,&stype);
167: if (!stype) {
168: PCFactorSetMatSolverPackage(pc,((PC_Factor*)dir)->fact->solvertype);
169: }
170: return(0);
171: }
175: static PetscErrorCode PCReset_Cholesky(PC pc)
176: {
177: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
181: if (!dir->inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
182: ISDestroy(&dir->row);
183: ISDestroy(&dir->col);
184: return(0);
185: }
189: static PetscErrorCode PCDestroy_Cholesky(PC pc)
190: {
191: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
195: PCReset_Cholesky(pc);
196: PetscFree(((PC_Factor*)dir)->ordering);
197: PetscFree(((PC_Factor*)dir)->solvertype);
198: PetscFree(pc->data);
199: return(0);
200: }
204: static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y)
205: {
206: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
210: if (dir->inplace) {
211: MatSolve(pc->pmat,x,y);
212: } else {
213: MatSolve(((PC_Factor*)dir)->fact,x,y);
214: }
215: return(0);
216: }
220: static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y)
221: {
222: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
226: if (dir->inplace) {
227: MatSolveTranspose(pc->pmat,x,y);
228: } else {
229: MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);
230: }
231: return(0);
232: }
234: /* -----------------------------------------------------------------------------------*/
238: static PetscErrorCode PCFactorSetUseInPlace_Cholesky(PC pc,PetscBool flg)
239: {
240: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
243: dir->inplace = flg;
244: return(0);
245: }
249: static PetscErrorCode PCFactorGetUseInPlace_Cholesky(PC pc,PetscBool *flg)
250: {
251: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
254: *flg = dir->inplace;
255: return(0);
256: }
258: /* -----------------------------------------------------------------------------------*/
262: /*@
263: PCFactorSetReuseOrdering - When similar matrices are factored, this
264: causes the ordering computed in the first factor to be used for all
265: following factors.
267: Logically Collective on PC
269: Input Parameters:
270: + pc - the preconditioner context
271: - flag - PETSC_TRUE to reuse else PETSC_FALSE
273: Options Database Key:
274: . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
276: Level: intermediate
278: .keywords: PC, levels, reordering, factorization, incomplete, LU
280: .seealso: PCFactorSetReuseFill()
281: @*/
282: PetscErrorCode PCFactorSetReuseOrdering(PC pc,PetscBool flag)
283: {
289: PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));
290: return(0);
291: }
293: /*MC
294: PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner
296: Options Database Keys:
297: + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
298: . -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like superlu
299: . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
300: . -pc_factor_fill <fill> - Sets fill amount
301: . -pc_factor_in_place - Activates in-place factorization
302: - -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
304: Notes: Not all options work for all matrix formats
306: Level: beginner
308: Concepts: Cholesky factorization, direct solver
310: Notes: Usually this will compute an "exact" solution in one iteration and does
311: not need a Krylov method (i.e. you can use -ksp_type preonly, or
312: KSPSetType(ksp,KSPPREONLY) for the Krylov method
314: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
315: PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
316: PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount()
317: PCFactorSetUseInPlace(), PCFactorGetUseInPlace(), PCFactorSetMatOrderingType()
319: M*/
323: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
324: {
326: PC_Cholesky *dir;
329: PetscNewLog(pc,&dir);
331: ((PC_Factor*)dir)->fact = 0;
332: dir->inplace = PETSC_FALSE;
334: MatFactorInfoInitialize(&((PC_Factor*)dir)->info);
336: ((PC_Factor*)dir)->factortype = MAT_FACTOR_CHOLESKY;
337: ((PC_Factor*)dir)->info.fill = 5.0;
338: ((PC_Factor*)dir)->info.shifttype = (PetscReal) MAT_SHIFT_NONE;
339: ((PC_Factor*)dir)->info.shiftamount = 0.0;
340: ((PC_Factor*)dir)->info.zeropivot = 100.0*PETSC_MACHINE_EPSILON;
341: ((PC_Factor*)dir)->info.pivotinblocks = 1.0;
343: dir->col = 0;
344: dir->row = 0;
346: PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)dir)->ordering);
347: dir->reusefill = PETSC_FALSE;
348: dir->reuseordering = PETSC_FALSE;
349: pc->data = (void*)dir;
351: pc->ops->destroy = PCDestroy_Cholesky;
352: pc->ops->reset = PCReset_Cholesky;
353: pc->ops->apply = PCApply_Cholesky;
354: pc->ops->applytranspose = PCApplyTranspose_Cholesky;
355: pc->ops->setup = PCSetUp_Cholesky;
356: pc->ops->setfromoptions = PCSetFromOptions_Cholesky;
357: pc->ops->view = PCView_Cholesky;
358: pc->ops->applyrichardson = 0;
359: pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor;
361: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C",PCFactorSetUpMatSolverPackage_Factor);
362: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatSolverPackage_C",PCFactorSetMatSolverPackage_Factor);
363: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetMatSolverPackage_C",PCFactorGetMatSolverPackage_Factor);
364: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetZeroPivot_C",PCFactorSetZeroPivot_Factor);
365: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Factor);
366: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftAmount_C",PCFactorSetShiftAmount_Factor);
367: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetFill_C",PCFactorSetFill_Factor);
368: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUseInPlace_C",PCFactorSetUseInPlace_Cholesky);
369: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetUseInPlace_C",PCFactorGetUseInPlace_Cholesky);
370: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatOrderingType_C",PCFactorSetMatOrderingType_Factor);
371: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C",PCFactorSetReuseOrdering_Cholesky);
372: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseFill_C",PCFactorSetReuseFill_Cholesky);
373: return(0);
374: }