Actual source code: pcksp.c
petsc-3.3-p7 2013-05-11
2: #include <petsc-private/pcimpl.h> /*I "petscpc.h" I*/
3: #include <petscksp.h> /*I "petscksp.h" I*/
5: typedef struct {
6: PetscBool use_true_matrix; /* use mat rather than pmat in inner linear solve */
7: KSP ksp;
8: PetscInt its; /* total number of iterations KSP uses */
9: } PC_KSP;
14: static PetscErrorCode PCKSPCreateKSP_KSP(PC pc)
15: {
17: const char *prefix;
18: PC_KSP *jac = (PC_KSP *)pc->data;
21: KSPCreate(((PetscObject)pc)->comm,&jac->ksp);
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: }
31: static PetscErrorCode PCApply_KSP(PC pc,Vec x,Vec y)
32: {
34: PetscInt its;
35: PC_KSP *jac = (PC_KSP*)pc->data;
38: KSPSetInitialGuessNonzero(jac->ksp,pc->nonzero_guess);
39: KSPSolve(jac->ksp,x,y);
40: KSPGetIterationNumber(jac->ksp,&its);
41: jac->its += its;
42: return(0);
43: }
47: static PetscErrorCode PCApplyTranspose_KSP(PC pc,Vec x,Vec y)
48: {
50: PetscInt its;
51: PC_KSP *jac = (PC_KSP*)pc->data;
54: KSPSolveTranspose(jac->ksp,x,y);
55: KSPGetIterationNumber(jac->ksp,&its);
56: jac->its += its;
57: return(0);
58: }
62: static PetscErrorCode PCSetUp_KSP(PC pc)
63: {
65: PC_KSP *jac = (PC_KSP*)pc->data;
66: Mat mat;
67: PetscBool A;
70: if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
71: KSPSetFromOptions(jac->ksp);
72: if (jac->use_true_matrix) mat = pc->mat;
73: else mat = pc->pmat;
75: KSPGetOperatorsSet(jac->ksp,&A,PETSC_NULL);
76: if (!A) {
77: KSPSetOperators(jac->ksp,mat,pc->pmat,pc->flag);
78: } else if (pc->flag != SAME_PRECONDITIONER) {
79: Mat Amat,Bmat;
80: KSPGetOperators(jac->ksp,&Amat,&Bmat,PETSC_NULL);
81: if (Amat == mat && Bmat == pc->pmat) {
82: /* The user has not replaced the matrices so we are expected to forward the update. This incorrectly diagnoses
83: * changed matrices at the top level as the user manually changing the inner matrices, but we have no way to
84: * identify that in this context. The longer term solution is to track matrix state internally.
85: */
86: KSPSetOperators(jac->ksp,mat,pc->pmat,pc->flag);
87: }
88: }
89: KSPSetUp(jac->ksp);
90: return(0);
91: }
93: /* Default destroy, if it has never been setup */
96: static PetscErrorCode PCReset_KSP(PC pc)
97: {
98: PC_KSP *jac = (PC_KSP*)pc->data;
102: if (jac->ksp) {KSPReset(jac->ksp);}
103: return(0);
104: }
108: static PetscErrorCode PCDestroy_KSP(PC pc)
109: {
110: PC_KSP *jac = (PC_KSP*)pc->data;
114: PCReset_KSP(pc);
115: KSPDestroy(&jac->ksp);
116: PetscFree(pc->data);
117: return(0);
118: }
122: static PetscErrorCode PCView_KSP(PC pc,PetscViewer viewer)
123: {
124: PC_KSP *jac = (PC_KSP*)pc->data;
126: PetscBool iascii;
129: if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
130: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
131: if (iascii) {
132: if (jac->use_true_matrix) {
133: PetscViewerASCIIPrintf(viewer,"Using true matrix (not preconditioner matrix) on inner solve\n");
134: }
135: PetscViewerASCIIPrintf(viewer,"KSP and PC on KSP preconditioner follow\n");
136: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
137: } else {
138: SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for this object",((PetscObject)viewer)->type_name);
139: }
140: PetscViewerASCIIPushTab(viewer);
141: KSPView(jac->ksp,viewer);
142: PetscViewerASCIIPopTab(viewer);
143: if (iascii) {
144: PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
145: }
146: return(0);
147: }
151: static PetscErrorCode PCSetFromOptions_KSP(PC pc)
152: {
154: PetscBool flg = PETSC_FALSE;
157: PetscOptionsHead("KSP preconditioner options");
158: PetscOptionsBool("-pc_ksp_true","Use true matrix to define inner linear system, not preconditioner matrix","PCKSPSetUseTrue",flg,&flg,PETSC_NULL);
159: if (flg) {
160: PCKSPSetUseTrue(pc);
161: }
162: PetscOptionsTail();
163: return(0);
164: }
166: /* ----------------------------------------------------------------------------------*/
168: EXTERN_C_BEGIN
171: PetscErrorCode PCKSPSetUseTrue_KSP(PC pc)
172: {
173: PC_KSP *jac;
176: jac = (PC_KSP*)pc->data;
177: jac->use_true_matrix = PETSC_TRUE;
178: return(0);
179: }
180: EXTERN_C_END
182: EXTERN_C_BEGIN
185: PetscErrorCode PCKSPGetKSP_KSP(PC pc,KSP *ksp)
186: {
187: PC_KSP *jac = (PC_KSP*)pc->data;
191: if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
192: *ksp = jac->ksp;
193: return(0);
194: }
195: EXTERN_C_END
199: /*@
200: PCKSPSetUseTrue - Sets a flag to indicate that the true matrix (rather than
201: the matrix used to define the preconditioner) is used to compute the
202: residual inside the inner solve.
204: Logically Collective on PC
206: Input Parameters:
207: . pc - the preconditioner context
209: Options Database Key:
210: . -pc_ksp_true - Activates PCKSPSetUseTrue()
212: Note:
213: For the common case in which the preconditioning and linear
214: system matrices are identical, this routine is unnecessary.
216: Level: advanced
218: .keywords: PC, KSP, set, true, local, flag
220: .seealso: PCSetOperators(), PCBJacobiSetUseTrueLocal()
221: @*/
222: PetscErrorCode PCKSPSetUseTrue(PC pc)
223: {
228: PetscTryMethod(pc,"PCKSPSetUseTrue_C",(PC),(pc));
229: return(0);
230: }
234: /*@
235: PCKSPGetKSP - Gets the KSP context for a KSP PC.
237: Not Collective but KSP returned is parallel if PC was parallel
239: Input Parameter:
240: . pc - the preconditioner context
242: Output Parameters:
243: . ksp - the PC solver
245: Notes:
246: You must call KSPSetUp() before calling PCKSPGetKSP().
248: Level: advanced
250: .keywords: PC, KSP, get, context
251: @*/
252: PetscErrorCode PCKSPGetKSP(PC pc,KSP *ksp)
253: {
259: PetscTryMethod(pc,"PCKSPGetKSP_C",(PC,KSP*),(pc,ksp));
260: return(0);
261: }
263: /* ----------------------------------------------------------------------------------*/
265: /*MC
266: PCKSP - Defines a preconditioner that can consist of any KSP solver.
267: This allows, for example, embedding a Krylov method inside a preconditioner.
269: Options Database Key:
270: . -pc_ksp_true - use the matrix that defines the linear system as the matrix for the
271: inner solver, otherwise by default it uses the matrix used to construct
272: the preconditioner (see PCSetOperators())
274: Level: intermediate
276: Concepts: inner iteration
278: Notes: Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
279: the incorrect answer) unless you use KSPFGMRES as the other Krylov method
281: Developer Notes: PCApply_KSP() uses the flag set by PCSetInitialGuessNonzero(), I think this is totally wrong, because it is then not
282: using this inner KSP as a preconditioner (that is a linear operator applied to some vector), it is actually just using
283: the inner KSP just like the outer KSP.
285: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
286: PCSHELL, PCCOMPOSITE, PCKSPUseTrue(), PCKSPGetKSP()
288: M*/
290: EXTERN_C_BEGIN
293: PetscErrorCode PCCreate_KSP(PC pc)
294: {
296: PC_KSP *jac;
299: PetscNewLog(pc,PC_KSP,&jac);
300: pc->ops->apply = PCApply_KSP;
301: pc->ops->applytranspose = PCApplyTranspose_KSP;
302: pc->ops->setup = PCSetUp_KSP;
303: pc->ops->reset = PCReset_KSP;
304: pc->ops->destroy = PCDestroy_KSP;
305: pc->ops->setfromoptions = PCSetFromOptions_KSP;
306: pc->ops->view = PCView_KSP;
307: pc->ops->applyrichardson = 0;
309: pc->data = (void*)jac;
310:
312: jac->use_true_matrix = PETSC_FALSE;
313: jac->its = 0;
315: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCKSPSetUseTrue_C","PCKSPSetUseTrue_KSP",
316: PCKSPSetUseTrue_KSP);
317: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCKSPGetKSP_C","PCKSPGetKSP_KSP",
318: PCKSPGetKSP_KSP);
320: return(0);
321: }
322: EXTERN_C_END