Actual source code: icc.c
1: #define PETSCKSP_DLL
3: /*
4: Defines a Cholesky factorization preconditioner for any Mat implementation.
5: Presently only provided for MPIRowbs format (i.e. BlockSolve).
6: */
8: #include src/ksp/pc/impls/factor/icc/icc.h
13: PetscErrorCode PCFactorSetZeroPivot_ICC(PC pc,PetscReal z)
14: {
15: PC_ICC *icc;
18: icc = (PC_ICC*)pc->data;
19: icc->info.zeropivot = z;
20: return(0);
21: }
27: PetscErrorCode PCFactorSetShiftNonzero_ICC(PC pc,PetscReal shift)
28: {
29: PC_ICC *dir;
32: dir = (PC_ICC*)pc->data;
33: if (shift == (PetscReal) PETSC_DECIDE) {
34: dir->info.shiftnz = 1.e-12;
35: } else {
36: dir->info.shiftnz = shift;
37: }
38: return(0);
39: }
45: PetscErrorCode PCFactorSetShiftPd_ICC(PC pc,PetscTruth shift)
46: {
47: PC_ICC *dir;
48:
50: dir = (PC_ICC*)pc->data;
51: if (shift) {
52: dir->info.shift_fraction = 0.0;
53: dir->info.shiftpd = 1.0;
54: } else {
55: dir->info.shiftpd = 0.0;
56: }
57: return(0);
58: }
64: PetscErrorCode PCFactorSetMatOrdering_ICC(PC pc,MatOrderingType ordering)
65: {
66: PC_ICC *dir = (PC_ICC*)pc->data;
68:
70: PetscStrfree(dir->ordering);
71: PetscStrallocpy(ordering,&dir->ordering);
72: return(0);
73: }
79: PetscErrorCode PCFactorSetFill_ICC(PC pc,PetscReal fill)
80: {
81: PC_ICC *dir;
84: dir = (PC_ICC*)pc->data;
85: dir->info.fill = fill;
86: return(0);
87: }
93: PetscErrorCode PCFactorSetLevels_ICC(PC pc,PetscInt levels)
94: {
95: PC_ICC *icc;
98: icc = (PC_ICC*)pc->data;
99: icc->info.levels = levels;
100: return(0);
101: }
106: static PetscErrorCode PCSetup_ICC(PC pc)
107: {
108: PC_ICC *icc = (PC_ICC*)pc->data;
109: IS perm,cperm;
111: MatInfo info;
114: MatGetOrdering(pc->pmat,icc->ordering,&perm,&cperm);
116: if (!pc->setupcalled) {
117: MatICCFactorSymbolic(pc->pmat,perm,&icc->info,&icc->fact);
118: } else if (pc->flag != SAME_NONZERO_PATTERN) {
119: MatDestroy(icc->fact);
120: MatICCFactorSymbolic(pc->pmat,perm,&icc->info,&icc->fact);
121: }
122: MatGetInfo(icc->fact,MAT_LOCAL,&info);
123: icc->actualfill = info.fill_ratio_needed;
125: ISDestroy(cperm);
126: ISDestroy(perm);
127: MatCholeskyFactorNumeric(pc->pmat,&icc->info,&icc->fact);
128: return(0);
129: }
133: static PetscErrorCode PCDestroy_ICC(PC pc)
134: {
135: PC_ICC *icc = (PC_ICC*)pc->data;
139: if (icc->fact) {MatDestroy(icc->fact);}
140: PetscStrfree(icc->ordering);
141: PetscFree(icc);
142: return(0);
143: }
147: static PetscErrorCode PCApply_ICC(PC pc,Vec x,Vec y)
148: {
149: PC_ICC *icc = (PC_ICC*)pc->data;
153: MatSolve(icc->fact,x,y);
154: return(0);
155: }
159: static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc,Vec x,Vec y)
160: {
162: PC_ICC *icc = (PC_ICC*)pc->data;
165: MatForwardSolve(icc->fact,x,y);
166: return(0);
167: }
171: static PetscErrorCode PCApplySymmetricRight_ICC(PC pc,Vec x,Vec y)
172: {
174: PC_ICC *icc = (PC_ICC*)pc->data;
177: MatBackwardSolve(icc->fact,x,y);
178: return(0);
179: }
183: static PetscErrorCode PCGetFactoredMatrix_ICC(PC pc,Mat *mat)
184: {
185: PC_ICC *icc = (PC_ICC*)pc->data;
188: *mat = icc->fact;
189: return(0);
190: }
194: static PetscErrorCode PCSetFromOptions_ICC(PC pc)
195: {
196: PC_ICC *icc = (PC_ICC*)pc->data;
197: char tname[256];
198: PetscTruth flg;
200: PetscFList ordlist;
203: MatOrderingRegisterAll(PETSC_NULL);
204: PetscOptionsHead("ICC Options");
205: PetscOptionsReal("-pc_factor_levels","levels of fill","PCFactorSetLevels",icc->info.levels,&icc->info.levels,&flg);
206: PetscOptionsReal("-pc_factor_fill","Expected fill in factorization","PCFactorSetFill",icc->info.fill,&icc->info.fill,&flg);
207: MatGetOrderingList(&ordlist);
208: PetscOptionsList("-pc_factor_mat_ordering_type","Reorder to reduce nonzeros in ICC","PCFactorSetMatOrdering",ordlist,icc->ordering,tname,256,&flg);
209: if (flg) {
210: PCFactorSetMatOrdering(pc,tname);
211: }
212: PetscOptionsName("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",&flg);
213: if (flg) {
214: PCFactorSetShiftNonzero(pc,(PetscReal)PETSC_DECIDE);
215: }
216: PetscOptionsReal("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",icc->info.shiftnz,&icc->info.shiftnz,0);
217: PetscOptionsName("-pc_factor_shift_positive_definite","Manteuffel shift applied to diagonal","PCFactorSetShift",&flg);
218: if (flg) {
219: PCFactorSetShiftPd(pc,PETSC_TRUE);
220: } else {
221: PCFactorSetShiftPd(pc,PETSC_FALSE);
222: }
223: PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",icc->info.zeropivot,&icc->info.zeropivot,0);
224:
225: PetscOptionsTail();
226: return(0);
227: }
231: static PetscErrorCode PCView_ICC(PC pc,PetscViewer viewer)
232: {
233: PC_ICC *icc = (PC_ICC*)pc->data;
235: PetscTruth isstring,iascii;
238: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
239: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
240: if (iascii) {
241: if (icc->info.levels == 1) {
242: PetscViewerASCIIPrintf(viewer," ICC: %D level of fill\n",(PetscInt)icc->info.levels);
243: } else {
244: PetscViewerASCIIPrintf(viewer," ICC: %D levels of fill\n",(PetscInt)icc->info.levels);
245: }
246: PetscViewerASCIIPrintf(viewer," ICC: factor fill ratio allocated %G\n",icc->info.fill);
247: if (icc->info.shiftpd) {PetscViewerASCIIPrintf(viewer," ICC: using Manteuffel shift\n");}
248: if (icc->fact) {
249: PetscViewerASCIIPrintf(viewer," ICC: factor fill ratio needed %G\n",icc->actualfill);
250: PetscViewerASCIIPrintf(viewer," Factored matrix follows\n");
251: PetscViewerASCIIPushTab(viewer);
252: PetscViewerASCIIPushTab(viewer);
253: PetscViewerASCIIPushTab(viewer);
254: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
255: MatView(icc->fact,viewer);
256: PetscViewerPopFormat(viewer);
257: PetscViewerASCIIPopTab(viewer);
258: PetscViewerASCIIPopTab(viewer);
259: PetscViewerASCIIPopTab(viewer);
260: }
261: } else if (isstring) {
262: PetscViewerStringSPrintf(viewer," lvls=%D",(PetscInt)icc->info.levels);
263: } else {
264: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCICC",((PetscObject)viewer)->type_name);
265: }
266: return(0);
267: }
269: /*MC
270: PCICC - Incomplete Cholesky factorization preconditioners.
272: Options Database Keys:
273: + -pc_factor_levels <k> - number of levels of fill for ICC(k)
274: . -pc_factor_in_place - only for ICC(0) with natural ordering, reuses the space of the matrix for
275: its factorization (overwrites original matrix)
276: . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
277: . -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
278: . -pc_factor_shift_nonzero <shift> - Sets shift amount or PETSC_DECIDE for the default
279: - -pc_factor_shift_positive_definite [PETSC_TRUE/PETSC_FALSE] - Activate/Deactivate PCFactorSetShiftPd(); the value
280: is optional with PETSC_TRUE being the default
282: Level: beginner
284: Concepts: incomplete Cholesky factorization
286: Notes: Only implemented for some matrix formats. Not implemented in parallel
288: For BAIJ matrices this implements a point block ICC.
290: The Manteuffel shift is only implemented for matrices with block size 1
292: By default, the Manteuffel is applied (for matrices with block size 1). Call PCFactorSetShiftPd(pc,PETSC_FALSE);
293: to turn off the shift.
296: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
297: PCFactorSetZeroPivot(), PCFactorSetShiftNonzero(), PCFactorSetShiftPd(),
298: PCFactorSetFill(), PCFactorSetMatOrdering(), PCFactorSetReuseOrdering(),
299: PCFactorSetLevels(),PCFactorSetShiftNonzero(),PCFactorSetShiftPd(),
301: M*/
306: PetscErrorCode PCCreate_ICC(PC pc)
307: {
309: PC_ICC *icc;
312: PetscNew(PC_ICC,&icc);
313: PetscLogObjectMemory(pc,sizeof(PC_ICC));
315: icc->fact = 0;
316: PetscStrallocpy(MATORDERING_NATURAL,&icc->ordering);
317: MatFactorInfoInitialize(&icc->info);
318: icc->info.levels = 0;
319: icc->info.fill = 1.0;
320: icc->implctx = 0;
322: icc->info.dtcol = PETSC_DEFAULT;
323: icc->info.shiftnz = 0.0;
324: icc->info.shiftpd = 1.0; /* true */
325: icc->info.shift_fraction = 0.0;
326: icc->info.zeropivot = 1.e-12;
327: pc->data = (void*)icc;
329: pc->ops->apply = PCApply_ICC;
330: pc->ops->setup = PCSetup_ICC;
331: pc->ops->destroy = PCDestroy_ICC;
332: pc->ops->setfromoptions = PCSetFromOptions_ICC;
333: pc->ops->view = PCView_ICC;
334: pc->ops->getfactoredmatrix = PCGetFactoredMatrix_ICC;
335: pc->ops->applysymmetricleft = PCApplySymmetricLeft_ICC;
336: pc->ops->applysymmetricright = PCApplySymmetricRight_ICC;
338: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_ICC",
339: PCFactorSetZeroPivot_ICC);
340: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftNonzero_C","PCFactorSetShiftNonzero_ICC",
341: PCFactorSetShiftNonzero_ICC);
342: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftPd_C","PCFactorSetShiftPd_ICC",
343: PCFactorSetShiftPd_ICC);
345: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetLevels_C","PCFactorSetLevels_ICC",
346: PCFactorSetLevels_ICC);
347: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_ICC",
348: PCFactorSetFill_ICC);
349: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrdering_C","PCFactorSetMatOrdering_ICC",
350: PCFactorSetMatOrdering_ICC);
351: return(0);
352: }