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: MatInfo info;
9: MatSolverType stype;
10: MatFactorError err;
11: PetscBool canuseordering;
13: pc->failedreason = PC_NOERROR;
15: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
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: MatFactorGetCanUseOrdering(((PC_Factor*)icc)->fact,&canuseordering);
21: if (canuseordering) {
22: PCFactorSetDefaultOrdering_Factor(pc);
23: MatGetOrdering(pc->pmat, ((PC_Factor*)icc)->ordering,&perm,&cperm);
24: }
25: MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);
26: } else if (pc->flag != SAME_NONZERO_PATTERN) {
27: PetscBool canuseordering;
28: MatDestroy(&((PC_Factor*)icc)->fact);
29: MatGetFactor(pc->pmat,((PC_Factor*)icc)->solvertype,MAT_FACTOR_ICC,&((PC_Factor*)icc)->fact);
30: MatFactorGetCanUseOrdering(((PC_Factor*)icc)->fact,&canuseordering);
31: if (canuseordering) {
32: PCFactorSetDefaultOrdering_Factor(pc);
33: MatGetOrdering(pc->pmat, ((PC_Factor*)icc)->ordering,&perm,&cperm);
34: }
35: MatICCFactorSymbolic(((PC_Factor*)icc)->fact,pc->pmat,perm,&((PC_Factor*)icc)->info);
36: }
37: MatGetInfo(((PC_Factor*)icc)->fact,MAT_LOCAL,&info);
38: icc->hdr.actualfill = info.fill_ratio_needed;
40: ISDestroy(&cperm);
41: ISDestroy(&perm);
43: MatFactorGetError(((PC_Factor*)icc)->fact,&err);
44: if (err) { /* FactorSymbolic() fails */
45: pc->failedreason = (PCFailedReason)err;
46: return 0;
47: }
49: MatCholeskyFactorNumeric(((PC_Factor*)icc)->fact,pc->pmat,&((PC_Factor*)icc)->info);
50: MatFactorGetError(((PC_Factor*)icc)->fact,&err);
51: if (err) { /* FactorNumeric() fails */
52: pc->failedreason = (PCFailedReason)err;
53: }
55: PCFactorGetMatSolverType(pc,&stype);
56: if (!stype) {
57: MatSolverType solverpackage;
58: MatFactorGetSolverType(((PC_Factor*)icc)->fact,&solverpackage);
59: PCFactorSetMatSolverType(pc,solverpackage);
60: }
61: return 0;
62: }
64: static PetscErrorCode PCReset_ICC(PC pc)
65: {
66: PC_ICC *icc = (PC_ICC*)pc->data;
68: MatDestroy(&((PC_Factor*)icc)->fact);
69: return 0;
70: }
72: static PetscErrorCode PCDestroy_ICC(PC pc)
73: {
74: PC_ICC *icc = (PC_ICC*)pc->data;
76: PCReset_ICC(pc);
77: PetscFree(((PC_Factor*)icc)->ordering);
78: PetscFree(((PC_Factor*)icc)->solvertype);
79: PetscFree(pc->data);
80: return 0;
81: }
83: static PetscErrorCode PCApply_ICC(PC pc,Vec x,Vec y)
84: {
85: PC_ICC *icc = (PC_ICC*)pc->data;
87: MatSolve(((PC_Factor*)icc)->fact,x,y);
88: return 0;
89: }
91: static PetscErrorCode PCMatApply_ICC(PC pc,Mat X,Mat Y)
92: {
93: PC_ICC *icc = (PC_ICC*)pc->data;
95: MatMatSolve(((PC_Factor*)icc)->fact,X,Y);
96: return 0;
97: }
99: static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc,Vec x,Vec y)
100: {
101: PC_ICC *icc = (PC_ICC*)pc->data;
103: MatForwardSolve(((PC_Factor*)icc)->fact,x,y);
104: return 0;
105: }
107: static PetscErrorCode PCApplySymmetricRight_ICC(PC pc,Vec x,Vec y)
108: {
109: PC_ICC *icc = (PC_ICC*)pc->data;
111: MatBackwardSolve(((PC_Factor*)icc)->fact,x,y);
112: return 0;
113: }
115: static PetscErrorCode PCSetFromOptions_ICC(PetscOptionItems *PetscOptionsObject,PC pc)
116: {
117: PC_ICC *icc = (PC_ICC*)pc->data;
118: PetscBool flg;
119: /* PetscReal dt[3];*/
121: PetscOptionsHead(PetscOptionsObject,"ICC Options");
122: PCSetFromOptions_Factor(PetscOptionsObject,pc);
124: PetscOptionsReal("-pc_factor_levels","levels of fill","PCFactorSetLevels",((PC_Factor*)icc)->info.levels,&((PC_Factor*)icc)->info.levels,&flg);
125: /*dt[0] = ((PC_Factor*)icc)->info.dt;
126: dt[1] = ((PC_Factor*)icc)->info.dtcol;
127: dt[2] = ((PC_Factor*)icc)->info.dtcount;
128: PetscInt dtmax = 3;
129: PetscOptionsRealArray("-pc_factor_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg);
130: if (flg) {
131: PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);
132: }
133: */
134: PetscOptionsTail();
135: return 0;
136: }
138: extern PetscErrorCode PCFactorSetDropTolerance_ILU(PC,PetscReal,PetscReal,PetscInt);
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: Notes:
153: Only implemented for some matrix formats. Not implemented in parallel.
155: For BAIJ matrices this implements a point block ICC.
157: The Manteuffel shift is only implemented for matrices with block size 1
159: By default, the Manteuffel is applied (for matrices with block size 1). Call PCFactorSetShiftType(pc,MAT_SHIFT_POSITIVE_DEFINITE);
160: to turn off the shift.
162: References:
163: . * - TONY F. CHAN AND HENK A. VAN DER VORST, Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS,
164: Chapter in Parallel Numerical Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
165: Science and Engineering, Kluwer.
167: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCSOR, MatOrderingType,
168: PCFactorSetZeroPivot(), PCFactorSetShiftType(), PCFactorSetShiftAmount(),
169: PCFactorSetFill(), PCFactorSetMatOrderingType(), PCFactorSetReuseOrdering(),
170: PCFactorSetLevels()
172: M*/
174: PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc)
175: {
176: PC_ICC *icc;
178: PetscNewLog(pc,&icc);
179: pc->data = (void*)icc;
180: PCFactorInitialize(pc, MAT_FACTOR_ICC);
182: ((PC_Factor*)icc)->info.fill = 1.0;
183: ((PC_Factor*)icc)->info.dtcol = PETSC_DEFAULT;
184: ((PC_Factor*)icc)->info.shifttype = (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE;
186: pc->ops->apply = PCApply_ICC;
187: pc->ops->matapply = PCMatApply_ICC;
188: pc->ops->applytranspose = PCApply_ICC;
189: pc->ops->setup = PCSetUp_ICC;
190: pc->ops->reset = PCReset_ICC;
191: pc->ops->destroy = PCDestroy_ICC;
192: pc->ops->setfromoptions = PCSetFromOptions_ICC;
193: pc->ops->view = PCView_Factor;
194: pc->ops->applysymmetricleft = PCApplySymmetricLeft_ICC;
195: pc->ops->applysymmetricright = PCApplySymmetricRight_ICC;
196: return 0;
197: }