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;
59: PetscCall(MatFactorGetSolverType(((PC_Factor *)icc)->fact, &solverpackage));
60: PetscCall(PCFactorSetMatSolverType(pc, solverpackage));
61: }
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
65: static PetscErrorCode PCReset_ICC(PC pc)
66: {
67: PC_ICC *icc = (PC_ICC *)pc->data;
69: PetscFunctionBegin;
70: PetscCall(MatDestroy(&((PC_Factor *)icc)->fact));
71: PetscFunctionReturn(PETSC_SUCCESS);
72: }
74: static PetscErrorCode PCDestroy_ICC(PC pc)
75: {
76: PC_ICC *icc = (PC_ICC *)pc->data;
78: PetscFunctionBegin;
79: PetscCall(PCReset_ICC(pc));
80: PetscCall(PetscFree(((PC_Factor *)icc)->ordering));
81: PetscCall(PetscFree(((PC_Factor *)icc)->solvertype));
82: PetscCall(PCFactorClearComposedFunctions(pc));
83: PetscCall(PetscFree(pc->data));
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: static PetscErrorCode PCApply_ICC(PC pc, Vec x, Vec y)
88: {
89: PC_ICC *icc = (PC_ICC *)pc->data;
91: PetscFunctionBegin;
92: PetscCall(MatSolve(((PC_Factor *)icc)->fact, x, y));
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: static PetscErrorCode PCMatApply_ICC(PC pc, Mat X, Mat Y)
97: {
98: PC_ICC *icc = (PC_ICC *)pc->data;
100: PetscFunctionBegin;
101: PetscCall(MatMatSolve(((PC_Factor *)icc)->fact, X, Y));
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc, Vec x, Vec y)
106: {
107: PC_ICC *icc = (PC_ICC *)pc->data;
109: PetscFunctionBegin;
110: PetscCall(MatForwardSolve(((PC_Factor *)icc)->fact, x, y));
111: PetscFunctionReturn(PETSC_SUCCESS);
112: }
114: static PetscErrorCode PCApplySymmetricRight_ICC(PC pc, Vec x, Vec y)
115: {
116: PC_ICC *icc = (PC_ICC *)pc->data;
118: PetscFunctionBegin;
119: PetscCall(MatBackwardSolve(((PC_Factor *)icc)->fact, x, y));
120: PetscFunctionReturn(PETSC_SUCCESS);
121: }
123: static PetscErrorCode PCSetFromOptions_ICC(PC pc, PetscOptionItems *PetscOptionsObject)
124: {
125: PC_ICC *icc = (PC_ICC *)pc->data;
126: PetscBool flg;
127: /* PetscReal dt[3];*/
129: PetscFunctionBegin;
130: PetscOptionsHeadBegin(PetscOptionsObject, "ICC Options");
131: PetscCall(PCSetFromOptions_Factor(pc, PetscOptionsObject));
133: PetscCall(PetscOptionsReal("-pc_factor_levels", "levels of fill", "PCFactorSetLevels", ((PC_Factor *)icc)->info.levels, &((PC_Factor *)icc)->info.levels, &flg));
134: /*dt[0] = ((PC_Factor*)icc)->info.dt;
135: dt[1] = ((PC_Factor*)icc)->info.dtcol;
136: dt[2] = ((PC_Factor*)icc)->info.dtcount;
137: PetscInt dtmax = 3;
138: PetscCall(PetscOptionsRealArray("-pc_factor_drop_tolerance","<dt,dtcol,maxrowcount>","PCFactorSetDropTolerance",dt,&dtmax,&flg));
139: if (flg) {
140: PetscCall(PCFactorSetDropTolerance(pc,dt[0],dt[1],(PetscInt)dt[2]));
141: }
142: */
143: PetscOptionsHeadEnd();
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: extern PetscErrorCode PCFactorSetDropTolerance_ILU(PC, PetscReal, PetscReal, PetscInt);
149: /*MC
150: PCICC - Incomplete Cholesky factorization preconditioners {cite}`chan1997approximate`
152: Options Database Keys:
153: + -pc_factor_levels <k> - number of levels of fill for ICC(k)
154: . -pc_factor_in_place - only for ICC(0) with natural ordering, reuses the space of the matrix for
155: its factorization (overwrites original matrix)
156: . -pc_factor_fill <nfill> - expected amount of fill in factored matrix compared to original matrix, nfill > 1
157: - -pc_factor_mat_ordering_type <natural,nd,1wd,rcm,qmd> - set the row/column ordering of the factored matrix
159: Level: beginner
161: Notes:
162: Only implemented for some matrix formats. Not implemented in parallel.
164: For `MATSEQBAIJ` matrices this implements a point block ICC.
166: By default, the Manteuffel shift {cite}`manteuffel1979shifted` is applied, for matrices with block size 1 only. Call `PCFactorSetShiftType`(pc,`MAT_SHIFT_POSITIVE_DEFINITE`);
167: to turn off the shift.
169: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSOR`, `MatOrderingType`, `PCILU`, `PCLU`, `PCCHOLESKY`,
170: `PCFactorSetZeroPivot()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`,
171: `PCFactorSetFill()`, `PCFactorSetMatOrderingType()`, `PCFactorSetReuseOrdering()`,
172: `PCFactorSetLevels()`
173: M*/
175: PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc)
176: {
177: PC_ICC *icc;
179: PetscFunctionBegin;
180: PetscCall(PetscNew(&icc));
181: pc->data = (void *)icc;
182: PetscCall(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: PetscFunctionReturn(PETSC_SUCCESS);
199: }