Actual source code: icc.c
petsc-3.7.5 2017-01-01
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;
12: Mat F;
13: const MatSolverPackage stype;
16: pc->failedreason = PC_NOERROR;
17: MatGetOrdering(pc->pmat, ((PC_Factor*)icc)->ordering,&perm,&cperm);
19: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
20: if (!pc->setupcalled) {
21: if (!((PC_Factor*)icc)->fact) {
22: MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);
23: }
24: MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);
25: } else if (pc->flag != SAME_NONZERO_PATTERN) {
26: MatDestroy(&((PC_Factor*)icc)->fact);
27: MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);
28: MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);
29: }
30: MatGetInfo(((PC_Factor*)icc)->fact,MAT_LOCAL,&info);
31: icc->actualfill = info.fill_ratio_needed;
33: ISDestroy(&cperm);
34: ISDestroy(&perm);
36: F = ((PC_Factor*)icc)->fact;
37: if (F->errortype) { /* FactorSymbolic() fails */
38: pc->failedreason = (PCFailedReason)F->errortype;
39: return(0);
40: }
41:
42: MatCholeskyFactorNumeric(((PC_Factor*)icc)->fact,pc->pmat,&((PC_Factor*)icc)->info);
43: if (F->errortype) { /* FactorNumeric() fails */
44: pc->failedreason = (PCFailedReason)F->errortype;
45: }
47: PCFactorGetMatSolverPackage(pc,&stype);
48: if (!stype) {
49: PCFactorSetMatSolverPackage(pc,((PC_Factor*)icc)->fact->solvertype);
50: }
51: return(0);
52: }
56: static PetscErrorCode PCReset_ICC(PC pc)
57: {
58: PC_ICC *icc = (PC_ICC*)pc->data;
62: MatDestroy(&((PC_Factor*)icc)->fact);
63: return(0);
64: }
68: static PetscErrorCode PCDestroy_ICC(PC pc)
69: {
70: PC_ICC *icc = (PC_ICC*)pc->data;
74: PCReset_ICC(pc);
75: PetscFree(((PC_Factor*)icc)->ordering);
76: PetscFree(((PC_Factor*)icc)->solvertype);
77: PetscFree(pc->data);
78: return(0);
79: }
83: static PetscErrorCode PCApply_ICC(PC pc,Vec x,Vec y)
84: {
85: PC_ICC *icc = (PC_ICC*)pc->data;
89: MatSolve(((PC_Factor*)icc)->fact,x,y);
90: return(0);
91: }
95: static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc,Vec x,Vec y)
96: {
98: PC_ICC *icc = (PC_ICC*)pc->data;
101: MatForwardSolve(((PC_Factor*)icc)->fact,x,y);
102: return(0);
103: }
107: static PetscErrorCode PCApplySymmetricRight_ICC(PC pc,Vec x,Vec y)
108: {
110: PC_ICC *icc = (PC_ICC*)pc->data;
113: MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);
114: return(0);
115: }
119: static PetscErrorCode PCSetFromOptions_ICC(PetscOptionItems *PetscOptionsObject,PC pc)
120: {
121: PC_ICC *icc = (PC_ICC*)pc->data;
122: PetscBool flg;
124: /* PetscReal dt[3];*/
127: PetscOptionsHead(PetscOptionsObject,"ICC Options");
128: PCSetFromOptions_Factor(PetscOptionsObject,pc);
130: PetscOptionsReal("-pc_factor_levels","levels of fill","PCFactorSetLevels",((PC_Factor*)icc)->info.levels,&((PC_Factor*)icc)->info.levels,&flg);
131: /*dt[0] = ((PC_Factor*)icc)->info.dt;
132: dt[1] = ((PC_Factor*)icc)->info.dtcol;
133: dt[2] = ((PC_Factor*)icc)->info.dtcount;
134: PetscInt dtmax = 3;
135: PetscOptionsRealArray("-pc_factor_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg);
136: if (flg) {
137: PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);
138: }
139: */
140: PetscOptionsTail();
141: return(0);
142: }
146: static PetscErrorCode PCView_ICC(PC pc,PetscViewer viewer)
147: {
151: PCView_Factor(pc,viewer);
152: return(0);
153: }
155: extern PetscErrorCode PCFactorSetDropTolerance_ILU(PC,PetscReal,PetscReal,PetscInt);
159: PetscErrorCode PCFactorGetUseInPlace_ICC(PC pc,PetscBool *flg)
160: {
162: *flg = PETSC_FALSE;
163: return(0);
164: }
166: /*MC
167: PCICC - Incomplete Cholesky factorization preconditioners.
169: Options Database Keys:
170: + -pc_factor_levels <k> - number of levels of fill for ICC(k)
171: . -pc_factor_in_place - only for ICC(0) with natural ordering, reuses the space of the matrix for
172: its factorization (overwrites original matrix)
173: . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
174: - -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
176: Level: beginner
178: Concepts: incomplete Cholesky factorization
180: Notes: Only implemented for some matrix formats. Not implemented in parallel.
182: For BAIJ matrices this implements a point block ICC.
184: The Manteuffel shift is only implemented for matrices with block size 1
186: By default, the Manteuffel is applied (for matrices with block size 1). Call PCFactorSetShiftType(pc,MAT_SHIFT_POSITIVE_DEFINITE);
187: to turn off the shift.
189: References:
190: . 1. - TONY F. CHAN AND HENK A. VAN DER VORST, Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS,
191: Chapter in Parallel Numerical Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
192: Science and Engineering, Kluwer.
195: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
196: PCFactorSetZeroPivot(), PCFactorSetShiftType(), PCFactorSetShiftAmount(),
197: PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(),
198: PCFactorSetLevels()
200: M*/
204: PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc)
205: {
207: PC_ICC *icc;
210: PetscNewLog(pc,&icc);
212: ((PC_Factor*)icc)->fact = 0;
214: PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)icc)->ordering);
215: MatFactorInfoInitialize(&((PC_Factor*)icc)->info);
217: ((PC_Factor*)icc)->factortype = MAT_FACTOR_ICC;
218: ((PC_Factor*)icc)->info.levels = 0.;
219: ((PC_Factor*)icc)->info.fill = 1.0;
220: icc->implctx = 0;
222: ((PC_Factor*)icc)->info.dtcol = PETSC_DEFAULT;
223: ((PC_Factor*)icc)->info.shifttype = (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE;
224: ((PC_Factor*)icc)->info.shiftamount = 100.0*PETSC_MACHINE_EPSILON;
225: ((PC_Factor*)icc)->info.zeropivot = 100.0*PETSC_MACHINE_EPSILON;
227: pc->data = (void*)icc;
228: pc->ops->apply = PCApply_ICC;
229: pc->ops->applytranspose = PCApply_ICC;
230: pc->ops->setup = PCSetUp_ICC;
231: pc->ops->reset = PCReset_ICC;
232: pc->ops->destroy = PCDestroy_ICC;
233: pc->ops->setfromoptions = PCSetFromOptions_ICC;
234: pc->ops->view = PCView_ICC;
235: pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor;
236: pc->ops->applysymmetricleft = PCApplySymmetricLeft_ICC;
237: pc->ops->applysymmetricright = PCApplySymmetricRight_ICC;
239: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetZeroPivot_C",PCFactorSetZeroPivot_Factor);
240: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetZeroPivot_C",PCFactorGetZeroPivot_Factor);
241: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Factor);
242: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetShiftType_C",PCFactorGetShiftType_Factor);
243: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftAmount_C",PCFactorSetShiftAmount_Factor);
244: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetShiftAmount_C",PCFactorGetShiftAmount_Factor);
245: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C",PCFactorSetUpMatSolverPackage_Factor);
246: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetMatSolverPackage_C",PCFactorGetMatSolverPackage_Factor);
247: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetLevels_C",PCFactorSetLevels_Factor);
248: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetLevels_C",PCFactorGetLevels_Factor);
249: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetFill_C",PCFactorSetFill_Factor);
250: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatOrderingType_C",PCFactorSetMatOrderingType_Factor);
251: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatSolverPackage_C",PCFactorSetMatSolverPackage_Factor);
252: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetDropTolerance_C",PCFactorSetDropTolerance_ILU);
253: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetUseInPlace_C",PCFactorGetUseInPlace_ICC);
254: return(0);
255: }