Actual source code: cholesky.c
1: #define PETSCKSP_DLL
3: /*
4: Defines a direct factorization preconditioner for any Mat implementation
5: Note: this need not be consided a preconditioner since it supplies
6: a direct solver.
7: */
8: #include private/pcimpl.h
10: typedef struct {
11: Mat fact; /* factored matrix */
12: PetscReal actualfill; /* actual fill in factor */
13: PetscTruth inplace; /* flag indicating in-place factorization */
14: IS row,col; /* index sets used for reordering */
15: MatOrderingType ordering; /* matrix ordering */
16: PetscTruth reuseordering; /* reuses previous reordering computed */
17: PetscTruth reusefill; /* reuse fill from previous Cholesky */
18: MatFactorInfo info;
19: } PC_Cholesky;
24: PetscErrorCode PCFactorSetZeroPivot_Cholesky(PC pc,PetscReal z)
25: {
26: PC_Cholesky *ch;
29: ch = (PC_Cholesky*)pc->data;
30: ch->info.zeropivot = z;
31: return(0);
32: }
38: PetscErrorCode PCFactorSetShiftNonzero_Cholesky(PC pc,PetscReal shift)
39: {
40: PC_Cholesky *dir;
43: dir = (PC_Cholesky*)pc->data;
44: if (shift == (PetscReal) PETSC_DECIDE) {
45: dir->info.shiftnz = 1.e-12;
46: } else {
47: dir->info.shiftnz = shift;
48: }
49: return(0);
50: }
56: PetscErrorCode PCFactorSetShiftPd_Cholesky(PC pc,PetscTruth shift)
57: {
58: PC_Cholesky *dir;
59:
61: dir = (PC_Cholesky*)pc->data;
62: if (shift) {
63: dir->info.shift_fraction = 0.0;
64: dir->info.shiftpd = 1.0;
65: } else {
66: dir->info.shiftpd = 0.0;
67: }
68: return(0);
69: }
75: PetscErrorCode PCFactorSetReuseOrdering_Cholesky(PC pc,PetscTruth flag)
76: {
77: PC_Cholesky *lu;
78:
80: lu = (PC_Cholesky*)pc->data;
81: lu->reuseordering = flag;
82: return(0);
83: }
89: PetscErrorCode PCFactorSetReuseFill_Cholesky(PC pc,PetscTruth flag)
90: {
91: PC_Cholesky *lu;
92:
94: lu = (PC_Cholesky*)pc->data;
95: lu->reusefill = flag;
96: return(0);
97: }
102: static PetscErrorCode PCSetFromOptions_Cholesky(PC pc)
103: {
104: PC_Cholesky *lu = (PC_Cholesky*)pc->data;
106: PetscTruth flg;
107: char tname[256];
108: PetscFList ordlist;
109:
111: MatOrderingRegisterAll(PETSC_NULL);
112: PetscOptionsHead("Cholesky options");
113: PetscOptionsName("-pc_factor_in_place","Form Cholesky in the same memory as the matrix","PCFactorSetUseInPlace",&flg);
114: if (flg) {
115: PCFactorSetUseInPlace(pc);
116: }
117: PetscOptionsReal("-pc_factor_fill","Expected non-zeros in Cholesky/non-zeros in matrix","PCFactorSetFill",lu->info.fill,&lu->info.fill,0);
118:
119: PetscOptionsName("-pc_factor_reuse_fill","Use fill from previous factorization","PCFactorSetReuseFill",&flg);
120: if (flg) {
121: PCFactorSetReuseFill(pc,PETSC_TRUE);
122: }
123: PetscOptionsName("-pc_factor_reuse_ordering","Reuse ordering from previous factorization","PCFactorSetReuseOrdering",&flg);
124: if (flg) {
125: PCFactorSetReuseOrdering(pc,PETSC_TRUE);
126: }
127:
128: MatGetOrderingList(&ordlist);
129: PetscOptionsList("-pc_factor_mat_ordering_type","Reordering to reduce nonzeros in Cholesky","PCFactorSetMatOrdering",ordlist,lu->ordering,tname,256,&flg);
130: if (flg) {
131: PCFactorSetMatOrdering(pc,tname);
132: }
133: PetscOptionsName("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",&flg);
134: if (flg) {
135: PCFactorSetShiftNonzero(pc,(PetscReal) PETSC_DECIDE);
136: }
137: PetscOptionsReal("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",lu->info.shiftnz,&lu->info.shiftnz,0);
138: PetscOptionsName("-pc_factor_shift_positive_definite","Manteuffel shift applied to diagonal","PCFactorSetShiftPd",&flg);
139: if (flg) {
140: PCFactorSetShiftPd(pc,PETSC_TRUE);
141: }
142: PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",lu->info.zeropivot,&lu->info.zeropivot,0);
144: PetscOptionsTail();
145: return(0);
146: }
150: static PetscErrorCode PCView_Cholesky(PC pc,PetscViewer viewer)
151: {
152: PC_Cholesky *lu = (PC_Cholesky*)pc->data;
154: PetscTruth iascii,isstring;
155:
157: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
158: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
159: if (iascii) {
160:
161: if (lu->inplace) {PetscViewerASCIIPrintf(viewer," Cholesky: in-place factorization\n");}
162: else {PetscViewerASCIIPrintf(viewer," Cholesky: out-of-place factorization\n");}
163: PetscViewerASCIIPrintf(viewer," matrix ordering: %s\n",lu->ordering);
164: if (lu->fact) {
165: PetscViewerASCIIPrintf(viewer," Cholesky: factor fill ratio needed %G\n",lu->actualfill);
166: PetscViewerASCIIPrintf(viewer," Factored matrix follows\n");
167: PetscViewerASCIIPushTab(viewer);
168: PetscViewerASCIIPushTab(viewer);
169: PetscViewerASCIIPushTab(viewer);
170: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
171: MatView(lu->fact,viewer);
172: PetscViewerPopFormat(viewer);
173: PetscViewerASCIIPopTab(viewer);
174: PetscViewerASCIIPopTab(viewer);
175: PetscViewerASCIIPopTab(viewer);
176: }
177: if (lu->reusefill) {PetscViewerASCIIPrintf(viewer," Reusing fill from past factorization\n");}
178: if (lu->reuseordering) {PetscViewerASCIIPrintf(viewer," Reusing reordering from past factorization\n");}
179: } else if (isstring) {
180: PetscViewerStringSPrintf(viewer," order=%s",lu->ordering);
181: } else {
182: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCCHOLESKY",((PetscObject)viewer)->type_name);
183: }
184: return(0);
185: }
189: static PetscErrorCode PCGetFactoredMatrix_Cholesky(PC pc,Mat *mat)
190: {
191: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
192:
194: if (!dir->fact) SETERRQ(PETSC_ERR_ORDER,"Matrix not yet factored; call after KSPSetUp() or PCSetUp()");
195: *mat = dir->fact;
196: return(0);
197: }
201: static PetscErrorCode PCSetUp_Cholesky(PC pc)
202: {
204: PetscTruth flg;
205: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
208: if (dir->reusefill && pc->setupcalled) dir->info.fill = dir->actualfill;
209:
210: if (dir->inplace) {
211: if (dir->row && dir->col && (dir->row != dir->col)) {
212: ISDestroy(dir->row);
213: dir->row = 0;
214: }
215: if (dir->col) {
216: ISDestroy(dir->col);
217: dir->col = 0;
218: }
219: MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);
220: if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */
221: ISDestroy(dir->col);
222: dir->col=0;
223: }
224: if (dir->row) {PetscLogObjectParent(pc,dir->row);}
225: MatCholeskyFactor(pc->pmat,dir->row,&dir->info);
226: dir->fact = pc->pmat;
227: } else {
228: MatInfo info;
229: if (!pc->setupcalled) {
230: MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);
231: /* check if dir->row == dir->col */
232: ISEqual(dir->row,dir->col,&flg);
233: if (!flg) SETERRQ(PETSC_ERR_ARG_INCOMP,"row and column permutations must equal");
234: ISDestroy(dir->col); /* only pass one ordering into CholeskyFactor */
235: dir->col=0;
237: PetscOptionsHasName(pc->prefix,"-pc_factor_nonzeros_along_diagonal",&flg);
238: if (flg) {
239: PetscReal tol = 1.e-10;
240: PetscOptionsGetReal(pc->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,PETSC_NULL);
241: MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
242: }
243: if (dir->row) {PetscLogObjectParent(pc,dir->row);}
244: MatCholeskyFactorSymbolic(pc->pmat,dir->row,&dir->info,&dir->fact);
245: MatGetInfo(dir->fact,MAT_LOCAL,&info);
246: dir->actualfill = info.fill_ratio_needed;
247: PetscLogObjectParent(pc,dir->fact);
248: } else if (pc->flag != SAME_NONZERO_PATTERN) {
249: if (!dir->reuseordering) {
250: if (dir->row && dir->col && (dir->row != dir->col)) {
251: ISDestroy(dir->row);
252: dir->row = 0;
253: }
254: if (dir->col) {
255: ISDestroy(dir->col);
256: dir->col =0;
257: }
258: MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);
259: if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */
260: ISDestroy(dir->col);
261: dir->col=0;
262: }
263: PetscOptionsHasName(pc->prefix,"-pc_factor_nonzeros_along_diagonal",&flg);
264: if (flg) {
265: PetscReal tol = 1.e-10;
266: PetscOptionsGetReal(pc->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,PETSC_NULL);
267: MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
268: }
269: if (dir->row) {PetscLogObjectParent(pc,dir->row);}
270: }
271: MatDestroy(dir->fact);
272: MatCholeskyFactorSymbolic(pc->pmat,dir->row,&dir->info,&dir->fact);
273: MatGetInfo(dir->fact,MAT_LOCAL,&info);
274: dir->actualfill = info.fill_ratio_needed;
275: PetscLogObjectParent(pc,dir->fact);
276: }
277: MatCholeskyFactorNumeric(pc->pmat,&dir->info,&dir->fact);
278: }
279: return(0);
280: }
284: static PetscErrorCode PCDestroy_Cholesky(PC pc)
285: {
286: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
290: if (!dir->inplace && dir->fact) {MatDestroy(dir->fact);}
291: if (dir->row) {ISDestroy(dir->row);}
292: if (dir->col) {ISDestroy(dir->col);}
293: PetscStrfree(dir->ordering);
294: PetscFree(dir);
295: return(0);
296: }
300: static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y)
301: {
302: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
304:
306: if (dir->inplace) {MatSolve(pc->pmat,x,y);}
307: else {MatSolve(dir->fact,x,y);}
308: return(0);
309: }
313: static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y)
314: {
315: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
319: if (dir->inplace) {MatSolveTranspose(pc->pmat,x,y);}
320: else {MatSolveTranspose(dir->fact,x,y);}
321: return(0);
322: }
324: /* -----------------------------------------------------------------------------------*/
329: PetscErrorCode PCFactorSetFill_Cholesky(PC pc,PetscReal fill)
330: {
331: PC_Cholesky *dir;
332:
334: dir = (PC_Cholesky*)pc->data;
335: dir->info.fill = fill;
336: return(0);
337: }
343: PetscErrorCode PCFactorSetUseInPlace_Cholesky(PC pc)
344: {
345: PC_Cholesky *dir;
348: dir = (PC_Cholesky*)pc->data;
349: dir->inplace = PETSC_TRUE;
350: return(0);
351: }
357: PetscErrorCode PCFactorSetMatOrdering_Cholesky(PC pc,MatOrderingType ordering)
358: {
359: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
363: PetscStrfree(dir->ordering);
364: PetscStrallocpy(ordering,&dir->ordering);
365: return(0);
366: }
369: /* -----------------------------------------------------------------------------------*/
373: /*@
374: PCFactorSetReuseOrdering - When similar matrices are factored, this
375: causes the ordering computed in the first factor to be used for all
376: following factors.
378: Collective on PC
380: Input Parameters:
381: + pc - the preconditioner context
382: - flag - PETSC_TRUE to reuse else PETSC_FALSE
384: Options Database Key:
385: . -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
387: Level: intermediate
389: .keywords: PC, levels, reordering, factorization, incomplete, LU
391: .seealso: PCFactorSetReuseFill()
392: @*/
393: PetscErrorCode PCFactorSetReuseOrdering(PC pc,PetscTruth flag)
394: {
395: PetscErrorCode ierr,(*f)(PC,PetscTruth);
399: PetscObjectQueryFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C",(void (**)(void))&f);
400: if (f) {
401: (*f)(pc,flag);
402: }
403: return(0);
404: }
408: /*@
409: PCFactorSetReuseFill - When matrices with same different nonzero structure are factored,
410: this causes later ones to use the fill ratio computed in the initial factorization.
412: Collective on PC
414: Input Parameters:
415: + pc - the preconditioner context
416: - flag - PETSC_TRUE to reuse else PETSC_FALSE
418: Options Database Key:
419: . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
421: Level: intermediate
423: .keywords: PC, levels, reordering, factorization, incomplete, Cholesky
425: .seealso: PCFactorSetReuseOrdering()
426: @*/
427: PetscErrorCode PCFactorSetReuseFill(PC pc,PetscTruth flag)
428: {
429: PetscErrorCode ierr,(*f)(PC,PetscTruth);
433: PetscObjectQueryFunction((PetscObject)pc,"PCFactorSetReuseFill_C",(void (**)(void))&f);
434: if (f) {
435: (*f)(pc,flag);
436: }
437: return(0);
438: }
440: /*MC
441: PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner
443: Options Database Keys:
444: + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
445: . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
446: . -pc_factor_fill <fill> - Sets fill amount
447: . -pc_factor_in_place - Activates in-place factorization
448: . -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
449: . -pc_factor_shift_nonzero <shift> - Sets shift amount or PETSC_DECIDE for the default
450: - -pc_factor_shift_positive_definite [PETSC_TRUE/PETSC_FALSE] - Activate/Deactivate PCFactorSetShiftPd(); the value
451: is optional with PETSC_TRUE being the default
453: Notes: Not all options work for all matrix formats
455: Level: beginner
457: Concepts: Cholesky factorization, direct solver
459: Notes: Usually this will compute an "exact" solution in one iteration and does
460: not need a Krylov method (i.e. you can use -ksp_type preonly, or
461: KSPSetType(ksp,KSPPREONLY) for the Krylov method
463: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
464: PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCGetFactoredMatrix(),
465: PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftPd(),
466: PCFactorSetUseInPlace(), PCFactorSetMatOrdering(),PCFactorSetShiftNonzero(),PCFactorSetShiftPd()
468: M*/
473: PetscErrorCode PCCreate_Cholesky(PC pc)
474: {
476: PC_Cholesky *dir;
479: PetscNew(PC_Cholesky,&dir);
480: PetscLogObjectMemory(pc,sizeof(PC_Cholesky));
482: dir->fact = 0;
483: dir->inplace = PETSC_FALSE;
484: MatFactorInfoInitialize(&dir->info);
485: dir->info.fill = 5.0;
486: dir->info.shiftnz = 0.0;
487: dir->info.shiftpd = 0.0; /* false */
488: dir->info.shift_fraction = 0.0;
489: dir->info.pivotinblocks = 1.0;
490: dir->col = 0;
491: dir->row = 0;
492: PetscStrallocpy(MATORDERING_NATURAL,&dir->ordering);
493: dir->reusefill = PETSC_FALSE;
494: dir->reuseordering = PETSC_FALSE;
495: pc->data = (void*)dir;
497: pc->ops->destroy = PCDestroy_Cholesky;
498: pc->ops->apply = PCApply_Cholesky;
499: pc->ops->applytranspose = PCApplyTranspose_Cholesky;
500: pc->ops->setup = PCSetUp_Cholesky;
501: pc->ops->setfromoptions = PCSetFromOptions_Cholesky;
502: pc->ops->view = PCView_Cholesky;
503: pc->ops->applyrichardson = 0;
504: pc->ops->getfactoredmatrix = PCGetFactoredMatrix_Cholesky;
506: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Cholesky",
507: PCFactorSetZeroPivot_Cholesky);
508: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftNonzero_C","PCFactorSetShiftNonzero_Cholesky",
509: PCFactorSetShiftNonzero_Cholesky);
510: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftPd_C","PCFactorSetShiftPd_Cholesky",
511: PCFactorSetShiftPd_Cholesky);
513: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Cholesky",
514: PCFactorSetFill_Cholesky);
515: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUseInPlace_C","PCFactorSetUseInPlace_Cholesky",
516: PCFactorSetUseInPlace_Cholesky);
517: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrdering_C","PCFactorSetMatOrdering_Cholesky",
518: PCFactorSetMatOrdering_Cholesky);
519: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseOrdering_C","PCFactorSetReuseOrdering_Cholesky",
520: PCFactorSetReuseOrdering_Cholesky);
521: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseFill_C","PCFactorSetReuseFill_Cholesky",
522: PCFactorSetReuseFill_Cholesky);
523: return(0);
524: }