Actual source code: cholesky.c
1: /*$Id: cholesky.c,v 1.12 2001/08/06 21:16:36 bsmith Exp $*/
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/sles/pc/pcimpl.h
9: typedef struct {
10: Mat fact; /* factored matrix */
11: PetscReal actualfill; /* actual fill in factor */
12: int inplace; /* flag indicating in-place factorization */
13: IS row,col; /* index sets used for reordering */
14: MatOrderingType ordering; /* matrix ordering */
15: PetscTruth reuseordering; /* reuses previous reordering computed */
16: PetscTruth reusefill; /* reuse fill from previous Cholesky */
17: MatCholeskyInfo info;
18: } PC_Cholesky;
20: EXTERN_C_BEGIN
21: #undef __FUNCT__
23: int PCCholeskySetReuseOrdering_Cholesky(PC pc,PetscTruth flag)
24: {
25: PC_Cholesky *lu;
26:
28: lu = (PC_Cholesky*)pc->data;
29: lu->reuseordering = flag;
30: return(0);
31: }
32: EXTERN_C_END
34: EXTERN_C_BEGIN
35: #undef __FUNCT__
37: int PCCholeskySetReuseFill_Cholesky(PC pc,PetscTruth flag)
38: {
39: PC_Cholesky *lu;
40:
42: lu = (PC_Cholesky*)pc->data;
43: lu->reusefill = flag;
44: return(0);
45: }
46: EXTERN_C_END
48: #undef __FUNCT__
50: static int PCSetFromOptions_Cholesky(PC pc)
51: {
52: PC_Cholesky *lu = (PC_Cholesky*)pc->data;
53: int ierr;
54: PetscTruth flg;
55: char tname[256];
56: PetscFList ordlist;
57:
59: MatOrderingRegisterAll(PETSC_NULL);
60: PetscOptionsHead("Cholesky options");
61: PetscOptionsName("-pc_cholesky_in_place","Form Cholesky in the same memory as the matrix","PCCholeskySetUseInPlace",&flg);
62: if (flg) {
63: PCCholeskySetUseInPlace(pc);
64: }
65: PetscOptionsReal("-pc_cholesky_fill","Expected non-zeros in Cholesky/non-zeros in matrix","PCCholeskySetFill",lu->info.fill,&lu->info.fill,0);
66:
67: PetscOptionsName("-pc_cholesky_reuse_fill","Use fill from previous factorization","PCCholeskySetReuseFill",&flg);
68: if (flg) {
69: PCCholeskySetReuseFill(pc,PETSC_TRUE);
70: }
71: PetscOptionsName("-pc_cholesky_reuse_ordering","Reuse ordering from previous factorization","PCCholeskySetReuseOrdering",&flg);
72: if (flg) {
73: PCCholeskySetReuseOrdering(pc,PETSC_TRUE);
74: }
75:
76: MatGetOrderingList(&ordlist);
77: PetscOptionsList("-pc_cholesky_mat_ordering_type","Reordering to reduce nonzeros in Cholesky","PCCholeskySetMatOrdering",ordlist,lu->ordering,tname,256,&flg);
78: if (flg) {
79: PCCholeskySetMatOrdering(pc,tname);
80: }
81:
82: PetscOptionsTail();
83: return(0);
84: }
86: #undef __FUNCT__
88: static int PCView_Cholesky(PC pc,PetscViewer viewer)
89: {
90: PC_Cholesky *lu = (PC_Cholesky*)pc->data;
91: int ierr;
92: PetscTruth isascii,isstring;
93:
95: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
96: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
97: if (isascii) {
98: MatInfo info;
99:
100: if (lu->inplace) {PetscViewerASCIIPrintf(viewer," Cholesky: in-place factorizationn");}
101: else {PetscViewerASCIIPrintf(viewer," Cholesky: out-of-place factorizationn");}
102: PetscViewerASCIIPrintf(viewer," matrix ordering: %sn",lu->ordering);
103: if (lu->fact) {
104: MatGetInfo(lu->fact,MAT_LOCAL,&info);
105: PetscViewerASCIIPrintf(viewer," Cholesky nonzeros %gn",info.nz_used);
106: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_FACTOR_INFO);
107: MatView(lu->fact,viewer);
108: PetscViewerPopFormat(viewer);
109: }
110: if (lu->reusefill) {PetscViewerASCIIPrintf(viewer," Reusing fill from past factorizationn");}
111: if (lu->reuseordering) {PetscViewerASCIIPrintf(viewer," Reusing reordering from past factorizationn");}
112: } else if (isstring) {
113: PetscViewerStringSPrintf(viewer," order=%s",lu->ordering);
114: } else {
115: SETERRQ1(1,"Viewer type %s not supported for PCCholesky",((PetscObject)viewer)->type_name);
116: }
117: return(0);
118: }
120: #undef __FUNCT__
122: static int PCGetFactoredMatrix_Cholesky(PC pc,Mat *mat)
123: {
124: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
125:
127: if (!dir->fact) SETERRQ(1,"Matrix not yet factored; call after SLESSetUp() or PCSetUp()");
128: *mat = dir->fact;
129: return(0);
130: }
132: #undef __FUNCT__
134: static int PCSetUp_Cholesky(PC pc)
135: {
136: int ierr;
137: PetscTruth flg;
138: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
141: if (dir->reusefill && pc->setupcalled) dir->info.fill = dir->actualfill;
142:
143: if (dir->inplace) {
144: if (dir->row && dir->col && (dir->row != dir->col)) {
145: ISDestroy(dir->row);
146: dir->row = 0;
147: }
148: if (dir->col) {
149: ISDestroy(dir->col);
150: dir->col = 0;
151: }
152: MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);
153: if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */
154: ISDestroy(dir->col);
155: dir->col=0;
156: }
157: if (dir->row) {PetscLogObjectParent(pc,dir->row);}
158: MatCholeskyFactor(pc->pmat,dir->row,dir->info.fill);
159: dir->fact = pc->pmat;
160: } else {
161: MatInfo info;
162: if (!pc->setupcalled) {
163: MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);
164: if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */
165: ISDestroy(dir->col);
166: dir->col=0;
167: }
168: PetscOptionsHasName(pc->prefix,"-pc_cholesky_nonzeros_along_diagonal",&flg);
169: if (flg) {
170: PetscReal tol = 1.e-10;
171: PetscOptionsGetReal(pc->prefix,"-pc_cholesky_nonzeros_along_diagonal",&tol,PETSC_NULL);
172: MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
173: }
174: if (dir->row) {PetscLogObjectParent(pc,dir->row);}
175: MatCholeskyFactorSymbolic(pc->pmat,dir->row,dir->info.fill,&dir->fact);
176: MatGetInfo(dir->fact,MAT_LOCAL,&info);
177: dir->actualfill = info.fill_ratio_needed;
178: PetscLogObjectParent(pc,dir->fact);
179: } else if (pc->flag != SAME_NONZERO_PATTERN) {
180: if (!dir->reuseordering) {
181: if (dir->row && dir->col && (dir->row != dir->col)) {
182: ISDestroy(dir->row);
183: dir->row = 0;
184: }
185: if (dir->col) {
186: ISDestroy(dir->col);
187: dir->col =0;
188: }
189: MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);
190: if (dir->col && (dir->row != dir->col)) { /* only use row ordering for SBAIJ */
191: ISDestroy(dir->col);
192: dir->col=0;
193: }
194: PetscOptionsHasName(pc->prefix,"-pc_cholesky_nonzeros_along_diagonal",&flg);
195: if (flg) {
196: PetscReal tol = 1.e-10;
197: PetscOptionsGetReal(pc->prefix,"-pc_cholesky_nonzeros_along_diagonal",&tol,PETSC_NULL);
198: MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
199: }
200: if (dir->row) {PetscLogObjectParent(pc,dir->row);}
201: }
202: MatDestroy(dir->fact);
203: MatCholeskyFactorSymbolic(pc->pmat,dir->row,dir->info.fill,&dir->fact);
204: MatGetInfo(dir->fact,MAT_LOCAL,&info);
205: dir->actualfill = info.fill_ratio_needed;
206: PetscLogObjectParent(pc,dir->fact);
207: }
208: MatCholeskyFactorNumeric(pc->pmat,&dir->fact);
209: }
210: return(0);
211: }
213: #undef __FUNCT__
215: static int PCDestroy_Cholesky(PC pc)
216: {
217: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
218: int ierr;
221: if (!dir->inplace && dir->fact) {MatDestroy(dir->fact);}
222: if (dir->row) {ISDestroy(dir->row);}
223: if (dir->col) {ISDestroy(dir->col);}
224: PetscStrfree(dir->ordering);
225: PetscFree(dir);
226: return(0);
227: }
229: #undef __FUNCT__
231: static int PCApply_Cholesky(PC pc,Vec x,Vec y)
232: {
233: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
234: int ierr;
235:
237: if (dir->inplace) {MatSolve(pc->pmat,x,y);}
238: else {MatSolve(dir->fact,x,y);}
239: return(0);
240: }
242: #undef __FUNCT__
244: static int PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y)
245: {
246: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
247: int ierr;
250: if (dir->inplace) {MatSolveTranspose(pc->pmat,x,y);}
251: else {MatSolveTranspose(dir->fact,x,y);}
252: return(0);
253: }
255: /* -----------------------------------------------------------------------------------*/
257: EXTERN_C_BEGIN
258: #undef __FUNCT__
260: int PCCholeskySetFill_Cholesky(PC pc,PetscReal fill)
261: {
262: PC_Cholesky *dir;
263:
265: dir = (PC_Cholesky*)pc->data;
266: dir->info.fill = fill;
267: return(0);
268: }
269: EXTERN_C_END
271: EXTERN_C_BEGIN
272: #undef __FUNCT__
274: int PCCholeskySetDamping_Cholesky(PC pc,PetscReal damping)
275: {
276: PC_Cholesky *dir;
277:
279: dir = (PC_Cholesky*)pc->data;
280: if (damping == (PetscReal) PETSC_DECIDE) {
281: dir->info.damping = 1.e-12;
282: } else {
283: dir->info.damping = damping;
284: }
285: return(0);
286: }
287: EXTERN_C_END
289: EXTERN_C_BEGIN
290: #undef __FUNCT__
292: int PCCholeskySetUseInPlace_Cholesky(PC pc)
293: {
294: PC_Cholesky *dir;
297: dir = (PC_Cholesky*)pc->data;
298: dir->inplace = 1;
299: return(0);
300: }
301: EXTERN_C_END
303: EXTERN_C_BEGIN
304: #undef __FUNCT__
306: int PCCholeskySetMatOrdering_Cholesky(PC pc,MatOrderingType ordering)
307: {
308: PC_Cholesky *dir = (PC_Cholesky*)pc->data;
309: int ierr;
312: PetscStrfree(dir->ordering);
313: PetscStrallocpy(ordering,&dir->ordering);
314: return(0);
315: }
316: EXTERN_C_END
318: /* -----------------------------------------------------------------------------------*/
320: #undef __FUNCT__
322: /*@
323: PCCholeskySetReuseOrdering - When similar matrices are factored, this
324: causes the ordering computed in the first factor to be used for all
325: following factors.
327: Collective on PC
329: Input Parameters:
330: + pc - the preconditioner context
331: - flag - PETSC_TRUE to reuse else PETSC_FALSE
333: Options Database Key:
334: . -pc_cholesky_reuse_ordering - Activate PCCholeskySetReuseOrdering()
336: Level: intermediate
338: .keywords: PC, levels, reordering, factorization, incomplete, LU
340: .seealso: PCCholeskySetReuseFill(), PCICholeskySetReuseOrdering(), PCICholeskyDTSetReuseFill()
341: @*/
342: int PCCholeskySetReuseOrdering(PC pc,PetscTruth flag)
343: {
344: int ierr,(*f)(PC,PetscTruth);
348: PetscObjectQueryFunction((PetscObject)pc,"PCCholeskySetReuseOrdering_C",(void (**)(void))&f);
349: if (f) {
350: (*f)(pc,flag);
351: }
352: return(0);
353: }
355: #undef __FUNCT__
357: /*@
358: PCCholeskySetReuseFill - When matrices with same nonzero structure are Cholesky factored,
359: this causes later ones to use the fill computed in the initial factorization.
361: Collective on PC
363: Input Parameters:
364: + pc - the preconditioner context
365: - flag - PETSC_TRUE to reuse else PETSC_FALSE
367: Options Database Key:
368: . -pc_cholesky_reuse_fill - Activates PCCholeskySetReuseFill()
370: Level: intermediate
372: .keywords: PC, levels, reordering, factorization, incomplete, Cholesky
374: .seealso: PCICholeskySetReuseOrdering(), PCCholeskySetReuseOrdering(), PCICholeskyDTSetReuseFill()
375: @*/
376: int PCCholeskySetReuseFill(PC pc,PetscTruth flag)
377: {
378: int ierr,(*f)(PC,PetscTruth);
382: PetscObjectQueryFunction((PetscObject)pc,"PCCholeskySetReuseFill_C",(void (**)(void))&f);
383: if (f) {
384: (*f)(pc,flag);
385: }
386: return(0);
387: }
389: #undef __FUNCT__
391: /*@
392: PCCholeskySetFill - Indicates the amount of fill you expect in the factored matrix,
393: fill = number nonzeros in factor/number nonzeros in original matrix.
395: Collective on PC
396:
397: Input Parameters:
398: + pc - the preconditioner context
399: - fill - amount of expected fill
401: Options Database Key:
402: . -pc_cholesky_fill <fill> - Sets fill amount
404: Level: intermediate
406: Note:
407: For sparse matrix factorizations it is difficult to predict how much
408: fill to expect. By running with the option -log_info PETSc will print the
409: actual amount of fill used; allowing you to set the value accurately for
410: future runs. Default PETSc uses a value of 5.0
412: .keywords: PC, set, factorization, direct, fill
414: .seealso: PCILUSetFill()
415: @*/
416: int PCCholeskySetFill(PC pc,PetscReal fill)
417: {
418: int ierr,(*f)(PC,PetscReal);
422: if (fill < 1.0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Fill factor cannot be less then 1.0");
423: PetscObjectQueryFunction((PetscObject)pc,"PCCholeskySetFill_C",(void (**)(void))&f);
424: if (f) {
425: (*f)(pc,fill);
426: }
427: return(0);
428: }
430: #undef __FUNCT__
432: /*@
433: PCCholeskySetDamping - Adds this quantity to the diagonal of the matrix during the
434: Cholesky numerical factorization.
436: Collective on PC
437:
438: Input Parameters:
439: + pc - the preconditioner context
440: - damping - amount of damping
442: Options Database Key:
443: . -pc_cholesky_damping <damping> - Sets damping amount
445: Level: intermediate
447: .keywords: PC, set, factorization, direct, fill
449: .seealso: PCICholeskySetFill(), PCILUSetDamp()
450: @*/
451: int PCCholeskySetDamping(PC pc,PetscReal damping)
452: {
453: int ierr,(*f)(PC,PetscReal);
457: PetscObjectQueryFunction((PetscObject)pc,"PCCholeskySetDamping_C",(void (**)(void))&f);
458: if (f) {
459: (*f)(pc,damping);
460: }
461: return(0);
462: }
464: #undef __FUNCT__
466: /*@
467: PCCholeskySetUseInPlace - Tells the system to do an in-place factorization.
468: For dense matrices, this enables the solution of much larger problems.
469: For sparse matrices the factorization cannot be done truly in-place
470: so this does not save memory during the factorization, but after the matrix
471: is factored, the original unfactored matrix is freed, thus recovering that
472: space.
474: Collective on PC
476: Input Parameters:
477: . pc - the preconditioner context
479: Options Database Key:
480: . -pc_cholesky_in_place - Activates in-place factorization
482: Notes:
483: PCCholeskySetUseInplace() can only be used with the KSP method KSPPREONLY or when
484: a different matrix is provided for the multiply and the preconditioner in
485: a call to SLESSetOperators().
486: This is because the Krylov space methods require an application of the
487: matrix multiplication, which is not possible here because the matrix has
488: been factored in-place, replacing the original matrix.
490: Level: intermediate
492: .keywords: PC, set, factorization, direct, inplace, in-place, Cholesky
494: .seealso: PCICholeskySetUseInPlace()
495: @*/
496: int PCCholeskySetUseInPlace(PC pc)
497: {
498: int ierr,(*f)(PC);
502: PetscObjectQueryFunction((PetscObject)pc,"PCCholeskySetUseInPlace_C",(void (**)(void))&f);
503: if (f) {
504: (*f)(pc);
505: }
506: return(0);
507: }
509: #undef __FUNCT__
511: /*@
512: PCCholeskySetMatOrdering - Sets the ordering routine (to reduce fill) to
513: be used it the Cholesky factorization.
515: Collective on PC
517: Input Parameters:
518: + pc - the preconditioner context
519: - ordering - the matrix ordering name, for example, MATORDERING_ND or MATORDERING_RCM
521: Options Database Key:
522: . -pc_cholesky_mat_ordering_type <nd,rcm,...> - Sets ordering routine
524: Level: intermediate
526: .seealso: PCICholeskySetMatOrdering()
527: @*/
528: int PCCholeskySetMatOrdering(PC pc,MatOrderingType ordering)
529: {
530: int ierr,(*f)(PC,MatOrderingType);
533: PetscObjectQueryFunction((PetscObject)pc,"PCCholeskySetMatOrdering_C",(void (**)(void))&f);
534: if (f) {
535: (*f)(pc,ordering);
536: }
537: return(0);
538: }
540: /*MC
541: PCCholesky - Uses a direct solver, based on Cholesky factorization, as a preconditioner
543: Options Database Keys:
544: + -pc_cholesky_reuse_ordering - Activate PCLUSetReuseOrdering()
545: . -pc_cholesky_reuse_fill - Activates PCLUSetReuseFill()
546: . -pc_cholesky_fill <fill> - Sets fill amount
547: . -pc_cholesky_damping <damping> - Sets damping amount
548: . -pc_cholesky_in_place - Activates in-place factorization
549: - -pc_cholesky_mat_ordering_type <nd,rcm,...> - Sets ordering routine
551: Notes: Not all options work for all matrix formats
553: Level: beginner
555: Concepts: Cholesky factorization, direct solver
557: Notes: Usually this will compute an "exact" solution in one iteration and does
558: not need a Krylov method (i.e. you can use -ksp_type preonly, or
559: KSPSetType(ksp,KSPPREONLY) for the Krylov method
561: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
562: PCILU, PCLU, PCICC, PCCholeskySetReuseOrdering(), PCCholeskySetReuseFill(), PCGetFactoredMatrix(),
563: PCCholeskySetFill(), PCCholeskySetDamping(), PCCholeskySetUseInPlace(), PCCholeskySetMatOrdering()
565: M*/
567: EXTERN_C_BEGIN
568: #undef __FUNCT__
570: int PCCreate_Cholesky(PC pc)
571: {
572: int ierr;
573: PC_Cholesky *dir;
576: PetscNew(PC_Cholesky,&dir);
577: PetscLogObjectMemory(pc,sizeof(PC_Cholesky));
579: dir->fact = 0;
580: dir->inplace = 0;
581: dir->info.fill = 5.0;
582: dir->info.damping = 0.0;
583: dir->info.pivotinblocks = 1.0;
584: dir->col = 0;
585: dir->row = 0;
586: PetscStrallocpy(MATORDERING_NATURAL,&dir->ordering);
587: dir->reusefill = PETSC_FALSE;
588: dir->reuseordering = PETSC_FALSE;
589: pc->data = (void*)dir;
591: pc->ops->destroy = PCDestroy_Cholesky;
592: pc->ops->apply = PCApply_Cholesky;
593: pc->ops->applytranspose = PCApplyTranspose_Cholesky;
594: pc->ops->setup = PCSetUp_Cholesky;
595: pc->ops->setfromoptions = PCSetFromOptions_Cholesky;
596: pc->ops->view = PCView_Cholesky;
597: pc->ops->applyrichardson = 0;
598: pc->ops->getfactoredmatrix = PCGetFactoredMatrix_Cholesky;
600: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCholeskySetFill_C","PCCholeskySetFill_Cholesky",
601: PCCholeskySetFill_Cholesky);
602: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCholeskySetDamping_C","PCCholeskySetDamping_Cholesky",
603: PCCholeskySetDamping_Cholesky);
604: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCholeskySetUseInPlace_C","PCCholeskySetUseInPlace_Cholesky",
605: PCCholeskySetUseInPlace_Cholesky);
606: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCholeskySetMatOrdering_C","PCCholeskySetMatOrdering_Cholesky",
607: PCCholeskySetMatOrdering_Cholesky);
608: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCholeskySetReuseOrdering_C","PCCholeskySetReuseOrdering_Cholesky",
609: PCCholeskySetReuseOrdering_Cholesky);
610: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCholeskySetReuseFill_C","PCCholeskySetReuseFill_Cholesky",
611: PCCholeskySetReuseFill_Cholesky);
612: return(0);
613: }
614: EXTERN_C_END