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 PCSetUp_KSP(PC pc)
87: {
88: PC_KSP *jac = (PC_KSP *)pc->data;
89: Mat mat;
91: PetscFunctionBegin;
92: if (!jac->ksp) {
93: PetscCall(PCKSPCreateKSP_KSP(pc));
94: PetscCall(KSPSetFromOptions(jac->ksp));
95: }
96: if (pc->useAmat) mat = pc->mat;
97: else mat = pc->pmat;
98: PetscCall(KSPSetOperators(jac->ksp, mat, pc->pmat));
99: PetscCall(KSPSetUp(jac->ksp));
100: PetscFunctionReturn(PETSC_SUCCESS);
101: }
103: /* Default destroy, if it has never been setup */
104: static PetscErrorCode PCReset_KSP(PC pc)
105: {
106: PC_KSP *jac = (PC_KSP *)pc->data;
108: PetscFunctionBegin;
109: PetscCall(KSPDestroy(&jac->ksp));
110: PetscFunctionReturn(PETSC_SUCCESS);
111: }
113: static PetscErrorCode PCDestroy_KSP(PC pc)
114: {
115: PC_KSP *jac = (PC_KSP *)pc->data;
117: PetscFunctionBegin;
118: PetscCall(KSPDestroy(&jac->ksp));
119: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCKSPGetKSP_C", NULL));
120: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCKSPSetKSP_C", NULL));
121: PetscCall(PetscFree(pc->data));
122: PetscFunctionReturn(PETSC_SUCCESS);
123: }
125: static PetscErrorCode PCView_KSP(PC pc, PetscViewer viewer)
126: {
127: PC_KSP *jac = (PC_KSP *)pc->data;
128: PetscBool iascii;
130: PetscFunctionBegin;
131: if (!jac->ksp) PetscCall(PCKSPCreateKSP_KSP(pc));
132: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
133: if (iascii) {
134: if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, " Using Amat (not Pmat) as operator on inner solve\n"));
135: PetscCall(PetscViewerASCIIPrintf(viewer, " KSP and PC on KSP preconditioner follow\n"));
136: PetscCall(PetscViewerASCIIPrintf(viewer, " ---------------------------------\n"));
137: }
138: PetscCall(PetscViewerASCIIPushTab(viewer));
139: PetscCall(KSPView(jac->ksp, viewer));
140: PetscCall(PetscViewerASCIIPopTab(viewer));
141: if (iascii) PetscCall(PetscViewerASCIIPrintf(viewer, " ---------------------------------\n"));
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
145: static PetscErrorCode PCKSPSetKSP_KSP(PC pc, KSP ksp)
146: {
147: PC_KSP *jac = (PC_KSP *)pc->data;
149: PetscFunctionBegin;
150: PetscCall(PetscObjectReference((PetscObject)ksp));
151: PetscCall(KSPDestroy(&jac->ksp));
152: jac->ksp = ksp;
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: /*@
157: PCKSPSetKSP - Sets the `KSP` context for a `PCKSP`.
159: Collective
161: Input Parameters:
162: + pc - the preconditioner context
163: - ksp - the `KSP` solver
165: Level: advanced
167: Notes:
168: The `PC` and the `KSP` must have the same communicator
170: This would rarely be used, the standard usage is to call `PCKSPGetKSP()` and then change options on that `KSP`
172: .seealso: [](ch_ksp), `PCKSP`, `PCKSPGetKSP()`
173: @*/
174: PetscErrorCode PCKSPSetKSP(PC pc, KSP ksp)
175: {
176: PetscFunctionBegin;
179: PetscCheckSameComm(pc, 1, ksp, 2);
180: PetscTryMethod(pc, "PCKSPSetKSP_C", (PC, KSP), (pc, ksp));
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: static PetscErrorCode PCKSPGetKSP_KSP(PC pc, KSP *ksp)
185: {
186: PC_KSP *jac = (PC_KSP *)pc->data;
188: PetscFunctionBegin;
189: if (!jac->ksp) PetscCall(PCKSPCreateKSP_KSP(pc));
190: *ksp = jac->ksp;
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
194: /*@
195: PCKSPGetKSP - Gets the `KSP` context for a `PCKSP`.
197: Not Collective but ksp returned is parallel if pc was parallel
199: Input Parameter:
200: . pc - the preconditioner context
202: Output Parameter:
203: . ksp - the `KSP` solver
205: Note:
206: If the `PC` is not a `PCKSP` object it raises an error
208: Level: advanced
210: .seealso: [](ch_ksp), `PCKSP`, `PCKSPSetKSP()`
211: @*/
212: PetscErrorCode PCKSPGetKSP(PC pc, KSP *ksp)
213: {
214: PetscFunctionBegin;
216: PetscAssertPointer(ksp, 2);
217: PetscUseMethod(pc, "PCKSPGetKSP_C", (PC, KSP *), (pc, ksp));
218: PetscFunctionReturn(PETSC_SUCCESS);
219: }
221: static PetscErrorCode PCSetFromOptions_KSP(PC pc, PetscOptionItems *PetscOptionsObject)
222: {
223: PC_KSP *jac = (PC_KSP *)pc->data;
225: PetscFunctionBegin;
226: PetscOptionsHeadBegin(PetscOptionsObject, "PC KSP options");
227: if (!jac->ksp) PetscCall(PCKSPCreateKSP_KSP(pc));
228: PetscCall(KSPSetFromOptions(jac->ksp));
229: PetscOptionsHeadEnd();
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
233: /*MC
234: PCKSP - Defines a preconditioner as any `KSP` solver. This allows, for example, embedding a Krylov method inside a preconditioner.
236: Options Database Key:
237: . -pc_use_amat - use the matrix that defines the linear system, Amat as the matrix for the
238: inner solver, otherwise by default it uses the matrix used to construct
239: the preconditioner, Pmat (see `PCSetOperators()`)
241: Level: intermediate
243: Note:
244: The application of an inexact Krylov solve is a nonlinear operation. Thus, performing a solve with `KSP` is,
245: in general, a nonlinear operation, so `PCKSP` is in general a nonlinear preconditioner.
246: 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.
248: Developer Note:
249: If the outer Krylov method has a nonzero initial guess it will compute a new residual based on that initial guess
250: 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
251: 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
252: input (current residual) vector, the action of the preconditioner cannot depend also on some other vector (the "initial guess"). For
253: `KSPRICHARDSON` it is possible to provide a `PCApplyRichardson_PCKSP()` that short circuits returning to the `KSP` object at each iteration to compute the
254: residual, see for example `PCApplyRichardson_SOR()`. We do not implement a `PCApplyRichardson_PCKSP()` because (1) using a `KSP` directly inside a Richardson
255: is not an efficient algorithm anyways and (2) implementing it for its > 1 would essentially require that we implement Richardson (reimplementing the
256: Richardson code) inside the `PCApplyRichardson_PCKSP()` leading to duplicate code.
258: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
259: `PCSHELL`, `PCCOMPOSITE`, `PCSetUseAmat()`, `PCKSPGetKSP()`, `KSPFGMRES`, `KSPGCR`, `KSPFCG`
260: M*/
262: PETSC_EXTERN PetscErrorCode PCCreate_KSP(PC pc)
263: {
264: PC_KSP *jac;
266: PetscFunctionBegin;
267: PetscCall(PetscNew(&jac));
268: pc->data = (void *)jac;
270: PetscCall(PetscMemzero(pc->ops, sizeof(struct _PCOps)));
271: pc->ops->apply = PCApply_KSP;
272: pc->ops->matapply = PCMatApply_KSP;
273: pc->ops->applytranspose = PCApplyTranspose_KSP;
274: pc->ops->setup = PCSetUp_KSP;
275: pc->ops->reset = PCReset_KSP;
276: pc->ops->destroy = PCDestroy_KSP;
277: pc->ops->setfromoptions = PCSetFromOptions_KSP;
278: pc->ops->view = PCView_KSP;
280: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCKSPGetKSP_C", PCKSPGetKSP_KSP));
281: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCKSPSetKSP_C", PCKSPSetKSP_KSP));
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }