Actual source code: pcksp.c
petsc-3.10.4 2019-02-26
1: #include <petsc/private/pcimpl.h>
2: #include <petscksp.h>
4: typedef struct {
5: KSP ksp;
6: PetscInt its; /* total number of iterations KSP uses */
7: } PC_KSP;
9: static PetscErrorCode PCKSPCreateKSP_KSP(PC pc)
10: {
12: const char *prefix;
13: PC_KSP *jac = (PC_KSP*)pc->data;
16: KSPCreate(PetscObjectComm((PetscObject)pc),&jac->ksp);
17: KSPSetErrorIfNotConverged(jac->ksp,pc->erroriffailure);
18: PetscObjectIncrementTabLevel((PetscObject)jac->ksp,(PetscObject)pc,1);
19: PCGetOptionsPrefix(pc,&prefix);
20: KSPSetOptionsPrefix(jac->ksp,prefix);
21: KSPAppendOptionsPrefix(jac->ksp,"ksp_");
22: return(0);
23: }
25: static PetscErrorCode PCApply_KSP(PC pc,Vec x,Vec y)
26: {
27: PetscErrorCode ierr;
28: PetscInt its;
29: PC_KSP *jac = (PC_KSP*)pc->data;
30: KSPConvergedReason reason;
33: KSPSolve(jac->ksp,x,y);
34: KSPGetConvergedReason(jac->ksp,&reason);
35: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
36: pc->failedreason = PC_SUBPC_ERROR;
37: }
38: KSPGetIterationNumber(jac->ksp,&its);
39: jac->its += its;
40: return(0);
41: }
43: static PetscErrorCode PCApplyTranspose_KSP(PC pc,Vec x,Vec y)
44: {
45: PetscErrorCode ierr;
46: PetscInt its;
47: PC_KSP *jac = (PC_KSP*)pc->data;
48: KSPConvergedReason reason;
51: KSPSolveTranspose(jac->ksp,x,y);
52: KSPGetConvergedReason(jac->ksp,&reason);
53: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
54: pc->failedreason = PC_SUBPC_ERROR;
55: }
56: KSPGetIterationNumber(jac->ksp,&its);
57: jac->its += its;
58: return(0);
59: }
61: static PetscErrorCode PCSetUp_KSP(PC pc)
62: {
64: PC_KSP *jac = (PC_KSP*)pc->data;
65: Mat mat;
68: if (!jac->ksp) {
69: PCKSPCreateKSP_KSP(pc);
70: KSPSetFromOptions(jac->ksp);
71: }
72: if (pc->useAmat) mat = pc->mat;
73: else mat = pc->pmat;
74: KSPSetOperators(jac->ksp,mat,pc->pmat);
75: KSPSetUp(jac->ksp);
76: return(0);
77: }
79: /* Default destroy, if it has never been setup */
80: static PetscErrorCode PCReset_KSP(PC pc)
81: {
82: PC_KSP *jac = (PC_KSP*)pc->data;
86: KSPDestroy(&jac->ksp);
87: return(0);
88: }
90: static PetscErrorCode PCDestroy_KSP(PC pc)
91: {
92: PC_KSP *jac = (PC_KSP*)pc->data;
96: KSPDestroy(&jac->ksp);
97: PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",NULL);
98: PetscObjectComposeFunction((PetscObject)pc,"PCKSPSetKSP_C",NULL);
99: PetscFree(pc->data);
100: return(0);
101: }
103: static PetscErrorCode PCView_KSP(PC pc,PetscViewer viewer)
104: {
105: PC_KSP *jac = (PC_KSP*)pc->data;
107: PetscBool iascii;
110: if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
111: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
112: if (iascii) {
113: if (pc->useAmat) {
114: PetscViewerASCIIPrintf(viewer," Using Amat (not Pmat) as operator on inner solve\n");
115: }
116: PetscViewerASCIIPrintf(viewer," KSP and PC on KSP preconditioner follow\n");
117: PetscViewerASCIIPrintf(viewer," ---------------------------------\n");
118: }
119: PetscViewerASCIIPushTab(viewer);
120: KSPView(jac->ksp,viewer);
121: PetscViewerASCIIPopTab(viewer);
122: if (iascii) {
123: PetscViewerASCIIPrintf(viewer," ---------------------------------\n");
124: }
125: return(0);
126: }
128: static PetscErrorCode PCKSPSetKSP_KSP(PC pc,KSP ksp)
129: {
130: PC_KSP *jac = (PC_KSP*)pc->data;
134: PetscObjectReference((PetscObject)ksp);
135: KSPDestroy(&jac->ksp);
136: jac->ksp = ksp;
137: return(0);
138: }
140: /*@
141: PCKSPSetKSP - Sets the KSP context for a KSP PC.
143: Collective on PC
145: Input Parameter:
146: + pc - the preconditioner context
147: - ksp - the KSP solver
149: Notes:
150: The PC and the KSP must have the same communicator
152: Level: advanced
154: .keywords: PC, KSP, set, context
155: @*/
156: PetscErrorCode PCKSPSetKSP(PC pc,KSP ksp)
157: {
164: PetscTryMethod(pc,"PCKSPSetKSP_C",(PC,KSP),(pc,ksp));
165: return(0);
166: }
168: static PetscErrorCode PCKSPGetKSP_KSP(PC pc,KSP *ksp)
169: {
170: PC_KSP *jac = (PC_KSP*)pc->data;
174: if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
175: *ksp = jac->ksp;
176: return(0);
177: }
179: /*@
180: PCKSPGetKSP - Gets the KSP context for a KSP PC.
182: Not Collective but KSP returned is parallel if PC was parallel
184: Input Parameter:
185: . pc - the preconditioner context
187: Output Parameters:
188: . ksp - the KSP solver
190: Notes:
191: You must call KSPSetUp() before calling PCKSPGetKSP().
193: If the PC is not a PCKSP object it raises an error
195: Level: advanced
197: .keywords: PC, KSP, get, context
198: @*/
199: PetscErrorCode PCKSPGetKSP(PC pc,KSP *ksp)
200: {
206: PetscUseMethod(pc,"PCKSPGetKSP_C",(PC,KSP*),(pc,ksp));
207: return(0);
208: }
210: static PetscErrorCode PCSetFromOptions_KSP(PetscOptionItems *PetscOptionsObject,PC pc)
211: {
212: PC_KSP *jac = (PC_KSP*)pc->data;
216: PetscOptionsHead(PetscOptionsObject,"PC KSP options");
217: if (jac->ksp) {
218: KSPSetFromOptions(jac->ksp);
219: }
220: PetscOptionsTail();
221: return(0);
222: }
224: /* ----------------------------------------------------------------------------------*/
226: /*MC
227: PCKSP - Defines a preconditioner that can consist of any KSP solver.
228: This allows, for example, embedding a Krylov method inside a preconditioner.
230: Options Database Key:
231: . -pc_use_amat - use the matrix that defines the linear system, Amat as the matrix for the
232: inner solver, otherwise by default it uses the matrix used to construct
233: the preconditioner, Pmat (see PCSetOperators())
235: Level: intermediate
237: Concepts: inner iteration
239: Notes:
240: Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
241: the incorrect answer) unless you use KSPFGMRES as the other Krylov method
243: Developer Notes:
244: If the outer Krylov method has a nonzero initial guess it will compute a new residual based on that initial guess
245: 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
246: 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
247: input (current residual) vector, the action of the preconditioner cannot depend also on some other vector (the "initial guess"). For
248: KSPRICHARDSON it is possible to provide a PCApplyRichardson_PCKSP() that short circuits returning to the KSP object at each iteration to compute the
249: residual, see for example PCApplyRichardson_SOR(). We do not implement a PCApplyRichardson_PCKSP() because (1) using a KSP directly inside a Richardson
250: is not an efficient algorithm anyways and (2) implementing it for its > 1 would essentially require that we implement Richardson (reimplementing the
251: Richardson code) inside the PCApplyRichardson_PCKSP() leading to duplicate code.
253: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
254: PCSHELL, PCCOMPOSITE, PCSetUseAmat(), PCKSPGetKSP()
256: M*/
258: PETSC_EXTERN PetscErrorCode PCCreate_KSP(PC pc)
259: {
261: PC_KSP *jac;
264: PetscNewLog(pc,&jac);
265: pc->data = (void*)jac;
267: PetscMemzero(pc->ops,sizeof(struct _PCOps));
268: pc->ops->apply = PCApply_KSP;
269: pc->ops->applytranspose = PCApplyTranspose_KSP;
270: pc->ops->setup = PCSetUp_KSP;
271: pc->ops->reset = PCReset_KSP;
272: pc->ops->destroy = PCDestroy_KSP;
273: pc->ops->setfromoptions = PCSetFromOptions_KSP;
274: pc->ops->view = PCView_KSP;
276: PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",PCKSPGetKSP_KSP);
277: PetscObjectComposeFunction((PetscObject)pc,"PCKSPSetKSP_C",PCKSPSetKSP_KSP);
278: return(0);
279: }