Actual source code: pcksp.c
petsc-3.14.0 2020-09-29
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: {
13: const char *prefix;
14: PC_KSP *jac = (PC_KSP*)pc->data;
15: DM dm;
18: KSPCreate(PetscObjectComm((PetscObject)pc),&jac->ksp);
19: KSPSetErrorIfNotConverged(jac->ksp,pc->erroriffailure);
20: PetscObjectIncrementTabLevel((PetscObject)jac->ksp,(PetscObject)pc,1);
21: PCGetOptionsPrefix(pc,&prefix);
22: KSPSetOptionsPrefix(jac->ksp,prefix);
23: KSPAppendOptionsPrefix(jac->ksp,"ksp_");
24: PCGetDM(pc, &dm);
25: if (dm) {
26: KSPSetDM(jac->ksp, dm);
27: KSPSetDMActive(jac->ksp, PETSC_FALSE);
28: }
29: return(0);
30: }
32: static PetscErrorCode PCApply_KSP(PC pc,Vec x,Vec y)
33: {
34: PetscErrorCode ierr;
35: PetscInt its;
36: PC_KSP *jac = (PC_KSP*)pc->data;
39: if (jac->ksp->presolve) {
40: VecCopy(x,y);
41: KSPSolve(jac->ksp,y,y);
42: } else {
43: KSPSolve(jac->ksp,x,y);
44: }
45: KSPCheckSolve(jac->ksp,pc,y);
46: KSPGetIterationNumber(jac->ksp,&its);
47: jac->its += its;
48: return(0);
49: }
51: static PetscErrorCode PCMatApply_KSP(PC pc,Mat X,Mat Y)
52: {
53: PetscErrorCode ierr;
54: PetscInt its;
55: PC_KSP *jac = (PC_KSP*)pc->data;
58: if (jac->ksp->presolve) {
59: MatCopy(X,Y,SAME_NONZERO_PATTERN);
60: KSPMatSolve(jac->ksp,Y,Y); /* TODO FIXME: this will fail since KSPMatSolve does not allow inplace solve yet */
61: } else {
62: KSPMatSolve(jac->ksp,X,Y);
63: }
64: KSPCheckSolve(jac->ksp,pc,NULL);
65: KSPGetIterationNumber(jac->ksp,&its);
66: jac->its += its;
67: return(0);
68: }
70: static PetscErrorCode PCApplyTranspose_KSP(PC pc,Vec x,Vec y)
71: {
72: PetscErrorCode ierr;
73: PetscInt its;
74: PC_KSP *jac = (PC_KSP*)pc->data;
77: if (jac->ksp->presolve) {
78: VecCopy(x,y);
79: KSPSolve(jac->ksp,y,y);
80: } else {
81: KSPSolveTranspose(jac->ksp,x,y);
82: }
83: KSPCheckSolve(jac->ksp,pc,y);
84: KSPGetIterationNumber(jac->ksp,&its);
85: jac->its += its;
86: return(0);
87: }
89: static PetscErrorCode PCSetUp_KSP(PC pc)
90: {
92: PC_KSP *jac = (PC_KSP*)pc->data;
93: Mat mat;
96: if (!jac->ksp) {
97: PCKSPCreateKSP_KSP(pc);
98: KSPSetFromOptions(jac->ksp);
99: }
100: if (pc->useAmat) mat = pc->mat;
101: else mat = pc->pmat;
102: KSPSetOperators(jac->ksp,mat,pc->pmat);
103: KSPSetUp(jac->ksp);
104: return(0);
105: }
107: /* Default destroy, if it has never been setup */
108: static PetscErrorCode PCReset_KSP(PC pc)
109: {
110: PC_KSP *jac = (PC_KSP*)pc->data;
114: KSPDestroy(&jac->ksp);
115: return(0);
116: }
118: static PetscErrorCode PCDestroy_KSP(PC pc)
119: {
120: PC_KSP *jac = (PC_KSP*)pc->data;
124: KSPDestroy(&jac->ksp);
125: PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",NULL);
126: PetscObjectComposeFunction((PetscObject)pc,"PCKSPSetKSP_C",NULL);
127: PetscFree(pc->data);
128: return(0);
129: }
131: static PetscErrorCode PCView_KSP(PC pc,PetscViewer viewer)
132: {
133: PC_KSP *jac = (PC_KSP*)pc->data;
135: PetscBool iascii;
138: if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
139: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
140: if (iascii) {
141: if (pc->useAmat) {
142: PetscViewerASCIIPrintf(viewer," Using Amat (not Pmat) as operator on inner solve\n");
143: }
144: PetscViewerASCIIPrintf(viewer," KSP and PC on KSP preconditioner follow\n");
145: PetscViewerASCIIPrintf(viewer," ---------------------------------\n");
146: }
147: PetscViewerASCIIPushTab(viewer);
148: KSPView(jac->ksp,viewer);
149: PetscViewerASCIIPopTab(viewer);
150: if (iascii) {
151: PetscViewerASCIIPrintf(viewer," ---------------------------------\n");
152: }
153: return(0);
154: }
156: static PetscErrorCode PCKSPSetKSP_KSP(PC pc,KSP ksp)
157: {
158: PC_KSP *jac = (PC_KSP*)pc->data;
162: PetscObjectReference((PetscObject)ksp);
163: KSPDestroy(&jac->ksp);
164: jac->ksp = ksp;
165: return(0);
166: }
168: /*@
169: PCKSPSetKSP - Sets the KSP context for a KSP PC.
171: Collective on PC
173: Input Parameter:
174: + pc - the preconditioner context
175: - ksp - the KSP solver
177: Notes:
178: The PC and the KSP must have the same communicator
180: Level: advanced
182: @*/
183: PetscErrorCode PCKSPSetKSP(PC pc,KSP ksp)
184: {
191: PetscTryMethod(pc,"PCKSPSetKSP_C",(PC,KSP),(pc,ksp));
192: return(0);
193: }
195: static PetscErrorCode PCKSPGetKSP_KSP(PC pc,KSP *ksp)
196: {
197: PC_KSP *jac = (PC_KSP*)pc->data;
201: if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
202: *ksp = jac->ksp;
203: return(0);
204: }
206: /*@
207: PCKSPGetKSP - Gets the KSP context for a KSP PC.
209: Not Collective but KSP returned is parallel if PC was parallel
211: Input Parameter:
212: . pc - the preconditioner context
214: Output Parameters:
215: . ksp - the KSP solver
217: Notes:
218: You must call KSPSetUp() before calling PCKSPGetKSP().
220: If the PC is not a PCKSP object it raises an error
222: Level: advanced
224: @*/
225: PetscErrorCode PCKSPGetKSP(PC pc,KSP *ksp)
226: {
232: PetscUseMethod(pc,"PCKSPGetKSP_C",(PC,KSP*),(pc,ksp));
233: return(0);
234: }
236: static PetscErrorCode PCSetFromOptions_KSP(PetscOptionItems *PetscOptionsObject,PC pc)
237: {
238: PC_KSP *jac = (PC_KSP*)pc->data;
242: PetscOptionsHead(PetscOptionsObject,"PC KSP options");
243: if (jac->ksp) {
244: KSPSetFromOptions(jac->ksp);
245: }
246: PetscOptionsTail();
247: return(0);
248: }
250: /* ----------------------------------------------------------------------------------*/
252: /*MC
253: PCKSP - Defines a preconditioner that can consist of any KSP solver.
254: This allows, for example, embedding a Krylov method inside a preconditioner.
256: Options Database Key:
257: . -pc_use_amat - use the matrix that defines the linear system, Amat as the matrix for the
258: inner solver, otherwise by default it uses the matrix used to construct
259: the preconditioner, Pmat (see PCSetOperators())
261: Level: intermediate
263: Notes:
264: Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
265: the incorrect answer) unless you use KSPFGMRES as the other Krylov method
267: Developer Notes:
268: If the outer Krylov method has a nonzero initial guess it will compute a new residual based on that initial guess
269: 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
270: 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
271: input (current residual) vector, the action of the preconditioner cannot depend also on some other vector (the "initial guess"). For
272: KSPRICHARDSON it is possible to provide a PCApplyRichardson_PCKSP() that short circuits returning to the KSP object at each iteration to compute the
273: residual, see for example PCApplyRichardson_SOR(). We do not implement a PCApplyRichardson_PCKSP() because (1) using a KSP directly inside a Richardson
274: is not an efficient algorithm anyways and (2) implementing it for its > 1 would essentially require that we implement Richardson (reimplementing the
275: Richardson code) inside the PCApplyRichardson_PCKSP() leading to duplicate code.
277: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
278: PCSHELL, PCCOMPOSITE, PCSetUseAmat(), PCKSPGetKSP()
280: M*/
282: PETSC_EXTERN PetscErrorCode PCCreate_KSP(PC pc)
283: {
285: PC_KSP *jac;
288: PetscNewLog(pc,&jac);
289: pc->data = (void*)jac;
291: PetscMemzero(pc->ops,sizeof(struct _PCOps));
292: pc->ops->apply = PCApply_KSP;
293: pc->ops->matapply = PCMatApply_KSP;
294: pc->ops->applytranspose = PCApplyTranspose_KSP;
295: pc->ops->setup = PCSetUp_KSP;
296: pc->ops->reset = PCReset_KSP;
297: pc->ops->destroy = PCDestroy_KSP;
298: pc->ops->setfromoptions = PCSetFromOptions_KSP;
299: pc->ops->view = PCView_KSP;
301: PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",PCKSPGetKSP_KSP);
302: PetscObjectComposeFunction((PetscObject)pc,"PCKSPSetKSP_C",PCKSPSetKSP_KSP);
303: return(0);
304: }