Actual source code: icc.c
2: #include <../src/ksp/pc/impls/factor/icc/icc.h>
4: static PetscErrorCode PCSetUp_ICC(PC pc)
5: {
6: PC_ICC *icc = (PC_ICC*)pc->data;
7: IS perm = NULL,cperm = NULL;
8: PetscErrorCode ierr;
9: MatInfo info;
10: MatSolverType stype;
11: MatFactorError err;
12: PetscBool canuseordering;
15: pc->failedreason = PC_NOERROR;
17: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
18: if (!pc->setupcalled) {
19: if (!((PC_Factor*)icc)->fact) {
20: MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);
21: }
22: MatFactorGetCanUseOrdering(((PC_Factor*)icc)->fact,&canuseordering);
23: if (canuseordering) {
24: PCFactorSetDefaultOrdering_Factor(pc);
25: MatGetOrdering(pc->pmat, ((PC_Factor*)icc)->ordering,&perm,&cperm);
26: }
27: MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);
28: } else if (pc->flag != SAME_NONZERO_PATTERN) {
29: PetscBool canuseordering;
30: MatDestroy(&((PC_Factor*)icc)->fact);
31: MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);
32: MatFactorGetCanUseOrdering(((PC_Factor*)icc)->fact,&canuseordering);
33: if (canuseordering) {
34: PCFactorSetDefaultOrdering_Factor(pc);
35: MatGetOrdering(pc->pmat, ((PC_Factor*)icc)->ordering,&perm,&cperm);
36: }
37: MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);
38: }
39: MatGetInfo(((PC_Factor*)icc)->fact,MAT_LOCAL,&info);
40: icc->hdr.actualfill = info.fill_ratio_needed;
42: ISDestroy(&cperm);
43: ISDestroy(&perm);
45: MatFactorGetError(((PC_Factor*)icc)->fact,&err);
46: if (err) { /* FactorSymbolic() fails */
47: pc->failedreason = (PCFailedReason)err;
48: return(0);
49: }
51: MatCholeskyFactorNumeric(((PC_Factor*)icc)->fact,pc->pmat,&((PC_Factor*)icc)->info);
52: MatFactorGetError(((PC_Factor*)icc)->fact,&err);
53: if (err) { /* FactorNumeric() fails */
54: pc->failedreason = (PCFailedReason)err;
55: }
57: PCFactorGetMatSolverType(pc,&stype);
58: if (!stype) {
59: MatSolverType solverpackage;
60: MatFactorGetSolverType(((PC_Factor*)icc)->fact,&solverpackage);
61: PCFactorSetMatSolverType(pc,solverpackage);
62: }
63: return(0);
64: }
66: static PetscErrorCode PCReset_ICC(PC pc)
67: {
68: PC_ICC *icc = (PC_ICC*)pc->data;
72: MatDestroy(&((PC_Factor*)icc)->fact);
73: return(0);
74: }
76: static PetscErrorCode PCDestroy_ICC(PC pc)
77: {
78: PC_ICC *icc = (PC_ICC*)pc->data;
82: PCReset_ICC(pc);
83: PetscFree(((PC_Factor*)icc)->ordering);
84: PetscFree(((PC_Factor*)icc)->solvertype);
85: PetscFree(pc->data);
86: return(0);
87: }
89: static PetscErrorCode PCApply_ICC(PC pc,Vec x,Vec y)
90: {
91: PC_ICC *icc = (PC_ICC*)pc->data;
95: MatSolve(((PC_Factor*)icc)->fact,x,y);
96: return(0);
97: }
99: static PetscErrorCode PCMatApply_ICC(PC pc,Mat X,Mat Y)
100: {
101: PC_ICC *icc = (PC_ICC*)pc->data;
105: MatMatSolve(((PC_Factor*)icc)->fact,X,Y);
106: return(0);
107: }
109: static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc,Vec x,Vec y)
110: {
112: PC_ICC *icc = (PC_ICC*)pc->data;
115: MatForwardSolve(((PC_Factor*)icc)->fact,x,y);
116: return(0);
117: }
119: static PetscErrorCode PCApplySymmetricRight_ICC(PC pc,Vec x,Vec y)
120: {
122: PC_ICC *icc = (PC_ICC*)pc->data;
125: MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);
126: return(0);
127: }
129: static PetscErrorCode PCSetFromOptions_ICC(PetscOptionItems *PetscOptionsObject,PC pc)
130: {
131: PC_ICC *icc = (PC_ICC*)pc->data;
132: PetscBool flg;
134: /* PetscReal dt[3];*/
137: PetscOptionsHead(PetscOptionsObject,"ICC Options");
138: PCSetFromOptions_Factor(PetscOptionsObject,pc);
140: PetscOptionsReal("-pc_factor_levels","levels of fill","PCFactorSetLevels",((PC_Factor*)icc)->info.levels,&((PC_Factor*)icc)->info.levels,&flg);
141: /*dt[0] = ((PC_Factor*)icc)->info.dt;
142: dt[1] = ((PC_Factor*)icc)->info.dtcol;
143: dt[2] = ((PC_Factor*)icc)->info.dtcount;
144: PetscInt dtmax = 3;
145: PetscOptionsRealArray("-pc_factor_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg);
146: if (flg) {
147: PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);
148: }
149: */
150: PetscOptionsTail();
151: return(0);
152: }
154: extern PetscErrorCode PCFactorSetDropTolerance_ILU(PC,PetscReal,PetscReal,PetscInt);
156: /*MC
157: PCICC - Incomplete Cholesky factorization preconditioners.
159: Options Database Keys:
160: + -pc_factor_levels <k> - number of levels of fill for ICC(k)
161: . -pc_factor_in_place - only for ICC(0) with natural ordering, reuses the space of the matrix for
162: its factorization (overwrites original matrix)
163: . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
164: - -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
166: Level: beginner
168: Notes:
169: Only implemented for some matrix formats. Not implemented in parallel.
171: For BAIJ matrices this implements a point block ICC.
173: The Manteuffel shift is only implemented for matrices with block size 1
175: By default, the Manteuffel is applied (for matrices with block size 1). Call PCFactorSetShiftType(pc,MAT_SHIFT_POSITIVE_DEFINITE);
176: to turn off the shift.
178: References:
179: . 1. - TONY F. CHAN AND HENK A. VAN DER VORST, Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS,
180: Chapter in Parallel Numerical Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
181: Science and Engineering, Kluwer.
183: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
184: PCFactorSetZeroPivot(), PCFactorSetShiftType(), PCFactorSetShiftAmount(),
185: PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(),
186: PCFactorSetLevels()
188: M*/
190: PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc)
191: {
193: PC_ICC *icc;
196: PetscNewLog(pc,&icc);
197: pc->data = (void*)icc;
198: PCFactorInitialize(pc, MAT_FACTOR_ICC);
200: ((PC_Factor*)icc)->info.fill = 1.0;
201: ((PC_Factor*)icc)->info.dtcol = PETSC_DEFAULT;
202: ((PC_Factor*)icc)->info.shifttype = (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE;
204: pc->ops->apply = PCApply_ICC;
205: pc->ops->matapply = PCMatApply_ICC;
206: pc->ops->applytranspose = PCApply_ICC;
207: pc->ops->setup = PCSetUp_ICC;
208: pc->ops->reset = PCReset_ICC;
209: pc->ops->destroy = PCDestroy_ICC;
210: pc->ops->setfromoptions = PCSetFromOptions_ICC;
211: pc->ops->view = PCView_Factor;
212: pc->ops->applysymmetricleft = PCApplySymmetricLeft_ICC;
213: pc->ops->applysymmetricright = PCApplySymmetricRight_ICC;
214: return(0);
215: }