1 #include <petsc/private/pcimpl.h>
2 #include <petsc/private/kspimpl.h>
3 #include <petscksp.h> /*I "petscksp.h" I*/
4
5 typedef struct {
6 KSP ksp;
7 PetscInt its; /* total number of iterations KSP uses */
8 } PC_KSP;
9
PCKSPCreateKSP_KSP(PC pc)10 static PetscErrorCode PCKSPCreateKSP_KSP(PC pc)
11 {
12 PetscErrorCode ierr;
13 const char *prefix;
14 PC_KSP *jac = (PC_KSP*)pc->data;
15 DM dm;
16
17 PetscFunctionBegin;
18 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&jac->ksp);CHKERRQ(ierr);
19 ierr = KSPSetErrorIfNotConverged(jac->ksp,pc->erroriffailure);CHKERRQ(ierr);
20 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->ksp,(PetscObject)pc,1);CHKERRQ(ierr);
21 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr);
22 ierr = KSPSetOptionsPrefix(jac->ksp,prefix);CHKERRQ(ierr);
23 ierr = KSPAppendOptionsPrefix(jac->ksp,"ksp_");CHKERRQ(ierr);
24 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr);
25 if (dm) {
26 ierr = KSPSetDM(jac->ksp, dm);CHKERRQ(ierr);
27 ierr = KSPSetDMActive(jac->ksp, PETSC_FALSE);CHKERRQ(ierr);
28 }
29 PetscFunctionReturn(0);
30 }
31
PCApply_KSP(PC pc,Vec x,Vec y)32 static PetscErrorCode PCApply_KSP(PC pc,Vec x,Vec y)
33 {
34 PetscErrorCode ierr;
35 PetscInt its;
36 PC_KSP *jac = (PC_KSP*)pc->data;
37
38 PetscFunctionBegin;
39 if (jac->ksp->presolve) {
40 ierr = VecCopy(x,y);CHKERRQ(ierr);
41 ierr = KSPSolve(jac->ksp,y,y);CHKERRQ(ierr);
42 } else {
43 ierr = KSPSolve(jac->ksp,x,y);CHKERRQ(ierr);
44 }
45 ierr = KSPCheckSolve(jac->ksp,pc,y);CHKERRQ(ierr);
46 ierr = KSPGetIterationNumber(jac->ksp,&its);CHKERRQ(ierr);
47 jac->its += its;
48 PetscFunctionReturn(0);
49 }
50
PCMatApply_KSP(PC pc,Mat X,Mat Y)51 static PetscErrorCode PCMatApply_KSP(PC pc,Mat X,Mat Y)
52 {
53 PetscErrorCode ierr;
54 PetscInt its;
55 PC_KSP *jac = (PC_KSP*)pc->data;
56
57 PetscFunctionBegin;
58 if (jac->ksp->presolve) {
59 ierr = MatCopy(X,Y,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
60 ierr = KSPMatSolve(jac->ksp,Y,Y);CHKERRQ(ierr); /* TODO FIXME: this will fail since KSPMatSolve does not allow inplace solve yet */
61 } else {
62 ierr = KSPMatSolve(jac->ksp,X,Y);CHKERRQ(ierr);
63 }
64 ierr = KSPCheckSolve(jac->ksp,pc,NULL);CHKERRQ(ierr);
65 ierr = KSPGetIterationNumber(jac->ksp,&its);CHKERRQ(ierr);
66 jac->its += its;
67 PetscFunctionReturn(0);
68 }
69
PCApplyTranspose_KSP(PC pc,Vec x,Vec y)70 static PetscErrorCode PCApplyTranspose_KSP(PC pc,Vec x,Vec y)
71 {
72 PetscErrorCode ierr;
73 PetscInt its;
74 PC_KSP *jac = (PC_KSP*)pc->data;
75
76 PetscFunctionBegin;
77 if (jac->ksp->presolve) {
78 ierr = VecCopy(x,y);CHKERRQ(ierr);
79 ierr = KSPSolve(jac->ksp,y,y);CHKERRQ(ierr);
80 } else {
81 ierr = KSPSolveTranspose(jac->ksp,x,y);CHKERRQ(ierr);
82 }
83 ierr = KSPCheckSolve(jac->ksp,pc,y);CHKERRQ(ierr);
84 ierr = KSPGetIterationNumber(jac->ksp,&its);CHKERRQ(ierr);
85 jac->its += its;
86 PetscFunctionReturn(0);
87 }
88
PCSetUp_KSP(PC pc)89 static PetscErrorCode PCSetUp_KSP(PC pc)
90 {
91 PetscErrorCode ierr;
92 PC_KSP *jac = (PC_KSP*)pc->data;
93 Mat mat;
94
95 PetscFunctionBegin;
96 if (!jac->ksp) {
97 ierr = PCKSPCreateKSP_KSP(pc);CHKERRQ(ierr);
98 ierr = KSPSetFromOptions(jac->ksp);CHKERRQ(ierr);
99 }
100 if (pc->useAmat) mat = pc->mat;
101 else mat = pc->pmat;
102 ierr = KSPSetOperators(jac->ksp,mat,pc->pmat);CHKERRQ(ierr);
103 ierr = KSPSetUp(jac->ksp);CHKERRQ(ierr);
104 PetscFunctionReturn(0);
105 }
106
107 /* Default destroy, if it has never been setup */
PCReset_KSP(PC pc)108 static PetscErrorCode PCReset_KSP(PC pc)
109 {
110 PC_KSP *jac = (PC_KSP*)pc->data;
111 PetscErrorCode ierr;
112
113 PetscFunctionBegin;
114 ierr = KSPDestroy(&jac->ksp);CHKERRQ(ierr);
115 PetscFunctionReturn(0);
116 }
117
PCDestroy_KSP(PC pc)118 static PetscErrorCode PCDestroy_KSP(PC pc)
119 {
120 PC_KSP *jac = (PC_KSP*)pc->data;
121 PetscErrorCode ierr;
122
123 PetscFunctionBegin;
124 ierr = KSPDestroy(&jac->ksp);CHKERRQ(ierr);
125 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",NULL);CHKERRQ(ierr);
126 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCKSPSetKSP_C",NULL);CHKERRQ(ierr);
127 ierr = PetscFree(pc->data);CHKERRQ(ierr);
128 PetscFunctionReturn(0);
129 }
130
PCView_KSP(PC pc,PetscViewer viewer)131 static PetscErrorCode PCView_KSP(PC pc,PetscViewer viewer)
132 {
133 PC_KSP *jac = (PC_KSP*)pc->data;
134 PetscErrorCode ierr;
135 PetscBool iascii;
136
137 PetscFunctionBegin;
138 if (!jac->ksp) {ierr = PCKSPCreateKSP_KSP(pc);CHKERRQ(ierr);}
139 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
140 if (iascii) {
141 if (pc->useAmat) {
142 ierr = PetscViewerASCIIPrintf(viewer," Using Amat (not Pmat) as operator on inner solve\n");CHKERRQ(ierr);
143 }
144 ierr = PetscViewerASCIIPrintf(viewer," KSP and PC on KSP preconditioner follow\n");CHKERRQ(ierr);
145 ierr = PetscViewerASCIIPrintf(viewer," ---------------------------------\n");CHKERRQ(ierr);
146 }
147 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
148 ierr = KSPView(jac->ksp,viewer);CHKERRQ(ierr);
149 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
150 if (iascii) {
151 ierr = PetscViewerASCIIPrintf(viewer," ---------------------------------\n");CHKERRQ(ierr);
152 }
153 PetscFunctionReturn(0);
154 }
155
PCKSPSetKSP_KSP(PC pc,KSP ksp)156 static PetscErrorCode PCKSPSetKSP_KSP(PC pc,KSP ksp)
157 {
158 PC_KSP *jac = (PC_KSP*)pc->data;
159 PetscErrorCode ierr;
160
161 PetscFunctionBegin;
162 ierr = PetscObjectReference((PetscObject)ksp);CHKERRQ(ierr);
163 ierr = KSPDestroy(&jac->ksp);CHKERRQ(ierr);
164 jac->ksp = ksp;
165 PetscFunctionReturn(0);
166 }
167
168 /*@
169 PCKSPSetKSP - Sets the KSP context for a KSP PC.
170
171 Collective on PC
172
173 Input Parameter:
174 + pc - the preconditioner context
175 - ksp - the KSP solver
176
177 Notes:
178 The PC and the KSP must have the same communicator
179
180 Level: advanced
181
182 @*/
PCKSPSetKSP(PC pc,KSP ksp)183 PetscErrorCode PCKSPSetKSP(PC pc,KSP ksp)
184 {
185 PetscErrorCode ierr;
186
187 PetscFunctionBegin;
188 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
189 PetscValidHeaderSpecific(ksp,KSP_CLASSID,2);
190 PetscCheckSameComm(pc,1,ksp,2);
191 ierr = PetscTryMethod(pc,"PCKSPSetKSP_C",(PC,KSP),(pc,ksp));CHKERRQ(ierr);
192 PetscFunctionReturn(0);
193 }
194
PCKSPGetKSP_KSP(PC pc,KSP * ksp)195 static PetscErrorCode PCKSPGetKSP_KSP(PC pc,KSP *ksp)
196 {
197 PC_KSP *jac = (PC_KSP*)pc->data;
198 PetscErrorCode ierr;
199
200 PetscFunctionBegin;
201 if (!jac->ksp) {ierr = PCKSPCreateKSP_KSP(pc);CHKERRQ(ierr);}
202 *ksp = jac->ksp;
203 PetscFunctionReturn(0);
204 }
205
206 /*@
207 PCKSPGetKSP - Gets the KSP context for a KSP PC.
208
209 Not Collective but KSP returned is parallel if PC was parallel
210
211 Input Parameter:
212 . pc - the preconditioner context
213
214 Output Parameters:
215 . ksp - the KSP solver
216
217 Notes:
218 You must call KSPSetUp() before calling PCKSPGetKSP().
219
220 If the PC is not a PCKSP object it raises an error
221
222 Level: advanced
223
224 @*/
PCKSPGetKSP(PC pc,KSP * ksp)225 PetscErrorCode PCKSPGetKSP(PC pc,KSP *ksp)
226 {
227 PetscErrorCode ierr;
228
229 PetscFunctionBegin;
230 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
231 PetscValidPointer(ksp,2);
232 ierr = PetscUseMethod(pc,"PCKSPGetKSP_C",(PC,KSP*),(pc,ksp));CHKERRQ(ierr);
233 PetscFunctionReturn(0);
234 }
235
PCSetFromOptions_KSP(PetscOptionItems * PetscOptionsObject,PC pc)236 static PetscErrorCode PCSetFromOptions_KSP(PetscOptionItems *PetscOptionsObject,PC pc)
237 {
238 PC_KSP *jac = (PC_KSP*)pc->data;
239 PetscErrorCode ierr;
240
241 PetscFunctionBegin;
242 ierr = PetscOptionsHead(PetscOptionsObject,"PC KSP options");CHKERRQ(ierr);
243 if (jac->ksp) {
244 ierr = KSPSetFromOptions(jac->ksp);CHKERRQ(ierr);
245 }
246 ierr = PetscOptionsTail();CHKERRQ(ierr);
247 PetscFunctionReturn(0);
248 }
249
250 /* ----------------------------------------------------------------------------------*/
251
252 /*MC
253 PCKSP - Defines a preconditioner that can consist of any KSP solver.
254 This allows, for example, embedding a Krylov method inside a preconditioner.
255
256 Options Database Key:
257 . -pc_use_amat - use the matrix that defines the linear system, Amat as the matrix for the
258 inner solver, otherwise by default it uses the matrix used to construct
259 the preconditioner, Pmat (see PCSetOperators())
260
261 Level: intermediate
262
263 Notes:
264 Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
265 the incorrect answer) unless you use KSPFGMRES as the other Krylov method
266
267 Developer Notes:
268 If the outer Krylov method has a nonzero initial guess it will compute a new residual based on that initial guess
269 and pass that as the right hand side into this KSP (and hence this KSP will always have a zero initial guess). For all outer Krylov methods
270 except Richardson this is neccessary since Krylov methods, even the flexible ones, need to "see" the result of the action of the preconditioner on the
271 input (current residual) vector, the action of the preconditioner cannot depend also on some other vector (the "initial guess"). For
272 KSPRICHARDSON it is possible to provide a PCApplyRichardson_PCKSP() that short circuits returning to the KSP object at each iteration to compute the
273 residual, see for example PCApplyRichardson_SOR(). We do not implement a PCApplyRichardson_PCKSP() because (1) using a KSP directly inside a Richardson
274 is not an efficient algorithm anyways and (2) implementing it for its > 1 would essentially require that we implement Richardson (reimplementing the
275 Richardson code) inside the PCApplyRichardson_PCKSP() leading to duplicate code.
276
277 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
278 PCSHELL, PCCOMPOSITE, PCSetUseAmat(), PCKSPGetKSP()
279
280 M*/
281
PCCreate_KSP(PC pc)282 PETSC_EXTERN PetscErrorCode PCCreate_KSP(PC pc)
283 {
284 PetscErrorCode ierr;
285 PC_KSP *jac;
286
287 PetscFunctionBegin;
288 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr);
289 pc->data = (void*)jac;
290
291 ierr = PetscMemzero(pc->ops,sizeof(struct _PCOps));CHKERRQ(ierr);
292 pc->ops->apply = PCApply_KSP;
293 pc->ops->matapply = PCMatApply_KSP;
294 pc->ops->applytranspose = PCApplyTranspose_KSP;
295 pc->ops->setup = PCSetUp_KSP;
296 pc->ops->reset = PCReset_KSP;
297 pc->ops->destroy = PCDestroy_KSP;
298 pc->ops->setfromoptions = PCSetFromOptions_KSP;
299 pc->ops->view = PCView_KSP;
300
301 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",PCKSPGetKSP_KSP);CHKERRQ(ierr);
302 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCKSPSetKSP_C",PCKSPSetKSP_KSP);CHKERRQ(ierr);
303 PetscFunctionReturn(0);
304 }
305