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: PetscFunctionBegin;
17: PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->ksp));
18: PetscCall(KSPSetNestLevel(jac->ksp, pc->kspnestlevel));
19: PetscCall(KSPSetErrorIfNotConverged(jac->ksp, pc->erroriffailure));
20: PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->ksp, (PetscObject)pc, 1));
21: PetscCall(PCGetOptionsPrefix(pc, &prefix));
22: PetscCall(KSPSetOptionsPrefix(jac->ksp, prefix));
23: PetscCall(KSPAppendOptionsPrefix(jac->ksp, "ksp_"));
24: PetscCall(PCGetDM(pc, &dm));
25: if (dm) {
26: PetscCall(KSPSetDM(jac->ksp, dm));
27: PetscCall(KSPSetDMActive(jac->ksp, PETSC_FALSE));
28: }
29: PetscFunctionReturn(PETSC_SUCCESS);
30: }
32: static PetscErrorCode PCApply_KSP(PC pc, Vec x, Vec y)
33: {
34: PetscInt its;
35: PC_KSP *jac = (PC_KSP *)pc->data;
37: PetscFunctionBegin;
38: if (jac->ksp->presolve) {
39: PetscCall(VecCopy(x, y));
40: PetscCall(KSPSolve(jac->ksp, y, y));
41: } else {
42: PetscCall(KSPSolve(jac->ksp, x, y));
43: }
44: PetscCall(KSPCheckSolve(jac->ksp, pc, y));
45: PetscCall(KSPGetIterationNumber(jac->ksp, &its));
46: jac->its += its;
47: PetscFunctionReturn(PETSC_SUCCESS);
48: }
50: static PetscErrorCode PCMatApply_KSP(PC pc, Mat X, Mat Y)
51: {
52: PetscInt its;
53: PC_KSP *jac = (PC_KSP *)pc->data;
55: PetscFunctionBegin;
56: if (jac->ksp->presolve) {
57: PetscCall(MatCopy(X, Y, SAME_NONZERO_PATTERN));
58: PetscCall(KSPMatSolve(jac->ksp, Y, Y)); /* TODO FIXME: this will fail since KSPMatSolve() does not allow inplace solve yet */
59: } else {
60: PetscCall(KSPMatSolve(jac->ksp, X, Y));
61: }
62: PetscCall(KSPCheckSolve(jac->ksp, pc, NULL));
63: PetscCall(KSPGetIterationNumber(jac->ksp, &its));
64: jac->its += its;
65: PetscFunctionReturn(PETSC_SUCCESS);
66: }
68: static PetscErrorCode PCApplyTranspose_KSP(PC pc, Vec x, Vec y)
69: {
70: PetscInt its;
71: PC_KSP *jac = (PC_KSP *)pc->data;
73: PetscFunctionBegin;
74: if (jac->ksp->presolve) {
75: PetscCall(VecCopy(x, y));
76: PetscCall(KSPSolve(jac->ksp, y, y));
77: } else {
78: PetscCall(KSPSolveTranspose(jac->ksp, x, y));
79: }
80: PetscCall(KSPCheckSolve(jac->ksp, pc, y));
81: PetscCall(KSPGetIterationNumber(jac->ksp, &its));
82: jac->its += its;
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: static PetscErrorCode PCMatApplyTranspose_KSP(PC pc, Mat X, Mat Y)
87: {
88: PetscInt its;
89: PC_KSP *jac = (PC_KSP *)pc->data;
91: PetscFunctionBegin;
92: if (jac->ksp->presolve) {
93: PetscCall(MatCopy(X, Y, SAME_NONZERO_PATTERN));
94: PetscCall(KSPMatSolveTranspose(jac->ksp, Y, Y)); /* TODO FIXME: this will fail since KSPMatSolveTranspose() does not allow inplace solve yet */
95: } else {
96: PetscCall(KSPMatSolveTranspose(jac->ksp, X, Y));
97: }
98: PetscCall(KSPCheckSolve(jac->ksp, pc, NULL));
99: PetscCall(KSPGetIterationNumber(jac->ksp, &its));
100: jac->its += its;
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: static PetscErrorCode PCSetUp_KSP(PC pc)
105: {
106: PC_KSP *jac = (PC_KSP *)pc->data;
107: Mat mat;
109: PetscFunctionBegin;
110: if (!jac->ksp) {
111: PetscCall(PCKSPCreateKSP_KSP(pc));
112: PetscCall(KSPSetFromOptions(jac->ksp));
113: }
114: if (pc->useAmat) mat = pc->mat;
115: else mat = pc->pmat;
116: PetscCall(KSPSetOperators(jac->ksp, mat, pc->pmat));
117: PetscCall(KSPSetUp(jac->ksp));
118: PetscFunctionReturn(PETSC_SUCCESS);
119: }
121: /* Default destroy, if it has never been setup */
122: static PetscErrorCode PCReset_KSP(PC pc)
123: {
124: PC_KSP *jac = (PC_KSP *)pc->data;
126: PetscFunctionBegin;
127: PetscCall(KSPDestroy(&jac->ksp));
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: static PetscErrorCode PCDestroy_KSP(PC pc)
132: {
133: PC_KSP *jac = (PC_KSP *)pc->data;
135: PetscFunctionBegin;
136: PetscCall(KSPDestroy(&jac->ksp));
137: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCKSPGetKSP_C", NULL));
138: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCKSPSetKSP_C", NULL));
139: PetscCall(PetscFree(pc->data));
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: static PetscErrorCode PCView_KSP(PC pc, PetscViewer viewer)
144: {
145: PC_KSP *jac = (PC_KSP *)pc->data;
146: PetscBool iascii;
148: PetscFunctionBegin;
149: if (!jac->ksp) PetscCall(PCKSPCreateKSP_KSP(pc));
150: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
151: if (iascii) {
152: if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, " Using Amat (not Pmat) as operator on inner solve\n"));
153: PetscCall(PetscViewerASCIIPrintf(viewer, " KSP and PC on KSP preconditioner follow\n"));
154: PetscCall(PetscViewerASCIIPrintf(viewer, " ---------------------------------\n"));
155: }
156: PetscCall(PetscViewerASCIIPushTab(viewer));
157: PetscCall(KSPView(jac->ksp, viewer));
158: PetscCall(PetscViewerASCIIPopTab(viewer));
159: if (iascii) PetscCall(PetscViewerASCIIPrintf(viewer, " ---------------------------------\n"));
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: static PetscErrorCode PCKSPSetKSP_KSP(PC pc, KSP ksp)
164: {
165: PC_KSP *jac = (PC_KSP *)pc->data;
167: PetscFunctionBegin;
168: PetscCall(PetscObjectReference((PetscObject)ksp));
169: PetscCall(KSPDestroy(&jac->ksp));
170: jac->ksp = ksp;
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: /*@
175: PCKSPSetKSP - Sets the `KSP` context for a `PCKSP`.
177: Collective
179: Input Parameters:
180: + pc - the preconditioner context
181: - ksp - the `KSP` solver
183: Level: advanced
185: Notes:
186: The `PC` and the `KSP` must have the same communicator
188: This would rarely be used, the standard usage is to call `PCKSPGetKSP()` and then change options on that `KSP`
190: .seealso: [](ch_ksp), `PCKSP`, `PCKSPGetKSP()`
191: @*/
192: PetscErrorCode PCKSPSetKSP(PC pc, KSP ksp)
193: {
194: PetscFunctionBegin;
197: PetscCheckSameComm(pc, 1, ksp, 2);
198: PetscTryMethod(pc, "PCKSPSetKSP_C", (PC, KSP), (pc, ksp));
199: PetscFunctionReturn(PETSC_SUCCESS);
200: }
202: static PetscErrorCode PCKSPGetKSP_KSP(PC pc, KSP *ksp)
203: {
204: PC_KSP *jac = (PC_KSP *)pc->data;
206: PetscFunctionBegin;
207: if (!jac->ksp) PetscCall(PCKSPCreateKSP_KSP(pc));
208: *ksp = jac->ksp;
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: /*@
213: PCKSPGetKSP - Gets the `KSP` context for a `PCKSP`.
215: Not Collective but ksp returned is parallel if pc was parallel
217: Input Parameter:
218: . pc - the preconditioner context
220: Output Parameter:
221: . ksp - the `KSP` solver
223: Note:
224: If the `PC` is not a `PCKSP` object it raises an error
226: Level: advanced
228: .seealso: [](ch_ksp), `PCKSP`, `PCKSPSetKSP()`
229: @*/
230: PetscErrorCode PCKSPGetKSP(PC pc, KSP *ksp)
231: {
232: PetscFunctionBegin;
234: PetscAssertPointer(ksp, 2);
235: PetscUseMethod(pc, "PCKSPGetKSP_C", (PC, KSP *), (pc, ksp));
236: PetscFunctionReturn(PETSC_SUCCESS);
237: }
239: static PetscErrorCode PCSetFromOptions_KSP(PC pc, PetscOptionItems PetscOptionsObject)
240: {
241: PC_KSP *jac = (PC_KSP *)pc->data;
243: PetscFunctionBegin;
244: PetscOptionsHeadBegin(PetscOptionsObject, "PC KSP options");
245: if (!jac->ksp) PetscCall(PCKSPCreateKSP_KSP(pc));
246: PetscCall(KSPSetFromOptions(jac->ksp));
247: PetscOptionsHeadEnd();
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
251: /*MC
252: PCKSP - Defines a preconditioner as any `KSP` solver. This allows, for example, embedding a Krylov method inside a preconditioner.
254: Options Database Key:
255: . -pc_use_amat - use the matrix that defines the linear system, Amat as the matrix for the
256: inner solver, otherwise by default it uses the matrix used to construct
257: the preconditioner, Pmat (see `PCSetOperators()`)
259: Level: intermediate
261: Note:
262: The application of an inexact Krylov solve is a nonlinear operation. Thus, performing a solve with `KSP` is,
263: in general, a nonlinear operation, so `PCKSP` is in general a nonlinear preconditioner.
264: 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.
266: Developer Note:
267: If the outer Krylov method has a nonzero initial guess it will compute a new residual based on that initial guess
268: 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
269: 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
270: input (current residual) vector, the action of the preconditioner cannot depend also on some other vector (the "initial guess"). For
271: `KSPRICHARDSON` it is possible to provide a `PCApplyRichardson_PCKSP()` that short circuits returning to the `KSP` object at each iteration to compute the
272: residual, see for example `PCApplyRichardson_SOR()`. We do not implement a `PCApplyRichardson_PCKSP()` because (1) using a `KSP` directly inside a Richardson
273: is not an efficient algorithm anyways and (2) implementing it for its > 1 would essentially require that we implement Richardson (reimplementing the
274: Richardson code) inside the `PCApplyRichardson_PCKSP()` leading to duplicate code.
276: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
277: `PCSHELL`, `PCCOMPOSITE`, `PCSetUseAmat()`, `PCKSPGetKSP()`, `KSPFGMRES`, `KSPGCR`, `KSPFCG`
278: M*/
280: PETSC_EXTERN PetscErrorCode PCCreate_KSP(PC pc)
281: {
282: PC_KSP *jac;
284: PetscFunctionBegin;
285: PetscCall(PetscNew(&jac));
286: pc->data = (void *)jac;
288: PetscCall(PetscMemzero(pc->ops, sizeof(struct _PCOps)));
289: pc->ops->apply = PCApply_KSP;
290: pc->ops->matapply = PCMatApply_KSP;
291: pc->ops->applytranspose = PCApplyTranspose_KSP;
292: pc->ops->matapplytranspose = PCMatApplyTranspose_KSP;
293: pc->ops->setup = PCSetUp_KSP;
294: pc->ops->reset = PCReset_KSP;
295: pc->ops->destroy = PCDestroy_KSP;
296: pc->ops->setfromoptions = PCSetFromOptions_KSP;
297: pc->ops->view = PCView_KSP;
299: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCKSPGetKSP_C", PCKSPGetKSP_KSP));
300: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCKSPSetKSP_C", PCKSPSetKSP_KSP));
301: PetscFunctionReturn(PETSC_SUCCESS);
302: }