Actual source code: pcksp.c
petsc-3.11.3 2019-06-26
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;
17: KSPCreate(PetscObjectComm((PetscObject)pc),&jac->ksp);
18: KSPSetErrorIfNotConverged(jac->ksp,pc->erroriffailure);
19: PetscObjectIncrementTabLevel((PetscObject)jac->ksp,(PetscObject)pc,1);
20: PCGetOptionsPrefix(pc,&prefix);
21: KSPSetOptionsPrefix(jac->ksp,prefix);
22: KSPAppendOptionsPrefix(jac->ksp,"ksp_");
23: return(0);
24: }
26: static PetscErrorCode PCApply_KSP(PC pc,Vec x,Vec y)
27: {
28: PetscErrorCode ierr;
29: PetscInt its;
30: PC_KSP *jac = (PC_KSP*)pc->data;
33: if (jac->ksp->presolve) {
34: VecCopy(x,y);
35: KSPSolve(jac->ksp,y,y);
36: } else {
37: KSPSolve(jac->ksp,x,y);
38: }
39: KSPCheckSolve(jac->ksp,pc,y);
40: KSPGetIterationNumber(jac->ksp,&its);
41: jac->its += its;
42: return(0);
43: }
45: static PetscErrorCode PCApplyTranspose_KSP(PC pc,Vec x,Vec y)
46: {
47: PetscErrorCode ierr;
48: PetscInt its;
49: PC_KSP *jac = (PC_KSP*)pc->data;
52: if (jac->ksp->presolve) {
53: VecCopy(x,y);
54: KSPSolve(jac->ksp,y,y);
55: } else {
56: KSPSolveTranspose(jac->ksp,x,y);
57: }
58: KSPCheckSolve(jac->ksp,pc,y);
59: KSPGetIterationNumber(jac->ksp,&its);
60: jac->its += its;
61: return(0);
62: }
64: static PetscErrorCode PCSetUp_KSP(PC pc)
65: {
67: PC_KSP *jac = (PC_KSP*)pc->data;
68: Mat mat;
71: if (!jac->ksp) {
72: PCKSPCreateKSP_KSP(pc);
73: KSPSetFromOptions(jac->ksp);
74: }
75: if (pc->useAmat) mat = pc->mat;
76: else mat = pc->pmat;
77: KSPSetOperators(jac->ksp,mat,pc->pmat);
78: KSPSetUp(jac->ksp);
79: return(0);
80: }
82: /* Default destroy, if it has never been setup */
83: static PetscErrorCode PCReset_KSP(PC pc)
84: {
85: PC_KSP *jac = (PC_KSP*)pc->data;
89: KSPDestroy(&jac->ksp);
90: return(0);
91: }
93: static PetscErrorCode PCDestroy_KSP(PC pc)
94: {
95: PC_KSP *jac = (PC_KSP*)pc->data;
99: KSPDestroy(&jac->ksp);
100: PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",NULL);
101: PetscObjectComposeFunction((PetscObject)pc,"PCKSPSetKSP_C",NULL);
102: PetscFree(pc->data);
103: return(0);
104: }
106: static PetscErrorCode PCView_KSP(PC pc,PetscViewer viewer)
107: {
108: PC_KSP *jac = (PC_KSP*)pc->data;
110: PetscBool iascii;
113: if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
114: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
115: if (iascii) {
116: if (pc->useAmat) {
117: PetscViewerASCIIPrintf(viewer," Using Amat (not Pmat) as operator on inner solve\n");
118: }
119: PetscViewerASCIIPrintf(viewer," KSP and PC on KSP preconditioner follow\n");
120: PetscViewerASCIIPrintf(viewer," ---------------------------------\n");
121: }
122: PetscViewerASCIIPushTab(viewer);
123: KSPView(jac->ksp,viewer);
124: PetscViewerASCIIPopTab(viewer);
125: if (iascii) {
126: PetscViewerASCIIPrintf(viewer," ---------------------------------\n");
127: }
128: return(0);
129: }
131: static PetscErrorCode PCKSPSetKSP_KSP(PC pc,KSP ksp)
132: {
133: PC_KSP *jac = (PC_KSP*)pc->data;
137: PetscObjectReference((PetscObject)ksp);
138: KSPDestroy(&jac->ksp);
139: jac->ksp = ksp;
140: return(0);
141: }
143: /*@
144: PCKSPSetKSP - Sets the KSP context for a KSP PC.
146: Collective on PC
148: Input Parameter:
149: + pc - the preconditioner context
150: - ksp - the KSP solver
152: Notes:
153: The PC and the KSP must have the same communicator
155: Level: advanced
157: .keywords: PC, KSP, set, context
158: @*/
159: PetscErrorCode PCKSPSetKSP(PC pc,KSP ksp)
160: {
167: PetscTryMethod(pc,"PCKSPSetKSP_C",(PC,KSP),(pc,ksp));
168: return(0);
169: }
171: static PetscErrorCode PCKSPGetKSP_KSP(PC pc,KSP *ksp)
172: {
173: PC_KSP *jac = (PC_KSP*)pc->data;
177: if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
178: *ksp = jac->ksp;
179: return(0);
180: }
182: /*@
183: PCKSPGetKSP - Gets the KSP context for a KSP PC.
185: Not Collective but KSP returned is parallel if PC was parallel
187: Input Parameter:
188: . pc - the preconditioner context
190: Output Parameters:
191: . ksp - the KSP solver
193: Notes:
194: You must call KSPSetUp() before calling PCKSPGetKSP().
196: If the PC is not a PCKSP object it raises an error
198: Level: advanced
200: .keywords: PC, KSP, get, context
201: @*/
202: PetscErrorCode PCKSPGetKSP(PC pc,KSP *ksp)
203: {
209: PetscUseMethod(pc,"PCKSPGetKSP_C",(PC,KSP*),(pc,ksp));
210: return(0);
211: }
213: static PetscErrorCode PCSetFromOptions_KSP(PetscOptionItems *PetscOptionsObject,PC pc)
214: {
215: PC_KSP *jac = (PC_KSP*)pc->data;
219: PetscOptionsHead(PetscOptionsObject,"PC KSP options");
220: if (jac->ksp) {
221: KSPSetFromOptions(jac->ksp);
222: }
223: PetscOptionsTail();
224: return(0);
225: }
227: /* ----------------------------------------------------------------------------------*/
229: /*MC
230: PCKSP - Defines a preconditioner that can consist of any KSP solver.
231: This allows, for example, embedding a Krylov method inside a preconditioner.
233: Options Database Key:
234: . -pc_use_amat - use the matrix that defines the linear system, Amat as the matrix for the
235: inner solver, otherwise by default it uses the matrix used to construct
236: the preconditioner, Pmat (see PCSetOperators())
238: Level: intermediate
240: Concepts: inner iteration
242: Notes:
243: Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
244: the incorrect answer) unless you use KSPFGMRES as the other Krylov method
246: Developer Notes:
247: If the outer Krylov method has a nonzero initial guess it will compute a new residual based on that initial guess
248: 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
249: 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
250: input (current residual) vector, the action of the preconditioner cannot depend also on some other vector (the "initial guess"). For
251: KSPRICHARDSON it is possible to provide a PCApplyRichardson_PCKSP() that short circuits returning to the KSP object at each iteration to compute the
252: residual, see for example PCApplyRichardson_SOR(). We do not implement a PCApplyRichardson_PCKSP() because (1) using a KSP directly inside a Richardson
253: is not an efficient algorithm anyways and (2) implementing it for its > 1 would essentially require that we implement Richardson (reimplementing the
254: Richardson code) inside the PCApplyRichardson_PCKSP() leading to duplicate code.
256: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
257: PCSHELL, PCCOMPOSITE, PCSetUseAmat(), PCKSPGetKSP()
259: M*/
261: PETSC_EXTERN PetscErrorCode PCCreate_KSP(PC pc)
262: {
264: PC_KSP *jac;
267: PetscNewLog(pc,&jac);
268: pc->data = (void*)jac;
270: PetscMemzero(pc->ops,sizeof(struct _PCOps));
271: pc->ops->apply = PCApply_KSP;
272: pc->ops->applytranspose = PCApplyTranspose_KSP;
273: pc->ops->setup = PCSetUp_KSP;
274: pc->ops->reset = PCReset_KSP;
275: pc->ops->destroy = PCDestroy_KSP;
276: pc->ops->setfromoptions = PCSetFromOptions_KSP;
277: pc->ops->view = PCView_KSP;
279: PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",PCKSPGetKSP_KSP);
280: PetscObjectComposeFunction((PetscObject)pc,"PCKSPSetKSP_C",PCKSPSetKSP_KSP);
281: return(0);
282: }