Actual source code: icc.c
1: #include <../src/ksp/pc/impls/factor/icc/icc.h>
3: static PetscErrorCode PCSetUp_ICC(PC pc)
4: {
5: PC_ICC *icc = (PC_ICC *)pc->data;
6: IS perm = NULL, cperm = NULL;
7: MatInfo info;
8: MatSolverType stype;
9: MatFactorError err;
10: PetscBool canuseordering;
11: const char *prefix;
13: PetscFunctionBegin;
14: pc->failedreason = PC_NOERROR;
16: PetscCall(PCGetOptionsPrefix(pc, &prefix));
17: PetscCall(MatSetOptionsPrefixFactor(pc->pmat, prefix));
19: PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
20: if (!pc->setupcalled) {
21: PetscCall(PCFactorSetUpMatSolverType(pc));
22: PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)icc)->fact, &canuseordering));
23: if (canuseordering) {
24: PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
25: PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)icc)->ordering, &perm, &cperm));
26: }
27: PetscCall(MatICCFactorSymbolic(((PC_Factor *)icc)->fact, pc->pmat, perm, &((PC_Factor *)icc)->info));
28: } else if (pc->flag != SAME_NONZERO_PATTERN) {
29: PetscCall(MatDestroy(&((PC_Factor *)icc)->fact));
30: PetscCall(PCFactorSetUpMatSolverType(pc));
31: PetscCall(MatFactorGetCanUseOrdering(((PC_Factor *)icc)->fact, &canuseordering));
32: if (canuseordering) {
33: PetscCall(PCFactorSetDefaultOrdering_Factor(pc));
34: PetscCall(MatGetOrdering(pc->pmat, ((PC_Factor *)icc)->ordering, &perm, &cperm));
35: }
36: PetscCall(MatICCFactorSymbolic(((PC_Factor *)icc)->fact, pc->pmat, perm, &((PC_Factor *)icc)->info));
37: }
38: PetscCall(MatGetInfo(((PC_Factor *)icc)->fact, MAT_LOCAL, &info));
39: icc->hdr.actualfill = info.fill_ratio_needed;
41: PetscCall(ISDestroy(&cperm));
42: PetscCall(ISDestroy(&perm));
44: PetscCall(MatFactorGetError(((PC_Factor *)icc)->fact, &err));
45: if (err) { /* FactorSymbolic() fails */
46: pc->failedreason = (PCFailedReason)err;
47: PetscFunctionReturn(PETSC_SUCCESS);
48: }
50: PetscCall(MatCholeskyFactorNumeric(((PC_Factor *)icc)->fact, pc->pmat, &((PC_Factor *)icc)->info));
51: PetscCall(MatFactorGetError(((PC_Factor *)icc)->fact, &err));
52: if (err) { /* FactorNumeric() fails */
53: pc->failedreason = (PCFailedReason)err;
54: }
56: PetscCall(PCFactorGetMatSolverType(pc, &stype));
57: if (!stype) {
58: MatSolverType solverpackage;
60: PetscCall(MatFactorGetSolverType(((PC_Factor *)icc)->fact, &solverpackage));
61: PetscCall(PCFactorSetMatSolverType(pc, solverpackage));
62: }
63: PetscFunctionReturn(PETSC_SUCCESS);
64: }
66: static PetscErrorCode PCReset_ICC(PC pc)
67: {
68: PC_ICC *icc = (PC_ICC *)pc->data;
70: PetscFunctionBegin;
71: PetscCall(MatDestroy(&((PC_Factor *)icc)->fact));
72: PetscFunctionReturn(PETSC_SUCCESS);
73: }
75: static PetscErrorCode PCDestroy_ICC(PC pc)
76: {
77: PC_ICC *icc = (PC_ICC *)pc->data;
79: PetscFunctionBegin;
80: PetscCall(PCReset_ICC(pc));
81: PetscCall(PetscFree(((PC_Factor *)icc)->ordering));
82: PetscCall(PetscFree(((PC_Factor *)icc)->solvertype));
83: PetscCall(PCFactorClearComposedFunctions(pc));
84: PetscCall(PetscFree(pc->data));
85: PetscFunctionReturn(PETSC_SUCCESS);
86: }
88: static PetscErrorCode PCApply_ICC(PC pc, Vec x, Vec y)
89: {
90: PC_ICC *icc = (PC_ICC *)pc->data;
92: PetscFunctionBegin;
93: PetscCall(MatSolve(((PC_Factor *)icc)->fact, x, y));
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
97: static PetscErrorCode PCMatApply_ICC(PC pc, Mat X, Mat Y)
98: {
99: PC_ICC *icc = (PC_ICC *)pc->data;
101: PetscFunctionBegin;
102: PetscCall(MatMatSolve(((PC_Factor *)icc)->fact, X, Y));
103: PetscFunctionReturn(PETSC_SUCCESS);
104: }
106: static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc, Vec x, Vec y)
107: {
108: PC_ICC *icc = (PC_ICC *)pc->data;
110: PetscFunctionBegin;
111: PetscCall(MatForwardSolve(((PC_Factor *)icc)->fact, x, y));
112: PetscFunctionReturn(PETSC_SUCCESS);
113: }
115: static PetscErrorCode PCApplySymmetricRight_ICC(PC pc, Vec x, Vec y)
116: {
117: PC_ICC *icc = (PC_ICC *)pc->data;
119: PetscFunctionBegin;
120: PetscCall(MatBackwardSolve(((PC_Factor *)icc)->fact, x, y));
121: PetscFunctionReturn(PETSC_SUCCESS);
122: }
124: static PetscErrorCode PCSetFromOptions_ICC(PC pc, PetscOptionItems PetscOptionsObject)
125: {
126: PC_ICC *icc = (PC_ICC *)pc->data;
127: PetscBool flg;
128: /* PetscReal dt[3];*/
130: PetscFunctionBegin;
131: PetscOptionsHeadBegin(PetscOptionsObject, "ICC Options");
132: PetscCall(PCSetFromOptions_Factor(pc, PetscOptionsObject));
134: PetscCall(PetscOptionsReal("-pc_factor_levels", "levels of fill", "PCFactorSetLevels", ((PC_Factor *)icc)->info.levels, &((PC_Factor *)icc)->info.levels, &flg));
135: /*dt[0] = ((PC_Factor*)icc)->info.dt;
136: dt[1] = ((PC_Factor*)icc)->info.dtcol;
137: dt[2] = ((PC_Factor*)icc)->info.dtcount;
138: PetscInt dtmax = 3;
139: PetscCall(PetscOptionsRealArray("-pc_factor_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg));
140: if (flg) PetscCall(PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]));
141: */
142: PetscOptionsHeadEnd();
143: PetscFunctionReturn(PETSC_SUCCESS);
144: }
146: /*MC
147: PCICC - Incomplete Cholesky factorization preconditioners {cite}`chan1997approximate`
149: Options Database Keys:
150: + -pc_factor_levels <k> - number of levels of fill for ICC(k)
151: . -pc_factor_in_place - only for ICC(0) with natural ordering, reuses the space of the matrix for
152: its factorization (overwrites original matrix)
153: . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
154: - -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
156: Level: beginner
158: Notes:
159: Only implemented for some matrix formats. Not implemented in parallel.
161: For `MATSEQBAIJ` matrices this implements a point block ICC.
163: By default, the Manteuffel shift {cite}`manteuffel1979shifted` is applied, for matrices with block size 1 only. Call `PCFactorSetShiftType`(pc,`MAT_SHIFT_POSITIVE_DEFINITE`);
164: to turn off the shift.
166: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSOR`, `MatOrderingType`, `PCILU`, `PCLU`, `PCCHOLESKY`,
167: `PCFactorSetZeroPivot()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`,
168: `PCFactorSetFill()`, `PCFactorSetMatOrderingType()`, `PCFactorSetReuseOrdering()`,
169: `PCFactorSetLevels()`
170: M*/
172: PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc)
173: {
174: PC_ICC *icc;
176: PetscFunctionBegin;
177: PetscCall(PetscNew(&icc));
178: pc->data = (void *)icc;
179: PetscCall(PCFactorInitialize(pc, MAT_FACTOR_ICC));
181: ((PC_Factor *)icc)->info.fill = 1.0;
182: ((PC_Factor *)icc)->info.dtcol = PETSC_DEFAULT;
183: ((PC_Factor *)icc)->info.shifttype = (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE;
185: pc->ops->apply = PCApply_ICC;
186: pc->ops->matapply = PCMatApply_ICC;
187: pc->ops->applytranspose = PCApply_ICC;
188: pc->ops->matapplytranspose = PCMatApply_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: PetscFunctionReturn(PETSC_SUCCESS);
197: }