Actual source code: pcksp.c
petsc-3.12.5 2020-03-29
1: #include <petsc/private/pcimpl.h>
2: #include <petsc/private/kspimpl.h>
3: #include <petscksp.h>
5: typedef struct {
6: KSP ksp;
7: PetscInt its; /* total number of iterations KSP uses */
8: } PC_KSP;
10: static PetscErrorCode PCKSPCreateKSP_KSP(PC pc)
11: {
13: const char *prefix;
14: PC_KSP *jac = (PC_KSP*)pc->data;
15: DM dm;
18: KSPCreate(PetscObjectComm((PetscObject)pc),&jac->ksp);
19: KSPSetErrorIfNotConverged(jac->ksp,pc->erroriffailure);
20: PetscObjectIncrementTabLevel((PetscObject)jac->ksp,(PetscObject)pc,1);
21: PCGetOptionsPrefix(pc,&prefix);
22: KSPSetOptionsPrefix(jac->ksp,prefix);
23: KSPAppendOptionsPrefix(jac->ksp,"ksp_");
24: PCGetDM(pc, &dm);
25: if (dm) {
26: KSPSetDM(jac->ksp, dm);
27: KSPSetDMActive(jac->ksp, PETSC_FALSE);
28: }
29: return(0);
30: }
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;
39: if (jac->ksp->presolve) {
40: VecCopy(x,y);
41: KSPSolve(jac->ksp,y,y);
42: } else {
43: KSPSolve(jac->ksp,x,y);
44: }
45: KSPCheckSolve(jac->ksp,pc,y);
46: KSPGetIterationNumber(jac->ksp,&its);
47: jac->its += its;
48: return(0);
49: }
51: static PetscErrorCode PCApplyTranspose_KSP(PC pc,Vec x,Vec y)
52: {
53: PetscErrorCode ierr;
54: PetscInt its;
55: PC_KSP *jac = (PC_KSP*)pc->data;
58: if (jac->ksp->presolve) {
59: VecCopy(x,y);
60: KSPSolve(jac->ksp,y,y);
61: } else {
62: KSPSolveTranspose(jac->ksp,x,y);
63: }
64: KSPCheckSolve(jac->ksp,pc,y);
65: KSPGetIterationNumber(jac->ksp,&its);
66: jac->its += its;
67: return(0);
68: }
70: static PetscErrorCode PCSetUp_KSP(PC pc)
71: {
73: PC_KSP *jac = (PC_KSP*)pc->data;
74: Mat mat;
77: if (!jac->ksp) {
78: PCKSPCreateKSP_KSP(pc);
79: KSPSetFromOptions(jac->ksp);
80: }
81: if (pc->useAmat) mat = pc->mat;
82: else mat = pc->pmat;
83: KSPSetOperators(jac->ksp,mat,pc->pmat);
84: KSPSetUp(jac->ksp);
85: return(0);
86: }
88: /* Default destroy, if it has never been setup */
89: static PetscErrorCode PCReset_KSP(PC pc)
90: {
91: PC_KSP *jac = (PC_KSP*)pc->data;
95: KSPDestroy(&jac->ksp);
96: return(0);
97: }
99: static PetscErrorCode PCDestroy_KSP(PC pc)
100: {
101: PC_KSP *jac = (PC_KSP*)pc->data;
105: KSPDestroy(&jac->ksp);
106: PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",NULL);
107: PetscObjectComposeFunction((PetscObject)pc,"PCKSPSetKSP_C",NULL);
108: PetscFree(pc->data);
109: return(0);
110: }
112: static PetscErrorCode PCView_KSP(PC pc,PetscViewer viewer)
113: {
114: PC_KSP *jac = (PC_KSP*)pc->data;
116: PetscBool iascii;
119: if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
120: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
121: if (iascii) {
122: if (pc->useAmat) {
123: PetscViewerASCIIPrintf(viewer," Using Amat (not Pmat) as operator on inner solve\n");
124: }
125: PetscViewerASCIIPrintf(viewer," KSP and PC on KSP preconditioner follow\n");
126: PetscViewerASCIIPrintf(viewer," ---------------------------------\n");
127: }
128: PetscViewerASCIIPushTab(viewer);
129: KSPView(jac->ksp,viewer);
130: PetscViewerASCIIPopTab(viewer);
131: if (iascii) {
132: PetscViewerASCIIPrintf(viewer," ---------------------------------\n");
133: }
134: return(0);
135: }
137: static PetscErrorCode PCKSPSetKSP_KSP(PC pc,KSP ksp)
138: {
139: PC_KSP *jac = (PC_KSP*)pc->data;
143: PetscObjectReference((PetscObject)ksp);
144: KSPDestroy(&jac->ksp);
145: jac->ksp = ksp;
146: return(0);
147: }
149: /*@
150: PCKSPSetKSP - Sets the KSP context for a KSP PC.
152: Collective on PC
154: Input Parameter:
155: + pc - the preconditioner context
156: - ksp - the KSP solver
158: Notes:
159: The PC and the KSP must have the same communicator
161: Level: advanced
163: @*/
164: PetscErrorCode PCKSPSetKSP(PC pc,KSP ksp)
165: {
172: PetscTryMethod(pc,"PCKSPSetKSP_C",(PC,KSP),(pc,ksp));
173: return(0);
174: }
176: static PetscErrorCode PCKSPGetKSP_KSP(PC pc,KSP *ksp)
177: {
178: PC_KSP *jac = (PC_KSP*)pc->data;
182: if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
183: *ksp = jac->ksp;
184: return(0);
185: }
187: /*@
188: PCKSPGetKSP - Gets the KSP context for a KSP PC.
190: Not Collective but KSP returned is parallel if PC was parallel
192: Input Parameter:
193: . pc - the preconditioner context
195: Output Parameters:
196: . ksp - the KSP solver
198: Notes:
199: You must call KSPSetUp() before calling PCKSPGetKSP().
201: If the PC is not a PCKSP object it raises an error
203: Level: advanced
205: @*/
206: PetscErrorCode PCKSPGetKSP(PC pc,KSP *ksp)
207: {
213: PetscUseMethod(pc,"PCKSPGetKSP_C",(PC,KSP*),(pc,ksp));
214: return(0);
215: }
217: static PetscErrorCode PCSetFromOptions_KSP(PetscOptionItems *PetscOptionsObject,PC pc)
218: {
219: PC_KSP *jac = (PC_KSP*)pc->data;
223: PetscOptionsHead(PetscOptionsObject,"PC KSP options");
224: if (jac->ksp) {
225: KSPSetFromOptions(jac->ksp);
226: }
227: PetscOptionsTail();
228: return(0);
229: }
231: /* ----------------------------------------------------------------------------------*/
233: /*MC
234: PCKSP - Defines a preconditioner that can consist of any KSP solver.
235: This allows, for example, embedding a Krylov method inside a preconditioner.
237: Options Database Key:
238: . -pc_use_amat - use the matrix that defines the linear system, Amat as the matrix for the
239: inner solver, otherwise by default it uses the matrix used to construct
240: the preconditioner, Pmat (see PCSetOperators())
242: Level: intermediate
244: Notes:
245: Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
246: the incorrect answer) unless you use KSPFGMRES as the other Krylov method
248: Developer Notes:
249: If the outer Krylov method has a nonzero initial guess it will compute a new residual based on that initial guess
250: 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
251: 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
252: input (current residual) vector, the action of the preconditioner cannot depend also on some other vector (the "initial guess"). For
253: KSPRICHARDSON it is possible to provide a PCApplyRichardson_PCKSP() that short circuits returning to the KSP object at each iteration to compute the
254: residual, see for example PCApplyRichardson_SOR(). We do not implement a PCApplyRichardson_PCKSP() because (1) using a KSP directly inside a Richardson
255: is not an efficient algorithm anyways and (2) implementing it for its > 1 would essentially require that we implement Richardson (reimplementing the
256: Richardson code) inside the PCApplyRichardson_PCKSP() leading to duplicate code.
258: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
259: PCSHELL, PCCOMPOSITE, PCSetUseAmat(), PCKSPGetKSP()
261: M*/
263: PETSC_EXTERN PetscErrorCode PCCreate_KSP(PC pc)
264: {
266: PC_KSP *jac;
269: PetscNewLog(pc,&jac);
270: pc->data = (void*)jac;
272: PetscMemzero(pc->ops,sizeof(struct _PCOps));
273: pc->ops->apply = PCApply_KSP;
274: pc->ops->applytranspose = PCApplyTranspose_KSP;
275: pc->ops->setup = PCSetUp_KSP;
276: pc->ops->reset = PCReset_KSP;
277: pc->ops->destroy = PCDestroy_KSP;
278: pc->ops->setfromoptions = PCSetFromOptions_KSP;
279: pc->ops->view = PCView_KSP;
281: PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",PCKSPGetKSP_KSP);
282: PetscObjectComposeFunction((PetscObject)pc,"PCKSPSetKSP_C",PCKSPSetKSP_KSP);
283: return(0);
284: }