Actual source code: pcksp.c
petsc-3.7.5 2017-01-01
2: #include <petsc/private/pcimpl.h>
3: #include <petscksp.h> /*I "petscksp.h" I*/
5: typedef struct {
6: KSP ksp;
7: PetscInt its; /* total number of iterations KSP uses */
8: } PC_KSP;
13: static PetscErrorCode PCKSPCreateKSP_KSP(PC pc)
14: {
16: const char *prefix;
17: PC_KSP *jac = (PC_KSP*)pc->data;
20: KSPCreate(PetscObjectComm((PetscObject)pc),&jac->ksp);
21: KSPSetErrorIfNotConverged(jac->ksp,pc->erroriffailure);
22: PetscObjectIncrementTabLevel((PetscObject)jac->ksp,(PetscObject)pc,1);
23: PCGetOptionsPrefix(pc,&prefix);
24: KSPSetOptionsPrefix(jac->ksp,prefix);
25: KSPAppendOptionsPrefix(jac->ksp,"ksp_");
26: return(0);
27: }
29: #include <petsc/private/kspimpl.h>
32: static PetscErrorCode PCApply_KSP(PC pc,Vec x,Vec y)
33: {
35: PetscInt its;
36: PC_KSP *jac = (PC_KSP*)pc->data;
39: KSPSolve(jac->ksp,x,y);
40: KSPGetIterationNumber(jac->ksp,&its);
41: jac->its += its;
42: if (jac->ksp->reason == KSP_DIVERGED_PCSETUP_FAILED) {
43: pc->failedreason = PC_SUBPC_ERROR;
44: }
45: return(0);
46: }
50: static PetscErrorCode PCApplyTranspose_KSP(PC pc,Vec x,Vec y)
51: {
53: PetscInt its;
54: PC_KSP *jac = (PC_KSP*)pc->data;
57: KSPSolveTranspose(jac->ksp,x,y);
58: KSPGetIterationNumber(jac->ksp,&its);
59: jac->its += its;
60: return(0);
61: }
65: static PetscErrorCode PCSetUp_KSP(PC pc)
66: {
68: PC_KSP *jac = (PC_KSP*)pc->data;
69: Mat mat;
72: if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
73: KSPSetFromOptions(jac->ksp);
74: if (pc->useAmat) mat = pc->mat;
75: else mat = pc->pmat;
77: KSPSetOperators(jac->ksp,mat,pc->pmat);
78: KSPSetUp(jac->ksp);
79: return(0);
80: }
82: /* Default destroy, if it has never been setup */
85: static PetscErrorCode PCReset_KSP(PC pc)
86: {
87: PC_KSP *jac = (PC_KSP*)pc->data;
91: if (jac->ksp) {KSPReset(jac->ksp);}
92: return(0);
93: }
97: static PetscErrorCode PCDestroy_KSP(PC pc)
98: {
99: PC_KSP *jac = (PC_KSP*)pc->data;
103: PCReset_KSP(pc);
104: KSPDestroy(&jac->ksp);
105: PetscFree(pc->data);
106: return(0);
107: }
111: static PetscErrorCode PCView_KSP(PC pc,PetscViewer viewer)
112: {
113: PC_KSP *jac = (PC_KSP*)pc->data;
115: PetscBool iascii;
118: if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
119: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
120: if (iascii) {
121: if (pc->useAmat) {
122: PetscViewerASCIIPrintf(viewer,"Using Amat (not Pmat) as operator on inner solve\n");
123: }
124: PetscViewerASCIIPrintf(viewer,"KSP and PC on KSP preconditioner follow\n");
125: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
126: }
127: PetscViewerASCIIPushTab(viewer);
128: KSPView(jac->ksp,viewer);
129: PetscViewerASCIIPopTab(viewer);
130: if (iascii) {
131: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
132: }
133: return(0);
134: }
138: static PetscErrorCode PCKSPGetKSP_KSP(PC pc,KSP *ksp)
139: {
140: PC_KSP *jac = (PC_KSP*)pc->data;
144: if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
145: *ksp = jac->ksp;
146: return(0);
147: }
151: /*@
152: PCKSPGetKSP - Gets the KSP context for a KSP PC.
154: Not Collective but KSP returned is parallel if PC was parallel
156: Input Parameter:
157: . pc - the preconditioner context
159: Output Parameters:
160: . ksp - the PC solver
162: Notes:
163: You must call KSPSetUp() before calling PCKSPGetKSP().
165: If the PC is not a PCKSP object then a NULL is returned
167: Level: advanced
169: .keywords: PC, KSP, get, context
170: @*/
171: PetscErrorCode PCKSPGetKSP(PC pc,KSP *ksp)
172: {
178: *ksp = NULL;
179: PetscTryMethod(pc,"PCKSPGetKSP_C",(PC,KSP*),(pc,ksp));
180: return(0);
181: }
183: /* ----------------------------------------------------------------------------------*/
185: /*MC
186: PCKSP - Defines a preconditioner that can consist of any KSP solver.
187: This allows, for example, embedding a Krylov method inside a preconditioner.
189: Options Database Key:
190: . -pc_use_amat - use the matrix that defines the linear system, Amat as the matrix for the
191: inner solver, otherwise by default it uses the matrix used to construct
192: the preconditioner, Pmat (see PCSetOperators())
194: Level: intermediate
196: Concepts: inner iteration
198: Notes: Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
199: the incorrect answer) unless you use KSPFGMRES as the other Krylov method
201: Developer Notes: If the outer Krylov method has a nonzero initial guess it will compute a new residual based on that initial guess
202: 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
203: 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
204: input (current residual) vector, the action of the preconditioner cannot depend also on some other vector (the "initial guess"). For
205: KSPRICHARDSON it is possible to provide a PCApplyRichardson_PCKSP() that short circuits returning to the KSP object at each iteration to compute the
206: residual, see for example PCApplyRichardson_SOR(). We do not implement a PCApplyRichardson_PCKSP() because (1) using a KSP directly inside a Richardson
207: is not an efficient algorithm anyways and (2) implementing it for its > 1 would essentially require that we implement Richardson (reimplementing the
208: Richardson code) inside the PCApplyRichardson_PCKSP() leading to duplicate code.
210: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
211: PCSHELL, PCCOMPOSITE, PCSetUseAmat(), PCKSPGetKSP()
213: M*/
217: PETSC_EXTERN PetscErrorCode PCCreate_KSP(PC pc)
218: {
220: PC_KSP *jac;
223: PetscNewLog(pc,&jac);
225: pc->ops->apply = PCApply_KSP;
226: pc->ops->applytranspose = PCApplyTranspose_KSP;
227: pc->ops->setup = PCSetUp_KSP;
228: pc->ops->reset = PCReset_KSP;
229: pc->ops->destroy = PCDestroy_KSP;
230: pc->ops->setfromoptions = 0;
231: pc->ops->view = PCView_KSP;
232: pc->ops->applyrichardson = 0;
234: pc->data = (void*)jac;
237: jac->its = 0;
238: PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",PCKSPGetKSP_KSP);
239: return(0);
240: }