Actual source code: icc.c
petsc-3.3-p7 2013-05-11
2: #include <../src/ksp/pc/impls/factor/icc/icc.h> /*I "petscpc.h" I*/
6: static PetscErrorCode PCSetup_ICC(PC pc)
7: {
8: PC_ICC *icc = (PC_ICC*)pc->data;
9: IS perm,cperm;
11: MatInfo info;
14: MatGetOrdering(pc->pmat, ((PC_Factor*)icc)->ordering,&perm,&cperm);
16: if (!pc->setupcalled) {
17: if (!((PC_Factor*)icc)->fact){
18: MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,& ((PC_Factor*)icc)->fact);
19: }
20: MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);
21: } else if (pc->flag != SAME_NONZERO_PATTERN) {
22: MatDestroy(&((PC_Factor*)icc)->fact);
23: MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);
24: MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);
25: }
26: MatGetInfo(((PC_Factor*)icc)->fact,MAT_LOCAL,&info);
27: icc->actualfill = info.fill_ratio_needed;
29: ISDestroy(&cperm);
30: ISDestroy(&perm);
31: MatCholeskyFactorNumeric(((PC_Factor*)icc)->fact,pc->pmat,&((PC_Factor*)icc)->info);
32: return(0);
33: }
37: static PetscErrorCode PCReset_ICC(PC pc)
38: {
39: PC_ICC *icc = (PC_ICC*)pc->data;
43: MatDestroy(&((PC_Factor*)icc)->fact);
44: return(0);
45: }
49: static PetscErrorCode PCDestroy_ICC(PC pc)
50: {
51: PC_ICC *icc = (PC_ICC*)pc->data;
55: PCReset_ICC(pc);
56: PetscFree(((PC_Factor*)icc)->ordering);
57: PetscFree(((PC_Factor*)icc)->solvertype);
58: PetscFree(pc->data);
59: return(0);
60: }
64: static PetscErrorCode PCApply_ICC(PC pc,Vec x,Vec y)
65: {
66: PC_ICC *icc = (PC_ICC*)pc->data;
70: MatSolve(((PC_Factor*)icc)->fact,x,y);
71: return(0);
72: }
76: static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc,Vec x,Vec y)
77: {
79: PC_ICC *icc = (PC_ICC*)pc->data;
82: MatForwardSolve(((PC_Factor*)icc)->fact,x,y);
83: return(0);
84: }
88: static PetscErrorCode PCApplySymmetricRight_ICC(PC pc,Vec x,Vec y)
89: {
91: PC_ICC *icc = (PC_ICC*)pc->data;
94: MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);
95: return(0);
96: }
100: static PetscErrorCode PCSetFromOptions_ICC(PC pc)
101: {
102: PC_ICC *icc = (PC_ICC*)pc->data;
103: PetscBool flg;
104: /* PetscReal dt[3];*/
108: PetscOptionsHead("ICC Options");
109: PCSetFromOptions_Factor(pc);
111: PetscOptionsReal("-pc_factor_levels","levels of fill","PCFactorSetLevels",((PC_Factor*)icc)->info.levels,&((PC_Factor*)icc)->info.levels,&flg);
112: /*dt[0] = ((PC_Factor*)icc)->info.dt;
113: dt[1] = ((PC_Factor*)icc)->info.dtcol;
114: dt[2] = ((PC_Factor*)icc)->info.dtcount;
115: PetscInt dtmax = 3;
116: PetscOptionsRealArray("-pc_factor_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg);
117: if (flg) {
118: PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);
119: }
120: */
121: PetscOptionsTail();
122: return(0);
123: }
127: static PetscErrorCode PCView_ICC(PC pc,PetscViewer viewer)
128: {
130:
132: PCView_Factor(pc,viewer);
133: return(0);
134: }
136: EXTERN_C_BEGIN
137: extern PetscErrorCode PCFactorSetDropTolerance_ILU(PC,PetscReal,PetscReal,PetscInt);
138: EXTERN_C_END
140: /*MC
141: PCICC - Incomplete Cholesky factorization preconditioners.
143: Options Database Keys:
144: + -pc_factor_levels <k> - number of levels of fill for ICC(k)
145: . -pc_factor_in_place - only for ICC(0) with natural ordering, reuses the space of the matrix for
146: its factorization (overwrites original matrix)
147: . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
148: - -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
150: Level: beginner
152: Concepts: incomplete Cholesky factorization
154: Notes: Only implemented for some matrix formats. Not implemented in parallel.
156: For BAIJ matrices this implements a point block ICC.
158: The Manteuffel shift is only implemented for matrices with block size 1
160: By default, the Manteuffel is applied (for matrices with block size 1). Call PCFactorSetShiftType(pc,MAT_SHIFT_POSITIVE_DEFINITE);
161: to turn off the shift.
163: References:
164: Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS, TONY F. CHAN AND HENK A. VAN DER VORST
165: http://igitur-archive.library.uu.nl/math/2001-0621-115821/proc.pdf chapter in Parallel Numerical
166: Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
167: Science and Engineering, Kluwer, pp. 167--202.
170: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
171: PCFactorSetZeroPivot(), PCFactorSetShiftType(), PCFactorSetShiftAmount(),
172: PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(),
173: PCFactorSetLevels()
175: M*/
177: EXTERN_C_BEGIN
180: PetscErrorCode PCCreate_ICC(PC pc)
181: {
183: PC_ICC *icc;
186: PetscNewLog(pc,PC_ICC,&icc);
188: ((PC_Factor*)icc)->fact = 0;
189: PetscStrallocpy(MATORDERINGNATURAL,&((PC_Factor*)icc)->ordering);
190: PetscStrallocpy(MATSOLVERPETSC,&((PC_Factor*)icc)->solvertype);
191: MatFactorInfoInitialize(&((PC_Factor*)icc)->info);
192: ((PC_Factor*)icc)->factortype = MAT_FACTOR_ICC;
193: ((PC_Factor*)icc)->info.levels = 0.;
194: ((PC_Factor*)icc)->info.fill = 1.0;
195: icc->implctx = 0;
197: ((PC_Factor*)icc)->info.dtcol = PETSC_DEFAULT;
198: ((PC_Factor*)icc)->info.shifttype = (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE;
199: ((PC_Factor*)icc)->info.shiftamount = 100.0*PETSC_MACHINE_EPSILON;
200: ((PC_Factor*)icc)->info.zeropivot = 100.0*PETSC_MACHINE_EPSILON;
202: pc->data = (void*)icc;
203: pc->ops->apply = PCApply_ICC;
204: pc->ops->applytranspose = PCApply_ICC;
205: pc->ops->setup = PCSetup_ICC;
206: pc->ops->reset = PCReset_ICC;
207: pc->ops->destroy = PCDestroy_ICC;
208: pc->ops->setfromoptions = PCSetFromOptions_ICC;
209: pc->ops->view = PCView_ICC;
210: pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor;
211: pc->ops->applysymmetricleft = PCApplySymmetricLeft_ICC;
212: pc->ops->applysymmetricright = PCApplySymmetricRight_ICC;
214: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C","PCFactorSetUpMatSolverPackage_Factor",
215: PCFactorSetUpMatSolverPackage_Factor);
216: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorGetMatSolverPackage_C","PCFactorGetMatSolverPackage_Factor",
217: PCFactorGetMatSolverPackage_Factor);
218: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Factor",
219: PCFactorSetZeroPivot_Factor);
220: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftType_C","PCFactorSetShiftType_Factor",
221: PCFactorSetShiftType_Factor);
222: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftAmount_C","PCFactorSetShiftAmount_Factor",
223: PCFactorSetShiftAmount_Factor);
224: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetLevels_C","PCFactorSetLevels_Factor",
225: PCFactorSetLevels_Factor);
226: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Factor",
227: PCFactorSetFill_Factor);
228: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Factor",
229: PCFactorSetMatOrderingType_Factor);
230: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatSolverPackage_C","PCFactorSetMatSolverPackage_Factor",
231: PCFactorSetMatSolverPackage_Factor);
232: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetDropTolerance_C","PCFactorSetDropTolerance_ILU",
233: PCFactorSetDropTolerance_ILU);
234: return(0);
235: }
236: EXTERN_C_END