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;
12: const char *prefix;
14: pc->failedreason = PC_NOERROR;
16: PCGetOptionsPrefix(pc, &prefix);
17: MatSetOptionsPrefixFactor(pc->pmat, prefix);
19: MatSetErrorIfFailure(pc->pmat, pc->erroriffailure);
20: if (!pc->setupcalled) {
21: if (!((PC_Factor *)icc)->fact) MatGetFactor(pc->pmat, ((PC_Factor *)icc)->solvertype, MAT_FACTOR_ICC, &((PC_Factor *)icc)->fact);
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;
70: MatDestroy(&((PC_Factor *)icc)->fact);
71: return 0;
72: }
74: static PetscErrorCode PCDestroy_ICC(PC pc)
75: {
76: PC_ICC *icc = (PC_ICC *)pc->data;
78: PCReset_ICC(pc);
79: PetscFree(((PC_Factor *)icc)->ordering);
80: PetscFree(((PC_Factor *)icc)->solvertype);
81: PCFactorClearComposedFunctions(pc);
82: PetscFree(pc->data);
83: return 0;
84: }
86: static PetscErrorCode PCApply_ICC(PC pc, Vec x, Vec y)
87: {
88: PC_ICC *icc = (PC_ICC *)pc->data;
90: MatSolve(((PC_Factor *)icc)->fact, x, y);
91: return 0;
92: }
94: static PetscErrorCode PCMatApply_ICC(PC pc, Mat X, Mat Y)
95: {
96: PC_ICC *icc = (PC_ICC *)pc->data;
98: MatMatSolve(((PC_Factor *)icc)->fact, X, Y);
99: return 0;
100: }
102: static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc, Vec x, Vec y)
103: {
104: PC_ICC *icc = (PC_ICC *)pc->data;
106: MatForwardSolve(((PC_Factor *)icc)->fact, x, y);
107: return 0;
108: }
110: static PetscErrorCode PCApplySymmetricRight_ICC(PC pc, Vec x, Vec y)
111: {
112: PC_ICC *icc = (PC_ICC *)pc->data;
114: MatBackwardSolve(((PC_Factor *)icc)->fact, x, y);
115: return 0;
116: }
118: static PetscErrorCode PCSetFromOptions_ICC(PC pc, PetscOptionItems *PetscOptionsObject)
119: {
120: PC_ICC *icc = (PC_ICC *)pc->data;
121: PetscBool flg;
122: /* PetscReal dt[3];*/
124: PetscOptionsHeadBegin(PetscOptionsObject, "ICC Options");
125: PCSetFromOptions_Factor(pc, PetscOptionsObject);
127: PetscOptionsReal("-pc_factor_levels", "levels of fill", "PCFactorSetLevels", ((PC_Factor *)icc)->info.levels, &((PC_Factor *)icc)->info.levels, &flg);
128: /*dt[0] = ((PC_Factor*)icc)->info.dt;
129: dt[1] = ((PC_Factor*)icc)->info.dtcol;
130: dt[2] = ((PC_Factor*)icc)->info.dtcount;
131: PetscInt dtmax = 3;
132: PetscOptionsRealArray("-pc_factor_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg);
133: if (flg) {
134: PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]);
135: }
136: */
137: PetscOptionsHeadEnd();
138: return 0;
139: }
141: extern PetscErrorCode PCFactorSetDropTolerance_ILU(PC, PetscReal, PetscReal, PetscInt);
143: /*MC
144: PCICC - Incomplete Cholesky factorization preconditioners.
146: Options Database Keys:
147: + -pc_factor_levels <k> - number of levels of fill for ICC(k)
148: . -pc_factor_in_place - only for ICC(0) with natural ordering, reuses the space of the matrix for
149: its factorization (overwrites original matrix)
150: . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
151: - -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
153: Level: beginner
155: Notes:
156: Only implemented for some matrix formats. Not implemented in parallel.
158: For `MATSEQBAIJ` matrices this implements a point block ICC.
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: The Manteuffel shift is only implemented for matrices with block size 1
165: References:
166: . * - TONY F. CHAN AND HENK A. VAN DER VORST, Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS,
167: Chapter in Parallel Numerical Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
168: Science and Engineering, Kluwer.
170: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSOR`, `MatOrderingType`, `PCILU`, `PCLU`, `PCCHOLESKY`,
171: `PCFactorSetZeroPivot()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`,
172: `PCFactorSetFill()`, `PCFactorSetMatOrderingType()`, `PCFactorSetReuseOrdering()`,
173: `PCFactorSetLevels()`
174: M*/
176: PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc)
177: {
178: PC_ICC *icc;
180: PetscNew(&icc);
181: pc->data = (void *)icc;
182: PCFactorInitialize(pc, MAT_FACTOR_ICC);
184: ((PC_Factor *)icc)->info.fill = 1.0;
185: ((PC_Factor *)icc)->info.dtcol = PETSC_DEFAULT;
186: ((PC_Factor *)icc)->info.shifttype = (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE;
188: pc->ops->apply = PCApply_ICC;
189: pc->ops->matapply = PCMatApply_ICC;
190: pc->ops->applytranspose = PCApply_ICC;
191: pc->ops->setup = PCSetUp_ICC;
192: pc->ops->reset = PCReset_ICC;
193: pc->ops->destroy = PCDestroy_ICC;
194: pc->ops->setfromoptions = PCSetFromOptions_ICC;
195: pc->ops->view = PCView_Factor;
196: pc->ops->applysymmetricleft = PCApplySymmetricLeft_ICC;
197: pc->ops->applysymmetricright = PCApplySymmetricRight_ICC;
198: return 0;
199: }