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