1 #include <../src/ksp/pc/impls/factor/icc/icc.h> /*I "petscpc.h" I*/
2
PCSetUp_ICC(PC pc)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;
12
13 PetscFunctionBegin;
14 pc->failedreason = PC_NOERROR;
15
16 PetscCall(PCGetOptionsPrefix(pc, &prefix));
17 PetscCall(MatSetOptionsPrefixFactor(pc->pmat, prefix));
18
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;
40
41 PetscCall(ISDestroy(&cperm));
42 PetscCall(ISDestroy(&perm));
43
44 PetscCall(MatFactorGetError(((PC_Factor *)icc)->fact, &err));
45 if (err) { /* FactorSymbolic() fails */
46 pc->failedreason = (PCFailedReason)err;
47 PetscFunctionReturn(PETSC_SUCCESS);
48 }
49
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 }
55
56 PetscCall(PCFactorGetMatSolverType(pc, &stype));
57 if (!stype) {
58 MatSolverType solverpackage;
59
60 PetscCall(MatFactorGetSolverType(((PC_Factor *)icc)->fact, &solverpackage));
61 PetscCall(PCFactorSetMatSolverType(pc, solverpackage));
62 }
63 PetscFunctionReturn(PETSC_SUCCESS);
64 }
65
PCReset_ICC(PC pc)66 static PetscErrorCode PCReset_ICC(PC pc)
67 {
68 PC_ICC *icc = (PC_ICC *)pc->data;
69
70 PetscFunctionBegin;
71 PetscCall(MatDestroy(&((PC_Factor *)icc)->fact));
72 PetscFunctionReturn(PETSC_SUCCESS);
73 }
74
PCDestroy_ICC(PC pc)75 static PetscErrorCode PCDestroy_ICC(PC pc)
76 {
77 PC_ICC *icc = (PC_ICC *)pc->data;
78
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 }
87
PCApply_ICC(PC pc,Vec x,Vec y)88 static PetscErrorCode PCApply_ICC(PC pc, Vec x, Vec y)
89 {
90 PC_ICC *icc = (PC_ICC *)pc->data;
91
92 PetscFunctionBegin;
93 PetscCall(MatSolve(((PC_Factor *)icc)->fact, x, y));
94 PetscFunctionReturn(PETSC_SUCCESS);
95 }
96
PCMatApply_ICC(PC pc,Mat X,Mat Y)97 static PetscErrorCode PCMatApply_ICC(PC pc, Mat X, Mat Y)
98 {
99 PC_ICC *icc = (PC_ICC *)pc->data;
100
101 PetscFunctionBegin;
102 PetscCall(MatMatSolve(((PC_Factor *)icc)->fact, X, Y));
103 PetscFunctionReturn(PETSC_SUCCESS);
104 }
105
PCApplySymmetricLeft_ICC(PC pc,Vec x,Vec y)106 static PetscErrorCode PCApplySymmetricLeft_ICC(PC pc, Vec x, Vec y)
107 {
108 PC_ICC *icc = (PC_ICC *)pc->data;
109
110 PetscFunctionBegin;
111 PetscCall(MatForwardSolve(((PC_Factor *)icc)->fact, x, y));
112 PetscFunctionReturn(PETSC_SUCCESS);
113 }
114
PCApplySymmetricRight_ICC(PC pc,Vec x,Vec y)115 static PetscErrorCode PCApplySymmetricRight_ICC(PC pc, Vec x, Vec y)
116 {
117 PC_ICC *icc = (PC_ICC *)pc->data;
118
119 PetscFunctionBegin;
120 PetscCall(MatBackwardSolve(((PC_Factor *)icc)->fact, x, y));
121 PetscFunctionReturn(PETSC_SUCCESS);
122 }
123
PCSetFromOptions_ICC(PC pc,PetscOptionItems PetscOptionsObject)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];*/
129
130 PetscFunctionBegin;
131 PetscOptionsHeadBegin(PetscOptionsObject, "ICC Options");
132 PetscCall(PCSetFromOptions_Factor(pc, PetscOptionsObject));
133
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 }
145
146 /*MC
147 PCICC - Incomplete Cholesky factorization preconditioners {cite}`chan1997approximate`
148
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
155
156 Level: beginner
157
158 Notes:
159 Only implemented for some matrix formats. Not implemented in parallel.
160
161 For `MATSEQBAIJ` matrices this implements a point block ICC.
162
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.
165
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*/
171
PCCreate_ICC(PC pc)172 PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc)
173 {
174 PC_ICC *icc;
175
176 PetscFunctionBegin;
177 PetscCall(PetscNew(&icc));
178 pc->data = (void *)icc;
179 PetscCall(PCFactorInitialize(pc, MAT_FACTOR_ICC));
180
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;
184
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 }
198