Actual source code: pcksp.c
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: {
12: const char *prefix;
13: PC_KSP *jac = (PC_KSP*)pc->data;
14: DM dm;
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: PCGetDM(pc, &dm);
23: if (dm) {
24: KSPSetDM(jac->ksp, dm);
25: KSPSetDMActive(jac->ksp, PETSC_FALSE);
26: }
27: return 0;
28: }
30: static PetscErrorCode PCApply_KSP(PC pc,Vec x,Vec y)
31: {
32: PetscInt its;
33: PC_KSP *jac = (PC_KSP*)pc->data;
35: if (jac->ksp->presolve) {
36: VecCopy(x,y);
37: KSPSolve(jac->ksp,y,y);
38: } else {
39: KSPSolve(jac->ksp,x,y);
40: }
41: KSPCheckSolve(jac->ksp,pc,y);
42: KSPGetIterationNumber(jac->ksp,&its);
43: jac->its += its;
44: return 0;
45: }
47: static PetscErrorCode PCMatApply_KSP(PC pc,Mat X,Mat Y)
48: {
49: PetscInt its;
50: PC_KSP *jac = (PC_KSP*)pc->data;
52: if (jac->ksp->presolve) {
53: MatCopy(X,Y,SAME_NONZERO_PATTERN);
54: KSPMatSolve(jac->ksp,Y,Y); /* TODO FIXME: this will fail since KSPMatSolve does not allow inplace solve yet */
55: } else {
56: KSPMatSolve(jac->ksp,X,Y);
57: }
58: KSPCheckSolve(jac->ksp,pc,NULL);
59: KSPGetIterationNumber(jac->ksp,&its);
60: jac->its += its;
61: return 0;
62: }
64: static PetscErrorCode PCApplyTranspose_KSP(PC pc,Vec x,Vec y)
65: {
66: PetscInt its;
67: PC_KSP *jac = (PC_KSP*)pc->data;
69: if (jac->ksp->presolve) {
70: VecCopy(x,y);
71: KSPSolve(jac->ksp,y,y);
72: } else {
73: KSPSolveTranspose(jac->ksp,x,y);
74: }
75: KSPCheckSolve(jac->ksp,pc,y);
76: KSPGetIterationNumber(jac->ksp,&its);
77: jac->its += its;
78: return 0;
79: }
81: static PetscErrorCode PCSetUp_KSP(PC pc)
82: {
83: PC_KSP *jac = (PC_KSP*)pc->data;
84: Mat mat;
86: if (!jac->ksp) {
87: PCKSPCreateKSP_KSP(pc);
88: KSPSetFromOptions(jac->ksp);
89: }
90: if (pc->useAmat) mat = pc->mat;
91: else mat = pc->pmat;
92: KSPSetOperators(jac->ksp,mat,pc->pmat);
93: KSPSetUp(jac->ksp);
94: return 0;
95: }
97: /* Default destroy, if it has never been setup */
98: static PetscErrorCode PCReset_KSP(PC pc)
99: {
100: PC_KSP *jac = (PC_KSP*)pc->data;
102: KSPDestroy(&jac->ksp);
103: return 0;
104: }
106: static PetscErrorCode PCDestroy_KSP(PC pc)
107: {
108: PC_KSP *jac = (PC_KSP*)pc->data;
110: KSPDestroy(&jac->ksp);
111: PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",NULL);
112: PetscObjectComposeFunction((PetscObject)pc,"PCKSPSetKSP_C",NULL);
113: PetscFree(pc->data);
114: return 0;
115: }
117: static PetscErrorCode PCView_KSP(PC pc,PetscViewer viewer)
118: {
119: PC_KSP *jac = (PC_KSP*)pc->data;
120: PetscBool iascii;
122: if (!jac->ksp) PCKSPCreateKSP_KSP(pc);
123: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
124: if (iascii) {
125: if (pc->useAmat) {
126: PetscViewerASCIIPrintf(viewer," Using Amat (not Pmat) as operator on inner solve\n");
127: }
128: PetscViewerASCIIPrintf(viewer," KSP and PC on KSP preconditioner follow\n");
129: PetscViewerASCIIPrintf(viewer," ---------------------------------\n");
130: }
131: PetscViewerASCIIPushTab(viewer);
132: KSPView(jac->ksp,viewer);
133: PetscViewerASCIIPopTab(viewer);
134: if (iascii) {
135: PetscViewerASCIIPrintf(viewer," ---------------------------------\n");
136: }
137: return 0;
138: }
140: static PetscErrorCode PCKSPSetKSP_KSP(PC pc,KSP ksp)
141: {
142: PC_KSP *jac = (PC_KSP*)pc->data;
144: PetscObjectReference((PetscObject)ksp);
145: KSPDestroy(&jac->ksp);
146: jac->ksp = ksp;
147: return 0;
148: }
150: /*@
151: PCKSPSetKSP - Sets the KSP context for a KSP PC.
153: Collective on PC
155: Input Parameters:
156: + pc - the preconditioner context
157: - ksp - the KSP solver
159: Notes:
160: The PC and the KSP must have the same communicator
162: Level: advanced
164: @*/
165: PetscErrorCode PCKSPSetKSP(PC pc,KSP ksp)
166: {
170: PetscTryMethod(pc,"PCKSPSetKSP_C",(PC,KSP),(pc,ksp));
171: return 0;
172: }
174: static PetscErrorCode PCKSPGetKSP_KSP(PC pc,KSP *ksp)
175: {
176: PC_KSP *jac = (PC_KSP*)pc->data;
178: if (!jac->ksp) PCKSPCreateKSP_KSP(pc);
179: *ksp = jac->ksp;
180: return 0;
181: }
183: /*@
184: PCKSPGetKSP - Gets the KSP context for a KSP PC.
186: Not Collective but KSP returned is parallel if PC was parallel
188: Input Parameter:
189: . pc - the preconditioner context
191: Output Parameters:
192: . ksp - the KSP solver
194: Notes:
195: You must call KSPSetUp() before calling PCKSPGetKSP().
197: If the PC is not a PCKSP object it raises an error
199: Level: advanced
201: @*/
202: PetscErrorCode PCKSPGetKSP(PC pc,KSP *ksp)
203: {
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;
214: PetscOptionsHead(PetscOptionsObject,"PC KSP options");
215: if (jac->ksp) {
216: KSPSetFromOptions(jac->ksp);
217: }
218: PetscOptionsTail();
219: return 0;
220: }
222: /* ----------------------------------------------------------------------------------*/
224: /*MC
225: PCKSP - Defines a preconditioner that can consist of any KSP solver.
226: This allows, for example, embedding a Krylov method inside a preconditioner.
228: Options Database Key:
229: . -pc_use_amat - use the matrix that defines the linear system, Amat as the matrix for the
230: inner solver, otherwise by default it uses the matrix used to construct
231: the preconditioner, Pmat (see PCSetOperators())
233: Level: intermediate
235: Notes:
236: The application of an inexact Krylov solve is a nonlinear operation. Thus, performing a solve with KSP is,
237: in general, a nonlinear operation, so PCKSP is in general a nonlinear preconditioner.
238: Thus, one can see divergence or an incorrect answer unless using a flexible Krylov method (e.g. KSPFGMRES, KSPGCR, or KSPFCG) for the outer Krylov solve.
240: Developer Notes:
241: If the outer Krylov method has a nonzero initial guess it will compute a new residual based on that initial guess
242: 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
243: 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
244: input (current residual) vector, the action of the preconditioner cannot depend also on some other vector (the "initial guess"). For
245: KSPRICHARDSON it is possible to provide a PCApplyRichardson_PCKSP() that short circuits returning to the KSP object at each iteration to compute the
246: residual, see for example PCApplyRichardson_SOR(). We do not implement a PCApplyRichardson_PCKSP() because (1) using a KSP directly inside a Richardson
247: is not an efficient algorithm anyways and (2) implementing it for its > 1 would essentially require that we implement Richardson (reimplementing the
248: Richardson code) inside the PCApplyRichardson_PCKSP() leading to duplicate code.
250: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
251: PCSHELL, PCCOMPOSITE, PCSetUseAmat(), PCKSPGetKSP()
253: M*/
255: PETSC_EXTERN PetscErrorCode PCCreate_KSP(PC pc)
256: {
257: PC_KSP *jac;
259: PetscNewLog(pc,&jac);
260: pc->data = (void*)jac;
262: PetscMemzero(pc->ops,sizeof(struct _PCOps));
263: pc->ops->apply = PCApply_KSP;
264: pc->ops->matapply = PCMatApply_KSP;
265: pc->ops->applytranspose = PCApplyTranspose_KSP;
266: pc->ops->setup = PCSetUp_KSP;
267: pc->ops->reset = PCReset_KSP;
268: pc->ops->destroy = PCDestroy_KSP;
269: pc->ops->setfromoptions = PCSetFromOptions_KSP;
270: pc->ops->view = PCView_KSP;
272: PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",PCKSPGetKSP_KSP);
273: PetscObjectComposeFunction((PetscObject)pc,"PCKSPSetKSP_C",PCKSPSetKSP_KSP);
274: return 0;
275: }