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: PetscFunctionBegin;
15: pc->failedreason = PC_NOERROR;
17: PetscCall(PCGetOptionsPrefix(pc, &prefix));
18: PetscCall(MatSetOptionsPrefixFactor(pc->pmat, prefix));
20: PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
21: if (!pc->setupcalled) {
22: if (!((PC_Factor *)icc)->fact) PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)icc)->solvertype, MAT_FACTOR_ICC, &((PC_Factor *)icc)->fact));
23: PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)icc)->fact, &canuseordering));
24: if (canuseordering) {
25: PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
26: PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)icc)->ordering, &perm, &cperm));
27: }
28: PetscCall(MatICCFactorSymbolic(((PC_Factor *)icc)->fact, pc->pmat, perm, &((PC_Factor *)icc)->info));
29: } else if (pc->flag != SAME_NONZERO_PATTERN) {
30: PetscBool canuseordering;
31: PetscCall(MatDestroy(&((PC_Factor *)icc)->fact));
32: PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)icc)->solvertype, MAT_FACTOR_ICC, &((PC_Factor *)icc)->fact));
33: PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)icc)->fact, &canuseordering));
34: if (canuseordering) {
35: PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
36: PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)icc)->ordering, &perm, &cperm));
37: }
38: PetscCall(MatICCFactorSymbolic(((PC_Factor *)icc)->fact, pc->pmat, perm, &((PC_Factor *)icc)->info));
39: }
40: PetscCall(MatGetInfo(((PC_Factor *)icc)->fact, MAT_LOCAL, &info));
41: icc->hdr.actualfill = info.fill_ratio_needed;
43: PetscCall(ISDestroy(&cperm));
44: PetscCall(ISDestroy(&perm));
46: PetscCall(MatFactorGetError(((PC_Factor *)icc)->fact, &err));
47: if (err) { /* FactorSymbolic() fails */
48: pc->failedreason = (PCFailedReason)err;
49: PetscFunctionReturn(PETSC_SUCCESS);
50: }
52: PetscCall(MatCholeskyFactorNumeric(((PC_Factor *)icc)->fact, pc->pmat, &((PC_Factor *)icc)->info));
53: PetscCall(MatFactorGetError(((PC_Factor *)icc)->fact, &err));
54: if (err) { /* FactorNumeric() fails */
55: pc->failedreason = (PCFailedReason)err;
56: }
58: PetscCall(PCFactorGetMatSolverType(pc, &stype));
59: if (!stype) {
60: MatSolverType solverpackage;
61: PetscCall(MatFactorGetSolverType(((PC_Factor *)icc)->fact, &solverpackage));
62: PetscCall(PCFactorSetMatSolverType(pc, solverpackage));
63: }
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: static PetscErrorCode PCReset_ICC(PC pc)
68: {
69: PC_ICC *icc = (PC_ICC *)pc->data;
71: PetscFunctionBegin;
72: PetscCall(MatDestroy(&((PC_Factor *)icc)->fact));
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: static PetscErrorCode PCDestroy_ICC(PC pc)
77: {
78: PC_ICC *icc = (PC_ICC *)pc->data;
80: PetscFunctionBegin;
81: PetscCall(PCReset_ICC(pc));
82: PetscCall(PetscFree(((PC_Factor *)icc)->ordering));
83: PetscCall(PetscFree(((PC_Factor *)icc)->solvertype));
84: PetscCall(PCFactorClearComposedFunctions(pc));
85: PetscCall(PetscFree(pc->data));
86: PetscFunctionReturn(PETSC_SUCCESS);
87: }
89: static PetscErrorCode PCApply_ICC(PC pc, Vec x, Vec y)
90: {
91: PC_ICC *icc = (PC_ICC *)pc->data;
93: PetscFunctionBegin;
94: PetscCall(MatSolve(((PC_Factor *)icc)->fact, x, y));
95: PetscFunctionReturn(PETSC_SUCCESS);
96: }
98: static PetscErrorCode PCMatApply_ICC(PC pc, Mat X, Mat Y)
99: {
100: PC_ICC *icc = (PC_ICC *)pc->data;
102: PetscFunctionBegin;
103: PetscCall(MatMatSolve(((PC_Factor *)icc)->fact, X, Y));
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
107: static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc, Vec x, Vec y)
108: {
109: PC_ICC *icc = (PC_ICC *)pc->data;
111: PetscFunctionBegin;
112: PetscCall(MatForwardSolve(((PC_Factor *)icc)->fact, x, y));
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: static PetscErrorCode PCApplySymmetricRight_ICC(PC pc, Vec x, Vec y)
117: {
118: PC_ICC *icc = (PC_ICC *)pc->data;
120: PetscFunctionBegin;
121: PetscCall(MatBackwardSolve(((PC_Factor *)icc)->fact, x, y));
122: PetscFunctionReturn(PETSC_SUCCESS);
123: }
125: static PetscErrorCode PCSetFromOptions_ICC(PC pc, PetscOptionItems *PetscOptionsObject)
126: {
127: PC_ICC *icc = (PC_ICC *)pc->data;
128: PetscBool flg;
129: /* PetscReal dt[3];*/
131: PetscFunctionBegin;
132: PetscOptionsHeadBegin(PetscOptionsObject, "ICC Options");
133: PetscCall(PCSetFromOptions_Factor(pc, PetscOptionsObject));
135: PetscCall(PetscOptionsReal("-pc_factor_levels", "levels of fill", "PCFactorSetLevels", ((PC_Factor *)icc)->info.levels, &((PC_Factor *)icc)->info.levels, &flg));
136: /*dt[0] = ((PC_Factor*)icc)->info.dt;
137: dt[1] = ((PC_Factor*)icc)->info.dtcol;
138: dt[2] = ((PC_Factor*)icc)->info.dtcount;
139: PetscInt dtmax = 3;
140: PetscCall(PetscOptionsRealArray("-pc_factor_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg));
141: if (flg) {
142: PetscCall(PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]));
143: }
144: */
145: PetscOptionsHeadEnd();
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: extern PetscErrorCode PCFactorSetDropTolerance_ILU(PC, PetscReal, PetscReal, PetscInt);
151: /*MC
152: PCICC - Incomplete Cholesky factorization preconditioners.
154: Options Database Keys:
155: + -pc_factor_levels <k> - number of levels of fill for ICC(k)
156: . -pc_factor_in_place - only for ICC(0) with natural ordering, reuses the space of the matrix for
157: its factorization (overwrites original matrix)
158: . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
159: - -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
161: Level: beginner
163: Notes:
164: Only implemented for some matrix formats. Not implemented in parallel.
166: For `MATSEQBAIJ` matrices this implements a point block ICC.
168: By default, the Manteuffel is applied (for matrices with block size 1). Call `PCFactorSetShiftType`(pc,`MAT_SHIFT_POSITIVE_DEFINITE`);
169: to turn off the shift.
171: The Manteuffel shift is only implemented for matrices with block size 1
173: References:
174: . * - TONY F. CHAN AND HENK A. VAN DER VORST, Review article: APPROXIMATE AND INCOMPLETE FACTORIZATIONS,
175: Chapter in Parallel Numerical Algorithms, edited by D. Keyes, A. Semah, V. Venkatakrishnan, ICASE/LaRC Interdisciplinary Series in
176: Science and Engineering, Kluwer.
178: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSOR`, `MatOrderingType`, `PCILU`, `PCLU`, `PCCHOLESKY`,
179: `PCFactorSetZeroPivot()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`,
180: `PCFactorSetFill()`, `PCFactorSetMatOrderingType()`, `PCFactorSetReuseOrdering()`,
181: `PCFactorSetLevels()`
182: M*/
184: PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc)
185: {
186: PC_ICC *icc;
188: PetscFunctionBegin;
189: PetscCall(PetscNew(&icc));
190: pc->data = (void *)icc;
191: PetscCall(PCFactorInitialize(pc, MAT_FACTOR_ICC));
193: ((PC_Factor *)icc)->info.fill = 1.0;
194: ((PC_Factor *)icc)->info.dtcol = PETSC_DEFAULT;
195: ((PC_Factor *)icc)->info.shifttype = (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE;
197: pc->ops->apply = PCApply_ICC;
198: pc->ops->matapply = PCMatApply_ICC;
199: pc->ops->applytranspose = PCApply_ICC;
200: pc->ops->setup = PCSetUp_ICC;
201: pc->ops->reset = PCReset_ICC;
202: pc->ops->destroy = PCDestroy_ICC;
203: pc->ops->setfromoptions = PCSetFromOptions_ICC;
204: pc->ops->view = PCView_Factor;
205: pc->ops->applysymmetricleft = PCApplySymmetricLeft_ICC;
206: pc->ops->applysymmetricright = PCApplySymmetricRight_ICC;
207: PetscFunctionReturn(PETSC_SUCCESS);
208: }