Actual source code: icc.c
petsc-3.7.3 2016-08-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: MatGetOrdering(pc->pmat, ((PC_Factor*)icc)->ordering,&perm,&cperm);
18: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
19: if (!pc->setupcalled) {
20: if (!((PC_Factor*)icc)->fact) {
21: MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);
22: }
23: MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);
24: } else if (pc->flag != SAME_NONZERO_PATTERN) {
25: MatDestroy(&((PC_Factor*)icc)->fact);
26: MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);
27: MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);
28: }
29: MatGetInfo(((PC_Factor*)icc)->fact,MAT_LOCAL,&info);
30: icc->actualfill = info.fill_ratio_needed;
32: ISDestroy(&cperm);
33: ISDestroy(&perm);
35: F = ((PC_Factor*)icc)->fact;
36: if (F->errortype) { /* FactorSymbolic() fails */
37: pc->failedreason = (PCFailedReason)F->errortype;
38: return(0);
39: }
40:
41: MatCholeskyFactorNumeric(((PC_Factor*)icc)->fact,pc->pmat,&((PC_Factor*)icc)->info);
42: if (F->errortype) { /* FactorNumeric() fails */
43: pc->failedreason = (PCFailedReason)F->errortype;
44: }
46: PCFactorGetMatSolverPackage(pc,&stype);
47: if (!stype) {
48: PCFactorSetMatSolverPackage(pc,((PC_Factor*)icc)->fact->solvertype);
49: }
50: return(0);
51: }
55: static PetscErrorCode PCReset_ICC(PC pc)
56: {
57: PC_ICC *icc = (PC_ICC*)pc->data;
61: MatDestroy(&((PC_Factor*)icc)->fact);
62: return(0);
63: }
67: static PetscErrorCode PCDestroy_ICC(PC pc)
68: {
69: PC_ICC *icc = (PC_ICC*)pc->data;
73: PCReset_ICC(pc);
74: PetscFree(((PC_Factor*)icc)->ordering);
75: PetscFree(((PC_Factor*)icc)->solvertype);
76: PetscFree(pc->data);
77: return(0);
78: }
82: static PetscErrorCode PCApply_ICC(PC pc,Vec x,Vec y)
83: {
84: PC_ICC *icc = (PC_ICC*)pc->data;
88: MatSolve(((PC_Factor*)icc)->fact,x,y);
89: return(0);
90: }
94: static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc,Vec x,Vec y)
95: {
97: PC_ICC *icc = (PC_ICC*)pc->data;
100: MatForwardSolve(((PC_Factor*)icc)->fact,x,y);
101: return(0);
102: }
106: static PetscErrorCode PCApplySymmetricRight_ICC(PC pc,Vec x,Vec y)
107: {
109: PC_ICC *icc = (PC_ICC*)pc->data;
112: MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);
113: return(0);
114: }
118: static PetscErrorCode PCSetFromOptions_ICC(PetscOptionItems *PetscOptionsObject,PC pc)
119: {
120: PC_ICC *icc = (PC_ICC*)pc->data;
121: PetscBool flg;
123: /* PetscReal dt[3];*/
126: PetscOptionsHead(PetscOptionsObject,"ICC Options");
127: PCSetFromOptions_Factor(PetscOptionsObject,pc);
129: PetscOptionsReal("-pc_factor_levels","levels of fill","PCFactorSetLevels",((PC_Factor*)icc)->info.levels,&((PC_Factor*)icc)->info.levels,&flg);
130: /*dt[0] = ((PC_Factor*)icc)->info.dt;
131: dt[1] = ((PC_Factor*)icc)->info.dtcol;
132: dt[2] = ((PC_Factor*)icc)->info.dtcount;
133: PetscInt dtmax = 3;
134: PetscOptionsRealArray("-pc_factor_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg);
135: if (flg) {
136: PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);
137: }
138: */
139: PetscOptionsTail();
140: return(0);
141: }
145: static PetscErrorCode PCView_ICC(PC pc,PetscViewer viewer)
146: {
150: PCView_Factor(pc,viewer);
151: return(0);
152: }
154: extern PetscErrorCode PCFactorSetDropTolerance_ILU(PC,PetscReal,PetscReal,PetscInt);
158: PetscErrorCode PCFactorGetUseInPlace_ICC(PC pc,PetscBool *flg)
159: {
161: *flg = PETSC_FALSE;
162: return(0);
163: }
165: /*MC
166: PCICC - Incomplete Cholesky factorization preconditioners.
168: Options Database Keys:
169: + -pc_factor_levels <k> - number of levels of fill for ICC(k)
170: . -pc_factor_in_place - only for ICC(0) with natural ordering, reuses the space of the matrix for
171: its factorization (overwrites original matrix)
172: . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
173: - -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
175: Level: beginner
177: Concepts: incomplete Cholesky factorization
179: Notes: Only implemented for some matrix formats. Not implemented in parallel.
181: For BAIJ matrices this implements a point block ICC.
183: The Manteuffel shift is only implemented for matrices with block size 1
185: By default, the Manteuffel is applied (for matrices with block size 1). Call PCFactorSetShiftType(pc,MAT_SHIFT_POSITIVE_DEFINITE);
186: to turn off the shift.
188: References:
189: . 1. - TONY F. CHAN AND HENK A. VAN DER VORST, Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS,
190: Chapter in Parallel Numerical Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
191: Science and Engineering, Kluwer.
194: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
195: PCFactorSetZeroPivot(), PCFactorSetShiftType(), PCFactorSetShiftAmount(),
196: PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(),
197: PCFactorSetLevels()
199: M*/
203: PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc)
204: {
206: PC_ICC *icc;
209: PetscNewLog(pc,&icc);
211: ((PC_Factor*)icc)->fact = 0;
213: PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)icc)->ordering);
214: MatFactorInfoInitialize(&((PC_Factor*)icc)->info);
216: ((PC_Factor*)icc)->factortype = MAT_FACTOR_ICC;
217: ((PC_Factor*)icc)->info.levels = 0.;
218: ((PC_Factor*)icc)->info.fill = 1.0;
219: icc->implctx = 0;
221: ((PC_Factor*)icc)->info.dtcol = PETSC_DEFAULT;
222: ((PC_Factor*)icc)->info.shifttype = (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE;
223: ((PC_Factor*)icc)->info.shiftamount = 100.0*PETSC_MACHINE_EPSILON;
224: ((PC_Factor*)icc)->info.zeropivot = 100.0*PETSC_MACHINE_EPSILON;
226: pc->data = (void*)icc;
227: pc->ops->apply = PCApply_ICC;
228: pc->ops->applytranspose = PCApply_ICC;
229: pc->ops->setup = PCSetUp_ICC;
230: pc->ops->reset = PCReset_ICC;
231: pc->ops->destroy = PCDestroy_ICC;
232: pc->ops->setfromoptions = PCSetFromOptions_ICC;
233: pc->ops->view = PCView_ICC;
234: pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor;
235: pc->ops->applysymmetricleft = PCApplySymmetricLeft_ICC;
236: pc->ops->applysymmetricright = PCApplySymmetricRight_ICC;
238: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C",PCFactorSetUpMatSolverPackage_Factor);
239: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetMatSolverPackage_C",PCFactorGetMatSolverPackage_Factor);
240: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetZeroPivot_C",PCFactorSetZeroPivot_Factor);
241: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Factor);
242: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftAmount_C",PCFactorSetShiftAmount_Factor);
243: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetLevels_C",PCFactorSetLevels_Factor);
244: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetLevels_C",PCFactorGetLevels_Factor);
245: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetFill_C",PCFactorSetFill_Factor);
246: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatOrderingType_C",PCFactorSetMatOrderingType_Factor);
247: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatSolverPackage_C",PCFactorSetMatSolverPackage_Factor);
248: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetDropTolerance_C",PCFactorSetDropTolerance_ILU);
249: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetUseInPlace_C",PCFactorGetUseInPlace_ICC);
250: return(0);
251: }