xref: /petsc/src/ksp/pc/impls/deflation/deflation.c (revision b30d477555b759647cb91caf69a8ec09a58e168d)
1 
2 /*  --------------------------------------------------------------------
3 
4      This file implements a Deflation preconditioner in PETSc as part of PC.
5      You can use this as a starting point for implementing your own
6      preconditioner that is not provided with PETSc. (You might also consider
7      just using PCSHELL)
8 
9      The following basic routines are required for each preconditioner.
10           PCCreate_XXX()          - Creates a preconditioner context
11           PCSetFromOptions_XXX()  - Sets runtime options
12           PCApply_XXX()           - Applies the preconditioner
13           PCDestroy_XXX()         - Destroys the preconditioner context
14      where the suffix "_XXX" denotes a particular implementation, in
15      this case we use _Deflation (e.g., PCCreate_Deflation, PCApply_Deflation).
16      These routines are actually called via the common user interface
17      routines PCCreate(), PCSetFromOptions(), PCApply(), and PCDestroy(),
18      so the application code interface remains identical for all
19      preconditioners.
20 
21      Another key routine is:
22           PCSetUp_XXX()           - Prepares for the use of a preconditioner
23      by setting data structures and options.   The interface routine PCSetUp()
24      is not usually called directly by the user, but instead is called by
25      PCApply() if necessary.
26 
27      Additional basic routines are:
28           PCView_XXX()            - Prints details of runtime options that
29                                     have actually been used.
30      These are called by application codes via the interface routines
31      PCView().
32 
33      The various types of solvers (preconditioners, Krylov subspace methods,
34      nonlinear solvers, timesteppers) are all organized similarly, so the
35      above description applies to these categories also.  One exception is
36      that the analogues of PCApply() for these components are KSPSolve(),
37      SNESSolve(), and TSSolve().
38 
39      Additional optional functionality unique to preconditioners is left and
40      right symmetric preconditioner application via PCApplySymmetricLeft()
41      and PCApplySymmetricRight().  The Deflation implementation is
42      PCApplySymmetricLeftOrRight_Deflation().
43 
44     -------------------------------------------------------------------- */
45 
46 /*
47    Include files needed for the Deflation preconditioner:
48      pcimpl.h - private include file intended for use by all preconditioners
49 */
50 
51 #include <../src/ksp/pc/impls/deflation/deflation.h> /*I "petscpc.h" I*/  /* includes for fortran wrappers */
52 
53 const char *const PCDeflationTypes[]    = {"INIT","PRE","POST","PCDeflationType","PC_DEFLATION_",0};
54 
55 static PetscErrorCode PCDeflationSetSpace_Deflation(PC pc,Mat W,PetscBool transpose)
56 {
57   PC_Deflation   *def = (PC_Deflation*)pc->data;
58   PetscErrorCode ierr;
59 
60   PetscFunctionBegin;
61   if (transpose) {
62     def->Wt = W;
63     def->W = NULL;
64   } else {
65     def->W = W;
66   }
67   ierr = PetscObjectReference((PetscObject)W);CHKERRQ(ierr);
68   PetscFunctionReturn(0);
69 }
70 
71 /* TODO create PCDeflationSetSpaceTranspose? */
72 /*@
73    PCDeflationSetSpace - Set deflation space matrix (or its transpose).
74 
75    Logically Collective on PC
76 
77    Input Parameters:
78 +  pc - the preconditioner context
79 .  W  - deflation matrix
80 -  tranpose - indicates that W is an explicit transpose of the deflation matrix
81 
82    Level: intermediate
83 
84 .seealso: PCDEFLATION
85 @*/
86 PetscErrorCode PCDeflationSetSpace(PC pc,Mat W,PetscBool transpose)
87 {
88   PetscErrorCode ierr;
89 
90   PetscFunctionBegin;
91   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
92   PetscValidHeaderSpecific(W,MAT_CLASSID,2);
93   PetscValidLogicalCollectiveBool(pc,transpose,3);
94   ierr = PetscTryMethod(pc,"PCDeflationSetSpace_C",(PC,Mat,PetscBool),(pc,W,transpose));CHKERRQ(ierr);
95   PetscFunctionReturn(0);
96 }
97 
98 static PetscErrorCode PCDeflationSetLvl_Deflation(PC pc,PetscInt current,PetscInt max)
99 {
100   PC_Deflation   *def = (PC_Deflation*)pc->data;
101 
102   PetscFunctionBegin;
103   def->nestedlvl = current;
104   def->maxnestedlvl = max;
105   PetscFunctionReturn(0);
106 }
107 
108 /*@
109    PCDeflationSetMaxLvl - Set maximum level of deflation.
110 
111    Logically Collective on PC
112 
113    Input Parameters:
114 +  pc  - the preconditioner context
115 .  max - maximum deflation level
116 
117    Level: intermediate
118 
119 .seealso: PCDEFLATION
120 @*/
121 PetscErrorCode PCDeflationSetMaxLvl(PC pc,PetscInt max)
122 {
123   PetscErrorCode ierr;
124 
125   PetscFunctionBegin;
126   PetscValidLogicalCollectiveInt(pc,max,2);
127   ierr = PetscTryMethod(pc,"PCDeflationSetLvl_C",(PC,PetscInt,PetscInt),(pc,0,max));CHKERRQ(ierr);
128   PetscFunctionReturn(0);
129 }
130 
131 /*
132   TODO CP corection?
133   x <- x + W*(W'*A*W)^{-1}*W'*r  = x + Q*r
134 */
135 static PetscErrorCode PCPreSolve_Deflation(PC pc,KSP ksp,Vec b, Vec x)
136 {
137   PC_Deflation     *def = (PC_Deflation*)pc->data;
138   Mat              A;
139   Vec              r,w1,w2;
140   PetscBool        nonzero;
141   PetscErrorCode   ierr;
142 
143   PetscFunctionBegin;
144   w1 = def->workcoarse[0];
145   w2 = def->workcoarse[1];
146   r  = def->work;
147   ierr = PCGetOperators(pc,NULL,&A);CHKERRQ(ierr);
148 
149   ierr = KSPGetInitialGuessNonzero(ksp,&nonzero);CHKERRQ(ierr);
150   ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr);
151   if (nonzero) {
152     ierr = MatMult(A,x,r);CHKERRQ(ierr);               /*    r  <- b - Ax              */
153     ierr = VecAYPX(r,-1.0,b);CHKERRQ(ierr);
154   } else {
155     ierr = VecCopy(b,r);CHKERRQ(ierr);                 /*    r  <- b (x is 0)          */
156   }
157 
158   ierr = MatMultTranspose(def->W,r,w1);CHKERRQ(ierr);  /*    w1 <- W'*r                */
159   ierr = KSPSolve(def->WtAWinv,w1,w2);CHKERRQ(ierr);   /*    w2 <- (W'*A*W)^{-1}*w1    */
160   ierr = MatMult(def->W,w2,r);CHKERRQ(ierr);           /*    r  <- W*w2                */
161   ierr = VecAYPX(x,1.0,r);CHKERRQ(ierr);
162   PetscFunctionReturn(0);
163 }
164 
165 /*
166   if (def->correct) {
167     z <- r - W*(W'*A*W)^{-1}*W'*(A*r -r) = (P-Q)*r
168   } else {
169     z <- r - W*(W'*A*W)^{-1}*W'*A*r = P*r
170   }
171 */
172 static PetscErrorCode PCApply_Deflation(PC pc,Vec r, Vec z)
173 {
174   PC_Deflation     *def = (PC_Deflation*)pc->data;
175   Mat              A;
176   Vec              u,w1,w2;
177   PetscErrorCode   ierr;
178 
179   PetscFunctionBegin;
180   w1 = def->workcoarse[0];
181   w2 = def->workcoarse[1];
182   u  = def->work;
183   ierr = PCGetOperators(pc,NULL,&A);CHKERRQ(ierr);
184 
185   if (!def->AW) {
186     ierr = MatMult(A,r,u);CHKERRQ(ierr);                       /*    u  <- A*r                 */
187     if (def->correct) ierr = VecAXPY(u,-1.0,r);CHKERRQ(ierr);  /*    u  <- A*r -r              */
188     ierr = MatMultTranspose(def->W,u,w1);CHKERRQ(ierr);        /*    w1 <- W'*u                */
189   } else {
190     ierr = MatMultTranspose(def->AW,u,w1);CHKERRQ(ierr);       /*    u  <- A*r                 */
191     if (def->correct) {
192       ierr = MatMultTranspose(def->W,r,w2);CHKERRQ(ierr);      /*    w2 <- W'*u                */
193       ierr = VecAXPY(w1,-1.0,w2);CHKERRQ(ierr);                /*    w1 <- w1 - w2             */
194     }
195   }
196   ierr = KSPSolve(def->WtAWinv,w1,w2);CHKERRQ(ierr);           /*    w2 <- (W'*A*W)^{-1}*w1    */
197   ierr = MatMult(def->W,w2,u);CHKERRQ(ierr);                   /*    u  <- W*w2                */
198   ierr = VecWAXPY(z,-1.0,u,r);CHKERRQ(ierr);                   /*    z  <- r - u               */
199   PetscFunctionReturn(0);
200 }
201 
202 static PetscErrorCode  PCDeflationSetType_Deflation(PC pc,PCDeflationType type)
203 {
204   PC_Deflation *def = (PC_Deflation*)pc->data;
205 
206   PetscFunctionBegin;
207   def->init = PETSC_FALSE;
208   def->pre = PETSC_FALSE;
209   if (type == PC_DEFLATION_POST) {
210     //pc->ops->postsolve = PCPostSolve_Deflation;
211     pc->ops->presolve = 0;
212   } else {
213     pc->ops->presolve = PCPreSolve_Deflation;
214     pc->ops->postsolve = 0;
215     if (type == PC_DEFLATION_INIT) {
216       def->init = PETSC_TRUE;
217       pc->ops->apply = 0;
218     } else {
219       def->pre  = PETSC_TRUE;
220     }
221   }
222   PetscFunctionReturn(0);
223 }
224 
225 /*@
226    PCDeflationSetType - Causes the deflation preconditioner to use only a special
227     initial gues or pre/post solve solution update
228 
229    Logically Collective on PC
230 
231    Input Parameters:
232 +  pc - the preconditioner context
233 -  type - PC_DEFLATION_PRE, PC_DEFLATION_INIT, PC_DEFLATION_POST
234 
235    Options Database Key:
236 .  -pc_deflation_type <pre,init,post>
237 
238    Level: intermediate
239 
240    Concepts: Deflation preconditioner
241 
242 .seealso: PCDeflationGetType()
243 @*/
244 PetscErrorCode  PCDeflationSetType(PC pc,PCDeflationType type)
245 {
246   PetscErrorCode ierr;
247 
248   PetscFunctionBegin;
249   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
250   ierr = PetscTryMethod(pc,"PCDeflationSetType_C",(PC,PCDeflationType),(pc,type));CHKERRQ(ierr);
251   PetscFunctionReturn(0);
252 }
253 
254 static PetscErrorCode  PCDeflationGetType_Deflation(PC pc,PCDeflationType *type)
255 {
256   PC_Deflation *def = (PC_Deflation*)pc->data;
257 
258   PetscFunctionBegin;
259   if (def->init) {
260     *type = PC_DEFLATION_INIT;
261   } else if (def->pre) {
262     *type = PC_DEFLATION_PRE;
263   } else {
264     *type = PC_DEFLATION_POST;
265   }
266   PetscFunctionReturn(0);
267 }
268 
269 /*@
270    PCDeflationGetType - Gets how the diagonal matrix is produced for the preconditioner
271 
272    Not Collective on PC
273 
274    Input Parameter:
275 .  pc - the preconditioner context
276 
277    Output Parameter:
278 -  type - PC_DEFLATION_PRE, PC_DEFLATION_INIT, PC_DEFLATION_POST
279 
280    Level: intermediate
281 
282    Concepts: Deflation preconditioner
283 
284 .seealso: PCDeflationSetType()
285 @*/
286 PetscErrorCode  PCDeflationGetType(PC pc,PCDeflationType *type)
287 {
288   PetscErrorCode ierr;
289 
290   PetscFunctionBegin;
291   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
292   ierr = PetscUseMethod(pc,"PCDeflationGetType_C",(PC,PCDeflationType*),(pc,type));CHKERRQ(ierr);
293   PetscFunctionReturn(0);
294 }
295 
296 static PetscErrorCode PCSetUp_Deflation(PC pc)
297 {
298   PC_Deflation     *def = (PC_Deflation*)pc->data;
299   KSP              innerksp;
300   PC               pcinner;
301   Mat              Amat,nextDef=NULL,*mats;
302   PetscInt         i,m,red,size,commsize;
303   PetscBool        match,flgspd,transp=PETSC_FALSE;
304   MatCompositeType ctype;
305   MPI_Comm         comm;
306   const char       *prefix;
307   PetscErrorCode   ierr;
308 
309   PetscFunctionBegin;
310   if (pc->setupcalled) PetscFunctionReturn(0);
311   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
312   ierr = PCGetOperators(pc,NULL,&Amat);CHKERRQ(ierr);
313 
314   /* compute a deflation space */
315   if (def->W || def->Wt) {
316     def->spacetype = PC_DEFLATION_SPACE_USER;
317   } else {
318     ierr = PCDeflationComputeSpace(pc);CHKERRQ(ierr);
319   }
320 
321   /* nested deflation */
322   if (def->W) {
323     ierr = PetscObjectTypeCompare((PetscObject)def->W,MATCOMPOSITE,&match);CHKERRQ(ierr);
324     if (match) {
325       ierr = MatCompositeGetType(def->W,&ctype);CHKERRQ(ierr);
326       ierr = MatCompositeGetNumberMat(def->W,&size);CHKERRQ(ierr);
327     }
328   } else {
329     ierr = MatCreateTranspose(def->Wt,&def->W);CHKERRQ(ierr);
330     ierr = PetscObjectTypeCompare((PetscObject)def->Wt,MATCOMPOSITE,&match);CHKERRQ(ierr);
331     if (match) {
332       ierr = MatCompositeGetType(def->Wt,&ctype);CHKERRQ(ierr);
333       ierr = MatCompositeGetNumberMat(def->Wt,&size);CHKERRQ(ierr);
334     }
335     transp = PETSC_TRUE;
336   }
337   if (match && ctype == MAT_COMPOSITE_MULTIPLICATIVE) {
338     ierr = PetscMalloc1(size,&mats);CHKERRQ(ierr);
339     if (!transp) {
340       for (i=0; i<size; i++) {
341         ierr = MatCompositeGetMat(def->W,i,&mats[i]);CHKERRQ(ierr);
342         //ierr = PetscObjectReference((PetscObject)mats[i]);CHKERRQ(ierr);
343       }
344       if (def->nestedlvl < def->maxnestedlvl) {
345         size -= 1;
346         ierr = MatDestroy(&def->W);CHKERRQ(ierr);
347         def->W = mats[size];
348         ierr = PetscObjectReference((PetscObject)mats[size]);CHKERRQ(ierr);
349         if (size > 1) {
350           ierr = MatCreateComposite(comm,size,mats,&nextDef);CHKERRQ(ierr);
351           ierr = MatCompositeSetType(nextDef,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr);
352         } else {
353           nextDef = mats[0];
354           ierr = PetscObjectReference((PetscObject)mats[0]);CHKERRQ(ierr);
355         }
356       } else {
357         /* ierr = MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT);CHKERRQ(ierr); */
358         ierr = MatCompositeMerge(def->W);CHKERRQ(ierr);
359       }
360     }
361     ierr = PetscFree(mats);CHKERRQ(ierr);
362   }
363 
364   /* setup coarse problem */
365   if (!def->WtAWinv) {
366     ierr = MatGetSize(def->W,NULL,&m);CHKERRQ(ierr); /* TODO works for W MatTranspose? */
367     if (!def->WtAW) {
368       /* TODO add implicit product version ? */
369       if (!def->AW) {
370         ierr = MatPtAP(Amat,def->W,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtAW);CHKERRQ(ierr);
371       } else {
372         ierr = MatTransposeMatMult(def->W,def->AW,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtAW);CHKERRQ(ierr);
373       }
374       /* TODO create MatInheritOption(Mat,MatOption) */
375       ierr = MatGetOption(Amat,MAT_SPD,&flgspd);CHKERRQ(ierr);
376       ierr = MatSetOption(def->WtAW,MAT_SPD,flgspd);CHKERRQ(ierr);
377 #if defined(PETSC_USE_DEBUG)
378       /* Check WtAW is not sigular */
379       PetscReal *norms;
380       ierr = PetscMalloc1(m,&norms);CHKERRQ(ierr);
381       ierr = MatGetColumnNorms(def->WtAW,NORM_INFINITY,norms);CHKERRQ(ierr);
382       for (i=0; i<m; i++) {
383         if (norms[i] < 100*PETSC_MACHINE_EPSILON) {
384           SETERRQ1(comm,PETSC_ERR_SUP,"Column %D of W is in kernel of A.",i);
385         }
386       }
387       ierr = PetscFree(norms);CHKERRQ(ierr);
388 #endif
389     } else {
390       ierr = MatGetOption(def->WtAW,MAT_SPD,&flgspd);CHKERRQ(ierr);
391     }
392     /* TODO use MATINV */
393     ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr);
394     ierr = KSPCreate(comm,&def->WtAWinv);CHKERRQ(ierr);
395     ierr = KSPSetOperators(def->WtAWinv,def->WtAW,def->WtAW);CHKERRQ(ierr);
396     ierr = KSPSetType(def->WtAWinv,KSPPREONLY);CHKERRQ(ierr);
397     ierr = KSPGetPC(def->WtAWinv,&pcinner);CHKERRQ(ierr);
398     ierr = PCSetType(pcinner,PCTELESCOPE);CHKERRQ(ierr);
399     /* ugly hack to not have overwritten PCTELESCOPE */
400     if (prefix) {
401       ierr = KSPSetOptionsPrefix(def->WtAWinv,prefix);CHKERRQ(ierr);
402       ierr = KSPAppendOptionsPrefix(def->WtAWinv,"tel_");CHKERRQ(ierr);
403     } else {
404       ierr = KSPSetOptionsPrefix(def->WtAWinv,"tel_");CHKERRQ(ierr);
405     }
406     ierr = PCSetFromOptions(pcinner);CHKERRQ(ierr);
407     /* Reduction factor choice */
408     red = def->reductionfact;
409     if (red < 0) {
410       ierr = MPI_Comm_size(comm,&commsize);CHKERRQ(ierr);
411       red  = ceil((float)commsize/ceil((float)m/commsize));
412       ierr = PetscObjectTypeCompareAny((PetscObject)(def->WtAW),&match,MATSEQDENSE,MATMPIDENSE,MATDENSE,"");CHKERRQ(ierr);
413       if (match) red = commsize;
414       ierr = PetscInfo1(pc,"Auto choosing reduction factor %D\n",red);CHKERRQ(ierr); /* TODO add level? */
415     }
416     ierr = PCTelescopeSetReductionFactor(pcinner,red);CHKERRQ(ierr);
417     ierr = PCSetUp(pcinner);CHKERRQ(ierr);
418     ierr = PCTelescopeGetKSP(pcinner,&innerksp);CHKERRQ(ierr);
419     if (innerksp) {
420       ierr = KSPGetPC(innerksp,&pcinner);CHKERRQ(ierr);
421       /* Setup KSP and PC */
422       if (nextDef) { /* next level for multilevel deflation */
423         /* set default KSPtype */
424         if (!def->ksptype) {
425           def->ksptype = KSPFGMRES;
426           if (flgspd) { /* SPD system */
427             def->ksptype = KSPFCG;
428           }
429         }
430         ierr = KSPSetType(innerksp,def->ksptype);CHKERRQ(ierr); /* TODO iherit from KSP */
431         ierr = PCSetType(pcinner,PCDEFLATION);CHKERRQ(ierr); /* TODO create coarse preconditinoner M_c = WtMW ? */
432         ierr = PCDeflationSetSpace(pcinner,nextDef,transp);CHKERRQ(ierr);
433         ierr = PCDeflationSetLvl_Deflation(pcinner,def->nestedlvl+1,def->maxnestedlvl);CHKERRQ(ierr);
434         /* inherit options TODO if not set */
435         ((PC_Deflation*)(pcinner))->ksptype = def->ksptype;
436         ((PC_Deflation*)(pcinner))->correct = def->correct;
437         ((PC_Deflation*)(pcinner))->adaptiveconv = def->adaptiveconv;
438         ((PC_Deflation*)(pcinner))->adaptiveconst = def->adaptiveconst;
439         ierr = MatDestroy(&nextDef);CHKERRQ(ierr);
440       } else { /* the last level */
441         ierr = KSPSetType(innerksp,KSPGMRES);CHKERRQ(ierr);
442         /* TODO Cholesky if flgspd? */
443         ierr = PCSetType(pcinner,PCLU);CHKERRQ(ierr);
444         //TODO remove explicit matSolverPackage
445         if (commsize == red) {
446           ierr = PCFactorSetMatSolverType(pcinner,MATSOLVERSUPERLU);CHKERRQ(ierr);
447         } else {
448           ierr = PCFactorSetMatSolverType(pcinner,MATSOLVERSUPERLU_DIST);CHKERRQ(ierr);
449         }
450       }
451       /* TODO use def_[lvl]_ if lvl > 0? */
452       if (prefix) {
453         ierr = KSPSetOptionsPrefix(innerksp,prefix);CHKERRQ(ierr);
454         ierr = KSPAppendOptionsPrefix(innerksp,"def_");CHKERRQ(ierr);
455       } else {
456         ierr = KSPSetOptionsPrefix(innerksp,"def_");CHKERRQ(ierr);
457       }
458     }
459     /* TODO: check if WtAWinv is KSP and move following from this if */
460     ierr = KSPSetFromOptions(def->WtAWinv);CHKERRQ(ierr);
461     //if (def->adaptiveconv) {
462     //  PetscReal *rnorm;
463     //  PetscNew(&rnorm);
464     //  ierr = KSPSetConvergenceTest(def->WtAWinv,KSPDCGConvergedAdaptive_DCG,rnorm,NULL);CHKERRQ(ierr);
465     //}
466     ierr = KSPSetUp(def->WtAWinv);CHKERRQ(ierr);
467   }
468 
469   /* create work vecs */
470   ierr = MatCreateVecs(Amat,NULL,&def->work);CHKERRQ(ierr);
471   ierr = KSPCreateVecs(def->WtAWinv,2,&def->workcoarse,0,NULL);CHKERRQ(ierr);
472   PetscFunctionReturn(0);
473 }
474 
475 static PetscErrorCode PCReset_Deflation(PC pc)
476 {
477   PC_Deflation      *def = (PC_Deflation*)pc->data;
478   PetscErrorCode    ierr;
479 
480   PetscFunctionBegin;
481   ierr = VecDestroy(&def->work);CHKERRQ(ierr);
482   ierr = VecDestroyVecs(2,&def->workcoarse);CHKERRQ(ierr);
483   ierr = MatDestroy(&def->W);CHKERRQ(ierr);
484   ierr = MatDestroy(&def->Wt);CHKERRQ(ierr);
485   ierr = MatDestroy(&def->AW);CHKERRQ(ierr);
486   ierr = MatDestroy(&def->WtAW);CHKERRQ(ierr);
487   ierr = KSPDestroy(&def->WtAWinv);CHKERRQ(ierr);
488   PetscFunctionReturn(0);
489 }
490 
491 /*
492    PCDestroy_Deflation - Destroys the private context for the Deflation preconditioner
493    that was created with PCCreate_Deflation().
494 
495    Input Parameter:
496 .  pc - the preconditioner context
497 
498    Application Interface Routine: PCDestroy()
499 */
500 static PetscErrorCode PCDestroy_Deflation(PC pc)
501 {
502   PetscErrorCode ierr;
503 
504   PetscFunctionBegin;
505   ierr = PCReset_Deflation(pc);CHKERRQ(ierr);
506   ierr = PetscFree(pc->data);CHKERRQ(ierr);
507   PetscFunctionReturn(0);
508 }
509 
510 static PetscErrorCode PCSetFromOptions_Deflation(PetscOptionItems *PetscOptionsObject,PC pc)
511 {
512   PC_Deflation      *jac = (PC_Deflation*)pc->data;
513   PetscErrorCode ierr;
514   PetscBool      flg;
515   PCDeflationType   deflt,type;
516 
517   PetscFunctionBegin;
518   ierr = PCDeflationGetType(pc,&deflt);CHKERRQ(ierr);
519   ierr = PetscOptionsHead(PetscOptionsObject,"Deflation options");CHKERRQ(ierr);
520   ierr = PetscOptionsEnum("-pc_jacobi_type","How to construct diagonal matrix","PCDeflationSetType",PCDeflationTypes,(PetscEnum)deflt,(PetscEnum*)&type,&flg);CHKERRQ(ierr);
521   if (flg) {
522     ierr = PCDeflationSetType(pc,type);CHKERRQ(ierr);
523   }
524   ierr = PetscOptionsTail();CHKERRQ(ierr);
525   PetscFunctionReturn(0);
526 }
527 
528 /*MC
529      PCDEFLATION - Deflation preconditioner shifts part of the spectrum to zero (deflates)
530      or to a predefined value
531 
532    Options Database Key:
533 +    -pc_deflation_type <init,pre,post> - selects approach to deflation (default: pre)
534 -    -pc_jacobi_abs - use the absolute value of the diagonal entry
535 
536    Level: beginner
537 
538   Notes:
539     todo
540 
541 .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
542            PCDeflationSetType(), PCDeflationSetSpace()
543 M*/
544 
545 PETSC_EXTERN PetscErrorCode PCCreate_Deflation(PC pc)
546 {
547   PC_Deflation   *def;
548   PetscErrorCode ierr;
549 
550   PetscFunctionBegin;
551   ierr     = PetscNewLog(pc,&def);CHKERRQ(ierr);
552   pc->data = (void*)def;
553 
554   def->init          = PETSC_FALSE;
555   def->pre           = PETSC_TRUE;
556   def->correct       = PETSC_FALSE;
557   def->truenorm      = PETSC_TRUE;
558   def->reductionfact = -1;
559   def->spacetype     = PC_DEFLATION_SPACE_HAAR;
560   def->spacesize     = 1;
561   def->extendsp      = PETSC_FALSE;
562   def->nestedlvl     = 0;
563   def->maxnestedlvl  = 0;
564   def->adaptiveconv  = PETSC_FALSE;
565   def->adaptiveconst = 1.0;
566 
567   /*
568       Set the pointers for the functions that are provided above.
569       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
570       are called, they will automatically call these functions.  Note we
571       choose not to provide a couple of these functions since they are
572       not needed.
573   */
574   pc->ops->apply               = PCApply_Deflation;
575   pc->ops->applytranspose      = PCApply_Deflation;
576   pc->ops->presolve            = PCPreSolve_Deflation;
577   pc->ops->postsolve           = 0;
578   pc->ops->setup               = PCSetUp_Deflation;
579   pc->ops->reset               = PCReset_Deflation;
580   pc->ops->destroy             = PCDestroy_Deflation;
581   pc->ops->setfromoptions      = PCSetFromOptions_Deflation;
582   pc->ops->view                = 0;
583   pc->ops->applyrichardson     = 0;
584   pc->ops->applysymmetricleft  = 0;
585   pc->ops->applysymmetricright = 0;
586 
587   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetType_C",PCDeflationSetType_Deflation);CHKERRQ(ierr);
588   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationGetType_C",PCDeflationGetType_Deflation);CHKERRQ(ierr);
589   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetSpace_C",PCDeflationSetSpace_Deflation);CHKERRQ(ierr);
590   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetLvl_C",PCDeflationSetLvl_Deflation);CHKERRQ(ierr);
591   PetscFunctionReturn(0);
592 }
593 
594