Actual source code: pcksp.c
petsc-3.4.5 2014-06-29
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: PetscObjectIncrementTabLevel((PetscObject)jac->ksp,(PetscObject)pc,1);
22: PCGetOptionsPrefix(pc,&prefix);
23: KSPSetOptionsPrefix(jac->ksp,prefix);
24: KSPAppendOptionsPrefix(jac->ksp,"ksp_");
25: return(0);
26: }
30: static PetscErrorCode PCApply_KSP(PC pc,Vec x,Vec y)
31: {
33: PetscInt its;
34: PC_KSP *jac = (PC_KSP*)pc->data;
37: KSPSetInitialGuessNonzero(jac->ksp,pc->nonzero_guess);
38: KSPSolve(jac->ksp,x,y);
39: KSPGetIterationNumber(jac->ksp,&its);
40: jac->its += its;
41: return(0);
42: }
46: static PetscErrorCode PCApplyTranspose_KSP(PC pc,Vec x,Vec y)
47: {
49: PetscInt its;
50: PC_KSP *jac = (PC_KSP*)pc->data;
53: KSPSolveTranspose(jac->ksp,x,y);
54: KSPGetIterationNumber(jac->ksp,&its);
55: jac->its += its;
56: return(0);
57: }
61: static PetscErrorCode PCSetUp_KSP(PC pc)
62: {
64: PC_KSP *jac = (PC_KSP*)pc->data;
65: Mat mat;
66: PetscBool A;
69: if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
70: KSPSetFromOptions(jac->ksp);
71: if (pc->useAmat) mat = pc->mat;
72: else mat = pc->pmat;
74: KSPGetOperatorsSet(jac->ksp,&A,NULL);
75: if (!A) {
76: KSPSetOperators(jac->ksp,mat,pc->pmat,pc->flag);
77: } else if (pc->flag != SAME_PRECONDITIONER) {
78: Mat Amat,Bmat;
79: KSPGetOperators(jac->ksp,&Amat,&Bmat,NULL);
80: if (Amat == mat && Bmat == pc->pmat) {
81: /* The user has not replaced the matrices so we are expected to forward the update. This incorrectly diagnoses
82: * changed matrices at the top level as the user manually changing the inner matrices, but we have no way to
83: * identify that in this context. The longer term solution is to track matrix state internally.
84: */
85: KSPSetOperators(jac->ksp,mat,pc->pmat,pc->flag);
86: }
87: }
88: KSPSetUp(jac->ksp);
89: return(0);
90: }
92: /* Default destroy, if it has never been setup */
95: static PetscErrorCode PCReset_KSP(PC pc)
96: {
97: PC_KSP *jac = (PC_KSP*)pc->data;
101: if (jac->ksp) {KSPReset(jac->ksp);}
102: return(0);
103: }
107: static PetscErrorCode PCDestroy_KSP(PC pc)
108: {
109: PC_KSP *jac = (PC_KSP*)pc->data;
113: PCReset_KSP(pc);
114: KSPDestroy(&jac->ksp);
115: PetscFree(pc->data);
116: return(0);
117: }
121: static PetscErrorCode PCView_KSP(PC pc,PetscViewer viewer)
122: {
123: PC_KSP *jac = (PC_KSP*)pc->data;
125: PetscBool iascii;
128: if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
129: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
130: if (iascii) {
131: if (pc->useAmat) {
132: PetscViewerASCIIPrintf(viewer,"Using Amat (not Pmat) as operator on inner solve\n");
133: }
134: PetscViewerASCIIPrintf(viewer,"KSP and PC on KSP preconditioner follow\n");
135: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
136: }
137: PetscViewerASCIIPushTab(viewer);
138: KSPView(jac->ksp,viewer);
139: PetscViewerASCIIPopTab(viewer);
140: if (iascii) {
141: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
142: }
143: return(0);
144: }
148: static PetscErrorCode PCKSPGetKSP_KSP(PC pc,KSP *ksp)
149: {
150: PC_KSP *jac = (PC_KSP*)pc->data;
154: if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
155: *ksp = jac->ksp;
156: return(0);
157: }
161: /*@
162: PCKSPGetKSP - Gets the KSP context for a KSP PC.
164: Not Collective but KSP returned is parallel if PC was parallel
166: Input Parameter:
167: . pc - the preconditioner context
169: Output Parameters:
170: . ksp - the PC solver
172: Notes:
173: You must call KSPSetUp() before calling PCKSPGetKSP().
175: Level: advanced
177: .keywords: PC, KSP, get, context
178: @*/
179: PetscErrorCode PCKSPGetKSP(PC pc,KSP *ksp)
180: {
186: *ksp = NULL;
187: PetscTryMethod(pc,"PCKSPGetKSP_C",(PC,KSP*),(pc,ksp));
188: return(0);
189: }
191: /* ----------------------------------------------------------------------------------*/
193: /*MC
194: PCKSP - Defines a preconditioner that can consist of any KSP solver.
195: This allows, for example, embedding a Krylov method inside a preconditioner.
197: Options Database Key:
198: . -pc_use_amat - use the matrix that defines the linear system, Amat as the matrix for the
199: inner solver, otherwise by default it uses the matrix used to construct
200: the preconditioner, Pmat (see PCSetOperators())
202: Level: intermediate
204: Concepts: inner iteration
206: Notes: Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
207: the incorrect answer) unless you use KSPFGMRES as the other Krylov method
209: Developer Notes: PCApply_KSP() uses the flag set by PCSetInitialGuessNonzero(), I think this is totally wrong, because it is then not
210: using this inner KSP as a preconditioner (that is a linear operator applied to some vector), it is actually just using
211: the inner KSP just like the outer KSP.
213: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
214: PCSHELL, PCCOMPOSITE, PCSetUseAmat(), PCKSPGetKSP()
216: M*/
220: PETSC_EXTERN PetscErrorCode PCCreate_KSP(PC pc)
221: {
223: PC_KSP *jac;
226: PetscNewLog(pc,PC_KSP,&jac);
228: pc->ops->apply = PCApply_KSP;
229: pc->ops->applytranspose = PCApplyTranspose_KSP;
230: pc->ops->setup = PCSetUp_KSP;
231: pc->ops->reset = PCReset_KSP;
232: pc->ops->destroy = PCDestroy_KSP;
233: pc->ops->setfromoptions = 0;
234: pc->ops->view = PCView_KSP;
235: pc->ops->applyrichardson = 0;
237: pc->data = (void*)jac;
240: jac->its = 0;
241: PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",PCKSPGetKSP_KSP);
242: return(0);
243: }