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) PetscViewerASCIIPrintf(viewer, " Using Amat (not Pmat) as operator on inner solve\n");
126: PetscViewerASCIIPrintf(viewer, " KSP and PC on KSP preconditioner follow\n");
127: PetscViewerASCIIPrintf(viewer, " ---------------------------------\n");
128: }
129: PetscViewerASCIIPushTab(viewer);
130: KSPView(jac->ksp, viewer);
131: PetscViewerASCIIPopTab(viewer);
132: if (iascii) PetscViewerASCIIPrintf(viewer, " ---------------------------------\n");
133: return 0;
134: }
136: static PetscErrorCode PCKSPSetKSP_KSP(PC pc, KSP ksp)
137: {
138: PC_KSP *jac = (PC_KSP *)pc->data;
140: PetscObjectReference((PetscObject)ksp);
141: KSPDestroy(&jac->ksp);
142: jac->ksp = ksp;
143: return 0;
144: }
146: /*@
147: PCKSPSetKSP - Sets the `KSP` context for a `PCKSP`.
149: Collective
151: Input Parameters:
152: + pc - the preconditioner context
153: - ksp - the `KSP` solver
155: Level: advanced
157: Notes:
158: The `PC` and the `KSP` must have the same communicator
160: This would rarely be used, the standard usage is to call `PCKSPGetKSP()` and then change options on that `KSP`
162: .seealso: `PCKSP`, `PCKSPGetKSP()`
163: @*/
164: PetscErrorCode PCKSPSetKSP(PC pc, KSP ksp)
165: {
169: PetscTryMethod(pc, "PCKSPSetKSP_C", (PC, KSP), (pc, ksp));
170: return 0;
171: }
173: static PetscErrorCode PCKSPGetKSP_KSP(PC pc, KSP *ksp)
174: {
175: 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 `PCKSP`.
185: Not Collective but ksp returned is parallel if pc was parallel
187: Input Parameter:
188: . pc - the preconditioner context
190: Output Parameter:
191: . ksp - the `KSP` solver
193: Note:
194: If the `PC` is not a `PCKSP` object it raises an error
196: Level: advanced
198: .seealso: `PCKSP`, `PCKSPSetKSP()`
199: @*/
200: PetscErrorCode PCKSPGetKSP(PC pc, KSP *ksp)
201: {
204: PetscUseMethod(pc, "PCKSPGetKSP_C", (PC, KSP *), (pc, ksp));
205: return 0;
206: }
208: static PetscErrorCode PCSetFromOptions_KSP(PC pc, PetscOptionItems *PetscOptionsObject)
209: {
210: PC_KSP *jac = (PC_KSP *)pc->data;
212: PetscOptionsHeadBegin(PetscOptionsObject, "PC KSP options");
213: if (jac->ksp) KSPSetFromOptions(jac->ksp);
214: PetscOptionsHeadEnd();
215: return 0;
216: }
218: /*MC
219: PCKSP - Defines a preconditioner as any `KSP` solver.
220: This allows, for example, embedding a Krylov method inside a preconditioner.
222: Options Database Key:
223: . -pc_use_amat - use the matrix that defines the linear system, Amat as the matrix for the
224: inner solver, otherwise by default it uses the matrix used to construct
225: the preconditioner, Pmat (see `PCSetOperators()`)
227: Level: intermediate
229: Note:
230: The application of an inexact Krylov solve is a nonlinear operation. Thus, performing a solve with `KSP` is,
231: in general, a nonlinear operation, so `PCKSP` is in general a nonlinear preconditioner.
232: 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.
234: Developer Note:
235: If the outer Krylov method has a nonzero initial guess it will compute a new residual based on that initial guess
236: 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
237: except Richardson this is necessary since Krylov methods, even the flexible ones, need to "see" the result of the action of the preconditioner on the
238: input (current residual) vector, the action of the preconditioner cannot depend also on some other vector (the "initial guess"). For
239: `KSPRICHARDSON` it is possible to provide a `PCApplyRichardson_PCKSP()` that short circuits returning to the `KSP` object at each iteration to compute the
240: residual, see for example `PCApplyRichardson_SOR()`. We do not implement a `PCApplyRichardson_PCKSP()` because (1) using a `KSP` directly inside a Richardson
241: is not an efficient algorithm anyways and (2) implementing it for its > 1 would essentially require that we implement Richardson (reimplementing the
242: Richardson code) inside the `PCApplyRichardson_PCKSP()` leading to duplicate code.
244: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
245: `PCSHELL`, `PCCOMPOSITE`, `PCSetUseAmat()`, `PCKSPGetKSP()`, `KSPFGMRES`, `KSPGCR`, `KSPFCG`
246: M*/
248: PETSC_EXTERN PetscErrorCode PCCreate_KSP(PC pc)
249: {
250: PC_KSP *jac;
252: PetscNew(&jac);
253: pc->data = (void *)jac;
255: PetscMemzero(pc->ops, sizeof(struct _PCOps));
256: pc->ops->apply = PCApply_KSP;
257: pc->ops->matapply = PCMatApply_KSP;
258: pc->ops->applytranspose = PCApplyTranspose_KSP;
259: pc->ops->setup = PCSetUp_KSP;
260: pc->ops->reset = PCReset_KSP;
261: pc->ops->destroy = PCDestroy_KSP;
262: pc->ops->setfromoptions = PCSetFromOptions_KSP;
263: pc->ops->view = PCView_KSP;
265: PetscObjectComposeFunction((PetscObject)pc, "PCKSPGetKSP_C", PCKSPGetKSP_KSP);
266: PetscObjectComposeFunction((PetscObject)pc, "PCKSPSetKSP_C", PCKSPSetKSP_KSP);
267: return 0;
268: }