xref: /petsc/src/ksp/pc/impls/deflation/deflation.c (revision b27c8092729a0e9494a7f32c95d604bc11dc3a17)
137eeb815SJakub Kruzik 
237eeb815SJakub Kruzik /*  --------------------------------------------------------------------
337eeb815SJakub Kruzik 
437eeb815SJakub Kruzik      This file implements a Deflation preconditioner in PETSc as part of PC.
537eeb815SJakub Kruzik      You can use this as a starting point for implementing your own
637eeb815SJakub Kruzik      preconditioner that is not provided with PETSc. (You might also consider
737eeb815SJakub Kruzik      just using PCSHELL)
837eeb815SJakub Kruzik 
937eeb815SJakub Kruzik      The following basic routines are required for each preconditioner.
1037eeb815SJakub Kruzik           PCCreate_XXX()          - Creates a preconditioner context
1137eeb815SJakub Kruzik           PCSetFromOptions_XXX()  - Sets runtime options
1237eeb815SJakub Kruzik           PCApply_XXX()           - Applies the preconditioner
1337eeb815SJakub Kruzik           PCDestroy_XXX()         - Destroys the preconditioner context
1437eeb815SJakub Kruzik      where the suffix "_XXX" denotes a particular implementation, in
1537eeb815SJakub Kruzik      this case we use _Deflation (e.g., PCCreate_Deflation, PCApply_Deflation).
1637eeb815SJakub Kruzik      These routines are actually called via the common user interface
1737eeb815SJakub Kruzik      routines PCCreate(), PCSetFromOptions(), PCApply(), and PCDestroy(),
1837eeb815SJakub Kruzik      so the application code interface remains identical for all
1937eeb815SJakub Kruzik      preconditioners.
2037eeb815SJakub Kruzik 
2137eeb815SJakub Kruzik      Another key routine is:
2237eeb815SJakub Kruzik           PCSetUp_XXX()           - Prepares for the use of a preconditioner
2337eeb815SJakub Kruzik      by setting data structures and options.   The interface routine PCSetUp()
2437eeb815SJakub Kruzik      is not usually called directly by the user, but instead is called by
2537eeb815SJakub Kruzik      PCApply() if necessary.
2637eeb815SJakub Kruzik 
2737eeb815SJakub Kruzik      Additional basic routines are:
2837eeb815SJakub Kruzik           PCView_XXX()            - Prints details of runtime options that
2937eeb815SJakub Kruzik                                     have actually been used.
3037eeb815SJakub Kruzik      These are called by application codes via the interface routines
3137eeb815SJakub Kruzik      PCView().
3237eeb815SJakub Kruzik 
3337eeb815SJakub Kruzik      The various types of solvers (preconditioners, Krylov subspace methods,
3437eeb815SJakub Kruzik      nonlinear solvers, timesteppers) are all organized similarly, so the
3537eeb815SJakub Kruzik      above description applies to these categories also.  One exception is
3637eeb815SJakub Kruzik      that the analogues of PCApply() for these components are KSPSolve(),
3737eeb815SJakub Kruzik      SNESSolve(), and TSSolve().
3837eeb815SJakub Kruzik 
3937eeb815SJakub Kruzik      Additional optional functionality unique to preconditioners is left and
4037eeb815SJakub Kruzik      right symmetric preconditioner application via PCApplySymmetricLeft()
4137eeb815SJakub Kruzik      and PCApplySymmetricRight().  The Deflation implementation is
4237eeb815SJakub Kruzik      PCApplySymmetricLeftOrRight_Deflation().
4337eeb815SJakub Kruzik 
4437eeb815SJakub Kruzik     -------------------------------------------------------------------- */
4537eeb815SJakub Kruzik 
4637eeb815SJakub Kruzik /*
4737eeb815SJakub Kruzik    Include files needed for the Deflation preconditioner:
4837eeb815SJakub Kruzik      pcimpl.h - private include file intended for use by all preconditioners
4937eeb815SJakub Kruzik */
5037eeb815SJakub Kruzik 
5137eeb815SJakub Kruzik #include <petsc/private/pcimpl.h>   /*I "petscpc.h" I*/
5237eeb815SJakub Kruzik 
5337eeb815SJakub Kruzik const char *const PCDeflationTypes[]    = {"INIT","PRE","POST","PCDeflationType","PC_DEFLATION_",0};
5437eeb815SJakub Kruzik 
5537eeb815SJakub Kruzik /*
5637eeb815SJakub Kruzik    Private context (data structure) for the deflation preconditioner.
5737eeb815SJakub Kruzik */
5837eeb815SJakub Kruzik typedef struct {
5937eeb815SJakub Kruzik   PetscBool init;            /* do only init step - error correction of direction is omitted */
6037eeb815SJakub Kruzik   PetscBool pre;             /* start with x0 being the solution in the deflation space */
6137eeb815SJakub Kruzik   PetscBool correct;         /* add CP (Qr) correction to descent direction */
6237eeb815SJakub Kruzik   PetscBool truenorm;
6337eeb815SJakub Kruzik   PetscBool adaptiveconv;
6437eeb815SJakub Kruzik   PetscReal adaptiveconst;
6537eeb815SJakub Kruzik   PetscInt  reductionfact;
6637eeb815SJakub Kruzik   Mat       W,Wt,AW,WtAW;    /* deflation space, coarse problem mats */
6737eeb815SJakub Kruzik   KSP       WtAWinv;         /* deflation coarse problem */
68409a357bSJakub Kruzik   KSPType   ksptype;
69f8bf229cSJakub Kruzik   Vec       work;
70f8bf229cSJakub Kruzik   Vec       *workcoarse;
7137eeb815SJakub Kruzik 
72409a357bSJakub Kruzik   PCDeflationSpaceType spacetype;
7337eeb815SJakub Kruzik   PetscInt             spacesize;
7437eeb815SJakub Kruzik   PetscInt             nestedlvl;
7537eeb815SJakub Kruzik   PetscInt             maxnestedlvl;
7637eeb815SJakub Kruzik   PetscBool            extendsp;
7737eeb815SJakub Kruzik } PC_Deflation;
7837eeb815SJakub Kruzik 
79e662bc50SJakub Kruzik static PetscErrorCode PCDeflationSetSpace_Deflation(PC pc,Mat W,PetscBool transpose)
80e662bc50SJakub Kruzik {
81e662bc50SJakub Kruzik   PC_Deflation   *def = (PC_Deflation*)pc->data;
82e662bc50SJakub Kruzik   PetscErrorCode ierr;
83e662bc50SJakub Kruzik 
84e662bc50SJakub Kruzik   PetscFunctionBegin;
85e662bc50SJakub Kruzik   if (transpose) {
86e662bc50SJakub Kruzik     def->Wt = W;
87e662bc50SJakub Kruzik     def->W = NULL;
88e662bc50SJakub Kruzik   } else {
89e662bc50SJakub Kruzik     def->W = W;
90e662bc50SJakub Kruzik   }
91e662bc50SJakub Kruzik   ierr = PetscObjectReference((PetscObject)W);CHKERRQ(ierr);
92e662bc50SJakub Kruzik   PetscFunctionReturn(0);
93e662bc50SJakub Kruzik }
94e662bc50SJakub Kruzik 
95e662bc50SJakub Kruzik /* TODO create PCDeflationSetSpaceTranspose? */
96e662bc50SJakub Kruzik /*@
97e662bc50SJakub Kruzik    PCDeflationSetSpace - Set deflation space matrix (or its transpose).
98e662bc50SJakub Kruzik 
99e662bc50SJakub Kruzik    Logically Collective on PC
100e662bc50SJakub Kruzik 
101e662bc50SJakub Kruzik    Input Parameters:
102e662bc50SJakub Kruzik +  pc - the preconditioner context
103e662bc50SJakub Kruzik .  W  - deflation matrix
104e662bc50SJakub Kruzik -  tranpose - indicates that W is an explicit transpose of the deflation matrix
105e662bc50SJakub Kruzik 
106e662bc50SJakub Kruzik    Level: intermediate
107e662bc50SJakub Kruzik 
108e662bc50SJakub Kruzik .seealso: PCDEFLATION
109e662bc50SJakub Kruzik @*/
110e662bc50SJakub Kruzik PetscErrorCode PCDeflationSetSpace(PC pc,Mat W,PetscBool transpose)
111e662bc50SJakub Kruzik {
112e662bc50SJakub Kruzik   PetscErrorCode ierr;
113e662bc50SJakub Kruzik 
114e662bc50SJakub Kruzik   PetscFunctionBegin;
115e662bc50SJakub Kruzik   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
116e662bc50SJakub Kruzik   PetscValidHeaderSpecific(W,MAT_CLASSID,2);
117e662bc50SJakub Kruzik   PetscValidLogicalCollectiveBool(pc,transpose,3);
118e662bc50SJakub Kruzik   ierr = PetscTryMethod(pc,"PCDeflationSetSpace_C",(PC,Mat,PetscBool),(pc,W,transpose));CHKERRQ(ierr);
119e662bc50SJakub Kruzik   PetscFunctionReturn(0);
120e662bc50SJakub Kruzik }
121e662bc50SJakub Kruzik 
1225c0c31c5SJakub Kruzik static PetscErrorCode PCDeflationSetLvl_Deflation(PC pc,PetscInt current,PetscInt max)
1235c0c31c5SJakub Kruzik {
1245c0c31c5SJakub Kruzik   PC_Deflation   *def = (PC_Deflation*)pc->data;
1255c0c31c5SJakub Kruzik 
1265c0c31c5SJakub Kruzik   PetscFunctionBegin;
1275c0c31c5SJakub Kruzik   def->nestedlvl = current;
1285c0c31c5SJakub Kruzik   def->maxnestedlvl = max;
1295c0c31c5SJakub Kruzik   PetscFunctionReturn(0);
1305c0c31c5SJakub Kruzik }
1315c0c31c5SJakub Kruzik 
1325c0c31c5SJakub Kruzik /*@
1335c0c31c5SJakub Kruzik    PCDeflationSetMaxLvl - Set maximum level of deflation.
1345c0c31c5SJakub Kruzik 
1355c0c31c5SJakub Kruzik    Logically Collective on PC
1365c0c31c5SJakub Kruzik 
1375c0c31c5SJakub Kruzik    Input Parameters:
1385c0c31c5SJakub Kruzik +  pc  - the preconditioner context
1395c0c31c5SJakub Kruzik .  max - maximum deflation level
1405c0c31c5SJakub Kruzik 
1415c0c31c5SJakub Kruzik    Level: intermediate
1425c0c31c5SJakub Kruzik 
1435c0c31c5SJakub Kruzik .seealso: PCDEFLATION
1445c0c31c5SJakub Kruzik @*/
1455c0c31c5SJakub Kruzik PetscErrorCode PCDeflationSetMaxLvl(PC pc,PetscInt max)
1465c0c31c5SJakub Kruzik {
1475c0c31c5SJakub Kruzik   PetscErrorCode ierr;
1485c0c31c5SJakub Kruzik 
1495c0c31c5SJakub Kruzik   PetscFunctionBegin;
1505c0c31c5SJakub Kruzik   PetscValidLogicalCollectiveInt(pc,max,2);
1515c0c31c5SJakub Kruzik   ierr = PetscTryMethod(pc,"PCDeflationSetLvl_C",(PC,PetscInt,PetscInt),(pc,0,max));CHKERRQ(ierr);
1525c0c31c5SJakub Kruzik   PetscFunctionReturn(0);
1535c0c31c5SJakub Kruzik }
154e662bc50SJakub Kruzik 
15537eeb815SJakub Kruzik /*
156*b27c8092SJakub Kruzik   TODO CP corection?
157*b27c8092SJakub Kruzik   x <- x + W*(W'*A*W)^{-1}*W'*r  = x + Q*r
158*b27c8092SJakub Kruzik */
159*b27c8092SJakub Kruzik static PetscErrorCode PCPreSolve_Deflation(PC pc,KSP ksp,Vec b, Vec x)
160*b27c8092SJakub Kruzik {
161*b27c8092SJakub Kruzik   PC_Deflation     *def = (PC_Deflation*)pc->data;
162*b27c8092SJakub Kruzik   Mat              A;
163*b27c8092SJakub Kruzik   Vec              r,w1,w2;
164*b27c8092SJakub Kruzik   PetscBool        nonzero;
165*b27c8092SJakub Kruzik   PetscErrorCode   ierr;
16637eeb815SJakub Kruzik 
167*b27c8092SJakub Kruzik   PetscFunctionBegin;
168*b27c8092SJakub Kruzik   w1 = def->workcoarse[0];
169*b27c8092SJakub Kruzik   w2 = def->workcoarse[1];
170*b27c8092SJakub Kruzik   r  = def->work;
171*b27c8092SJakub Kruzik   ierr = PCGetOperators(pc,NULL,&A);CHKERRQ(ierr);
17237eeb815SJakub Kruzik 
173*b27c8092SJakub Kruzik   ierr = KSPGetInitialGuessNonzero(ksp,&nonzero);CHKERRQ(ierr);
174*b27c8092SJakub Kruzik   ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr);
175*b27c8092SJakub Kruzik   if (nonzero) {
176*b27c8092SJakub Kruzik     ierr = MatMult(A,x,r);CHKERRQ(ierr);               /*    r  <- b - Ax              */
177*b27c8092SJakub Kruzik     ierr = VecAYPX(r,-1.0,b);CHKERRQ(ierr);
178*b27c8092SJakub Kruzik   } else {
179*b27c8092SJakub Kruzik     ierr = VecCopy(b,r);CHKERRQ(ierr);                 /*    r  <- b (x is 0)          */
180*b27c8092SJakub Kruzik   }
181*b27c8092SJakub Kruzik 
182*b27c8092SJakub Kruzik   ierr = MatMultTranspose(def->W,r,w1);CHKERRQ(ierr);  /*    w1 <- W'*r                */
183*b27c8092SJakub Kruzik   ierr = KSPSolve(def->WtAWinv,w1,w2);CHKERRQ(ierr);   /*    w2 <- (W'*A*W)^{-1}*w1    */
184*b27c8092SJakub Kruzik   ierr = MatMult(def->W,w2,r);CHKERRQ(ierr);           /*    r  <- W*w2                */
185*b27c8092SJakub Kruzik   ierr = VecAYPX(x,1.0,r);CHKERRQ(ierr);
186*b27c8092SJakub Kruzik   PetscFunctionReturn(0);
187*b27c8092SJakub Kruzik }
18837eeb815SJakub Kruzik 
189f8bf229cSJakub Kruzik /*
190f8bf229cSJakub Kruzik   if (def->correct) {
191f8bf229cSJakub Kruzik     z <- r - W*(W'*A*W)^{-1}*W'*(A*r -r) = (P-Q)*r
192f8bf229cSJakub Kruzik   } else {
193f8bf229cSJakub Kruzik     z <- r - W*(W'*A*W)^{-1}*W'*A*r = P*r
194f8bf229cSJakub Kruzik   }
19537eeb815SJakub Kruzik */
196f8bf229cSJakub Kruzik static PetscErrorCode PCApply_Deflation(PC pc,Vec r, Vec z)
197f8bf229cSJakub Kruzik {
198f8bf229cSJakub Kruzik   PC_Deflation     *def = (PC_Deflation*)pc->data;
199f8bf229cSJakub Kruzik   Mat              A;
200f8bf229cSJakub Kruzik   Vec              u,w1,w2;
201f8bf229cSJakub Kruzik   PetscErrorCode   ierr;
202f8bf229cSJakub Kruzik 
203f8bf229cSJakub Kruzik   PetscFunctionBegin;
204f8bf229cSJakub Kruzik   w1 = def->workcoarse[0];
205f8bf229cSJakub Kruzik   w2 = def->workcoarse[1];
206f8bf229cSJakub Kruzik   u  = def->work;
207f8bf229cSJakub Kruzik   ierr = PCGetOperators(pc,NULL,&A);CHKERRQ(ierr);
208f8bf229cSJakub Kruzik 
209f8bf229cSJakub Kruzik   if (!def->AW) {
210f8bf229cSJakub Kruzik     ierr = MatMult(A,r,u);CHKERRQ(ierr);                       /*    u  <- A*r                 */
211f8bf229cSJakub Kruzik     if (def->correct) ierr = VecAXPY(u,-1.0,r);CHKERRQ(ierr);  /*    u  <- A*r -r              */
212f8bf229cSJakub Kruzik     ierr = MatMultTranspose(def->W,u,w1);CHKERRQ(ierr);        /*    w1 <- W'*u                */
213f8bf229cSJakub Kruzik   } else {
214f8bf229cSJakub Kruzik     ierr = MatMultTranspose(def->AW,u,w1);CHKERRQ(ierr);       /*    u  <- A*r                 */
215f8bf229cSJakub Kruzik     if (def->correct) {
216f8bf229cSJakub Kruzik       ierr = MatMultTranspose(def->W,r,w2);CHKERRQ(ierr);      /*    w2 <- W'*u                */
217f8bf229cSJakub Kruzik       ierr = VecAXPY(w1,-1.0,w2);CHKERRQ(ierr);                /*    w1 <- w1 - w2             */
218f8bf229cSJakub Kruzik     }
219f8bf229cSJakub Kruzik   }
220f8bf229cSJakub Kruzik   ierr = KSPSolve(def->WtAWinv,w1,w2);CHKERRQ(ierr);           /*    w2 <- (W'*A*W)^{-1}*w1    */
221f8bf229cSJakub Kruzik   ierr = MatMult(def->W,w2,u);CHKERRQ(ierr);                   /*    u  <- W*w2                */
222f8bf229cSJakub Kruzik   ierr = VecWAXPY(z,-1.0,u,r);CHKERRQ(ierr);                   /*    z  <- r - u               */
223f8bf229cSJakub Kruzik   PetscFunctionReturn(0);
224f8bf229cSJakub Kruzik }
225f8bf229cSJakub Kruzik 
226*b27c8092SJakub Kruzik static PetscErrorCode  PCDeflationSetType_Deflation(PC pc,PCDeflationType type)
227*b27c8092SJakub Kruzik {
228*b27c8092SJakub Kruzik   PC_Deflation *def = (PC_Deflation*)pc->data;
229*b27c8092SJakub Kruzik 
230*b27c8092SJakub Kruzik   PetscFunctionBegin;
231*b27c8092SJakub Kruzik   def->init = PETSC_FALSE;
232*b27c8092SJakub Kruzik   def->pre = PETSC_FALSE;
233*b27c8092SJakub Kruzik   if (type == PC_DEFLATION_POST) {
234*b27c8092SJakub Kruzik     //pc->ops->postsolve = PCPostSolve_Deflation;
235*b27c8092SJakub Kruzik     pc->ops->presolve = 0;
236*b27c8092SJakub Kruzik   } else {
237*b27c8092SJakub Kruzik     pc->ops->presolve = PCPreSolve_Deflation;
238*b27c8092SJakub Kruzik     pc->ops->postsolve = 0;
239*b27c8092SJakub Kruzik     if (type == PC_DEFLATION_INIT) {
240*b27c8092SJakub Kruzik       def->init = PETSC_TRUE;
241*b27c8092SJakub Kruzik       pc->ops->apply = 0;
242*b27c8092SJakub Kruzik     } else {
243*b27c8092SJakub Kruzik       def->pre  = PETSC_TRUE;
244*b27c8092SJakub Kruzik     }
245*b27c8092SJakub Kruzik   }
246*b27c8092SJakub Kruzik   PetscFunctionReturn(0);
247*b27c8092SJakub Kruzik }
248*b27c8092SJakub Kruzik 
249*b27c8092SJakub Kruzik /*@
250*b27c8092SJakub Kruzik    PCDeflationSetType - Causes the deflation preconditioner to use only a special
251*b27c8092SJakub Kruzik     initial gues or pre/post solve solution update
252*b27c8092SJakub Kruzik 
253*b27c8092SJakub Kruzik    Logically Collective on PC
254*b27c8092SJakub Kruzik 
255*b27c8092SJakub Kruzik    Input Parameters:
256*b27c8092SJakub Kruzik +  pc - the preconditioner context
257*b27c8092SJakub Kruzik -  type - PC_DEFLATION_PRE, PC_DEFLATION_INIT, PC_DEFLATION_POST
258*b27c8092SJakub Kruzik 
259*b27c8092SJakub Kruzik    Options Database Key:
260*b27c8092SJakub Kruzik .  -pc_deflation_type <pre,init,post>
261*b27c8092SJakub Kruzik 
262*b27c8092SJakub Kruzik    Level: intermediate
263*b27c8092SJakub Kruzik 
264*b27c8092SJakub Kruzik    Concepts: Deflation preconditioner
265*b27c8092SJakub Kruzik 
266*b27c8092SJakub Kruzik .seealso: PCDeflationGetType()
267*b27c8092SJakub Kruzik @*/
268*b27c8092SJakub Kruzik PetscErrorCode  PCDeflationSetType(PC pc,PCDeflationType type)
269*b27c8092SJakub Kruzik {
270*b27c8092SJakub Kruzik   PetscErrorCode ierr;
271*b27c8092SJakub Kruzik 
272*b27c8092SJakub Kruzik   PetscFunctionBegin;
273*b27c8092SJakub Kruzik   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
274*b27c8092SJakub Kruzik   ierr = PetscTryMethod(pc,"PCDeflationSetType_C",(PC,PCDeflationType),(pc,type));CHKERRQ(ierr);
275*b27c8092SJakub Kruzik   PetscFunctionReturn(0);
276*b27c8092SJakub Kruzik }
277*b27c8092SJakub Kruzik 
278*b27c8092SJakub Kruzik static PetscErrorCode  PCDeflationGetType_Deflation(PC pc,PCDeflationType *type)
279*b27c8092SJakub Kruzik {
280*b27c8092SJakub Kruzik   PC_Deflation *def = (PC_Deflation*)pc->data;
281*b27c8092SJakub Kruzik 
282*b27c8092SJakub Kruzik   PetscFunctionBegin;
283*b27c8092SJakub Kruzik   if (def->init) {
284*b27c8092SJakub Kruzik     *type = PC_DEFLATION_INIT;
285*b27c8092SJakub Kruzik   } else if (def->pre) {
286*b27c8092SJakub Kruzik     *type = PC_DEFLATION_PRE;
287*b27c8092SJakub Kruzik   } else {
288*b27c8092SJakub Kruzik     *type = PC_DEFLATION_POST;
289*b27c8092SJakub Kruzik   }
290*b27c8092SJakub Kruzik   PetscFunctionReturn(0);
291*b27c8092SJakub Kruzik }
292*b27c8092SJakub Kruzik 
293*b27c8092SJakub Kruzik /*@
294*b27c8092SJakub Kruzik    PCDeflationGetType - Gets how the diagonal matrix is produced for the preconditioner
295*b27c8092SJakub Kruzik 
296*b27c8092SJakub Kruzik    Not Collective on PC
297*b27c8092SJakub Kruzik 
298*b27c8092SJakub Kruzik    Input Parameter:
299*b27c8092SJakub Kruzik .  pc - the preconditioner context
300*b27c8092SJakub Kruzik 
301*b27c8092SJakub Kruzik    Output Parameter:
302*b27c8092SJakub Kruzik -  type - PC_DEFLATION_PRE, PC_DEFLATION_INIT, PC_DEFLATION_POST
303*b27c8092SJakub Kruzik 
304*b27c8092SJakub Kruzik    Level: intermediate
305*b27c8092SJakub Kruzik 
306*b27c8092SJakub Kruzik    Concepts: Deflation preconditioner
307*b27c8092SJakub Kruzik 
308*b27c8092SJakub Kruzik .seealso: PCDeflationSetType()
309*b27c8092SJakub Kruzik @*/
310*b27c8092SJakub Kruzik PetscErrorCode  PCDeflationGetType(PC pc,PCDeflationType *type)
311*b27c8092SJakub Kruzik {
312*b27c8092SJakub Kruzik   PetscErrorCode ierr;
313*b27c8092SJakub Kruzik 
314*b27c8092SJakub Kruzik   PetscFunctionBegin;
315*b27c8092SJakub Kruzik   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
316*b27c8092SJakub Kruzik   ierr = PetscUseMethod(pc,"PCDeflationGetType_C",(PC,PCDeflationType*),(pc,type));CHKERRQ(ierr);
317*b27c8092SJakub Kruzik   PetscFunctionReturn(0);
318*b27c8092SJakub Kruzik }
319*b27c8092SJakub Kruzik 
32037eeb815SJakub Kruzik static PetscErrorCode PCSetUp_Deflation(PC pc)
32137eeb815SJakub Kruzik {
322409a357bSJakub Kruzik   PC_Deflation     *def = (PC_Deflation*)pc->data;
323409a357bSJakub Kruzik   KSP              innerksp;
324409a357bSJakub Kruzik   PC               pcinner;
325409a357bSJakub Kruzik   Mat              Amat,nextDef=NULL,*mats;
326409a357bSJakub Kruzik   PetscInt         i,m,red,size,commsize;
327409a357bSJakub Kruzik   PetscBool        match,flgspd,transp=PETSC_FALSE;
328409a357bSJakub Kruzik   MatCompositeType ctype;
329409a357bSJakub Kruzik   MPI_Comm         comm;
330409a357bSJakub Kruzik   const char       *prefix;
33137eeb815SJakub Kruzik   PetscErrorCode   ierr;
33237eeb815SJakub Kruzik 
33337eeb815SJakub Kruzik   PetscFunctionBegin;
334409a357bSJakub Kruzik   if (pc->setupcalled) PetscFunctionReturn(0);
335409a357bSJakub Kruzik   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
336f8bf229cSJakub Kruzik   ierr = PCGetOperators(pc,NULL,&Amat);CHKERRQ(ierr);
337f8bf229cSJakub Kruzik 
338f8bf229cSJakub Kruzik   /* compute a deflation space */
339409a357bSJakub Kruzik   if (def->W || def->Wt) {
340409a357bSJakub Kruzik     def->spacetype = PC_DEFLATION_SPACE_USER;
341409a357bSJakub Kruzik   } else {
342409a357bSJakub Kruzik     //ierr = KSPDCGComputeDeflationSpace(ksp);CHKERRQ(ierr);
34337eeb815SJakub Kruzik   }
34437eeb815SJakub Kruzik 
345409a357bSJakub Kruzik   /* nested deflation */
346409a357bSJakub Kruzik   if (def->W) {
347409a357bSJakub Kruzik     ierr = PetscObjectTypeCompare((PetscObject)def->W,MATCOMPOSITE,&match);CHKERRQ(ierr);
348409a357bSJakub Kruzik     if (match) {
349409a357bSJakub Kruzik       ierr = MatCompositeGetType(def->W,&ctype);CHKERRQ(ierr);
350409a357bSJakub Kruzik       ierr = MatCompositeGetNumberMat(def->W,&size);CHKERRQ(ierr);
35137eeb815SJakub Kruzik     }
352409a357bSJakub Kruzik   } else {
353409a357bSJakub Kruzik     ierr = MatCreateTranspose(def->Wt,&def->W);CHKERRQ(ierr);
354409a357bSJakub Kruzik     ierr = PetscObjectTypeCompare((PetscObject)def->Wt,MATCOMPOSITE,&match);CHKERRQ(ierr);
355409a357bSJakub Kruzik     if (match) {
356409a357bSJakub Kruzik       ierr = MatCompositeGetType(def->Wt,&ctype);CHKERRQ(ierr);
357409a357bSJakub Kruzik       ierr = MatCompositeGetNumberMat(def->Wt,&size);CHKERRQ(ierr);
358409a357bSJakub Kruzik     }
359409a357bSJakub Kruzik     transp = PETSC_TRUE;
360409a357bSJakub Kruzik   }
361409a357bSJakub Kruzik   if (match && ctype == MAT_COMPOSITE_MULTIPLICATIVE) {
362409a357bSJakub Kruzik     ierr = PetscMalloc1(size,&mats);CHKERRQ(ierr);
363409a357bSJakub Kruzik     if (!transp) {
364409a357bSJakub Kruzik       for (i=0; i<size; i++) {
365409a357bSJakub Kruzik         ierr = MatCompositeGetMat(def->W,i,&mats[i]);CHKERRQ(ierr);
366409a357bSJakub Kruzik         //ierr = PetscObjectReference((PetscObject)mats[i]);CHKERRQ(ierr);
367409a357bSJakub Kruzik       }
368409a357bSJakub Kruzik       if (def->nestedlvl < def->maxnestedlvl) {
369409a357bSJakub Kruzik         size -= 1;
370409a357bSJakub Kruzik         ierr = MatDestroy(&def->W);CHKERRQ(ierr);
371409a357bSJakub Kruzik         def->W = mats[size];
372409a357bSJakub Kruzik         ierr = PetscObjectReference((PetscObject)mats[size]);CHKERRQ(ierr);
373409a357bSJakub Kruzik         if (size > 1) {
374409a357bSJakub Kruzik           ierr = MatCreateComposite(comm,size,mats,&nextDef);CHKERRQ(ierr);
375409a357bSJakub Kruzik           ierr = MatCompositeSetType(nextDef,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr);
376409a357bSJakub Kruzik         } else {
377409a357bSJakub Kruzik           nextDef = mats[0];
378409a357bSJakub Kruzik           ierr = PetscObjectReference((PetscObject)mats[0]);CHKERRQ(ierr);
379409a357bSJakub Kruzik         }
380409a357bSJakub Kruzik       } else {
381409a357bSJakub Kruzik         /* ierr = MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT);CHKERRQ(ierr); */
382409a357bSJakub Kruzik         ierr = MatCompositeMerge(def->W);CHKERRQ(ierr);
383409a357bSJakub Kruzik       }
384409a357bSJakub Kruzik     }
385409a357bSJakub Kruzik     ierr = PetscFree(mats);CHKERRQ(ierr);
386409a357bSJakub Kruzik   }
387409a357bSJakub Kruzik 
388409a357bSJakub Kruzik   /* setup coarse problem */
389409a357bSJakub Kruzik   if (!def->WtAWinv) {
390409a357bSJakub Kruzik     ierr = MatGetSize(def->W,NULL,&m);CHKERRQ(ierr); /* TODO works for W MatTranspose? */
391409a357bSJakub Kruzik     if (!def->WtAW) {
392409a357bSJakub Kruzik       /* TODO add implicit product version ? */
393409a357bSJakub Kruzik       if (!def->AW) {
394409a357bSJakub Kruzik         ierr = MatPtAP(Amat,def->W,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtAW);CHKERRQ(ierr);
395409a357bSJakub Kruzik       } else {
396409a357bSJakub Kruzik         ierr = MatTransposeMatMult(def->W,def->AW,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtAW);CHKERRQ(ierr);
397409a357bSJakub Kruzik       }
398409a357bSJakub Kruzik       /* TODO create MatInheritOption(Mat,MatOption) */
399409a357bSJakub Kruzik       ierr = MatGetOption(Amat,MAT_SPD,&flgspd);CHKERRQ(ierr);
400409a357bSJakub Kruzik       ierr = MatSetOption(def->WtAW,MAT_SPD,flgspd);CHKERRQ(ierr);
401409a357bSJakub Kruzik #if defined(PETSC_USE_DEBUG)
402409a357bSJakub Kruzik       /* Check WtAW is not sigular */
403409a357bSJakub Kruzik       PetscReal *norms;
404409a357bSJakub Kruzik       ierr = PetscMalloc1(m,&norms);CHKERRQ(ierr);
405409a357bSJakub Kruzik       ierr = MatGetColumnNorms(def->WtAW,NORM_INFINITY,norms);CHKERRQ(ierr);
406409a357bSJakub Kruzik       for (i=0; i<m; i++) {
407409a357bSJakub Kruzik         if (norms[i] < 100*PETSC_MACHINE_EPSILON) {
408409a357bSJakub Kruzik           SETERRQ1(comm,PETSC_ERR_SUP,"Column %D of W is in kernel of A.",i);
409409a357bSJakub Kruzik         }
410409a357bSJakub Kruzik       }
411409a357bSJakub Kruzik       ierr = PetscFree(norms);CHKERRQ(ierr);
412409a357bSJakub Kruzik #endif
413409a357bSJakub Kruzik     } else {
414409a357bSJakub Kruzik       ierr = MatGetOption(def->WtAW,MAT_SPD,&flgspd);CHKERRQ(ierr);
415409a357bSJakub Kruzik     }
416409a357bSJakub Kruzik     /* TODO use MATINV */
417409a357bSJakub Kruzik     ierr = KSPCreate(comm,&def->WtAWinv);CHKERRQ(ierr);
418409a357bSJakub Kruzik     ierr = KSPSetOperators(def->WtAWinv,def->WtAW,def->WtAW);CHKERRQ(ierr);
419409a357bSJakub Kruzik     ierr = KSPSetType(def->WtAWinv,KSPPREONLY);CHKERRQ(ierr);
420409a357bSJakub Kruzik     ierr = KSPGetPC(def->WtAWinv,&pcinner);CHKERRQ(ierr);
421409a357bSJakub Kruzik     ierr = PCSetType(pcinner,PCTELESCOPE);CHKERRQ(ierr);
422409a357bSJakub Kruzik     /* Reduction factor choice */
423409a357bSJakub Kruzik     red = def->reductionfact;
424409a357bSJakub Kruzik     if (red < 0) {
425409a357bSJakub Kruzik       ierr = MPI_Comm_size(comm,&commsize);CHKERRQ(ierr);
426409a357bSJakub Kruzik       red  = ceil((float)commsize/ceil((float)m/commsize));
427409a357bSJakub Kruzik       ierr = PetscObjectTypeCompareAny((PetscObject)(def->WtAW),&match,MATSEQDENSE,MATMPIDENSE,MATDENSE,"");CHKERRQ(ierr);
428409a357bSJakub Kruzik       if (match) red = commsize;
429409a357bSJakub Kruzik       ierr = PetscInfo1(pc,"Auto choosing reduction factor %D\n",red);CHKERRQ(ierr); /* TODO add level? */
430409a357bSJakub Kruzik     }
431409a357bSJakub Kruzik     ierr = PCTelescopeSetReductionFactor(pcinner,red);CHKERRQ(ierr);
432409a357bSJakub Kruzik     ierr = PCTelescopeGetKSP(pcinner,&innerksp);CHKERRQ(ierr);
433409a357bSJakub Kruzik     ierr = KSPGetPC(innerksp,&pcinner);CHKERRQ(ierr);
434409a357bSJakub Kruzik     /* Setup KSP and PC */
435409a357bSJakub Kruzik     if (nextDef) { /* next level for multilevel deflation */
436409a357bSJakub Kruzik       /* set default KSPtype */
437409a357bSJakub Kruzik       if (!def->ksptype) {
438409a357bSJakub Kruzik         def->ksptype = KSPFGMRES;
439409a357bSJakub Kruzik         if (flgspd) { /* SPD system */
440409a357bSJakub Kruzik           def->ksptype = KSPFCG;
441409a357bSJakub Kruzik         }
442409a357bSJakub Kruzik       }
443409a357bSJakub Kruzik       ierr = KSPSetType(innerksp,def->ksptype);CHKERRQ(ierr); /* TODO iherit from KSP */
444409a357bSJakub Kruzik       ierr = PCSetType(pcinner,PCDEFLATION);CHKERRQ(ierr); /* TODO create coarse preconditinoner M_c = WtMW ? */
445409a357bSJakub Kruzik       ierr = PCDeflationSetSpace(pcinner,nextDef,transp);CHKERRQ(ierr);
4465c0c31c5SJakub Kruzik       ierr = PCDeflationSetLvl_Deflation(pcinner,def->nestedlvl+1,def->maxnestedlvl);CHKERRQ(ierr);
447409a357bSJakub Kruzik       /* inherit options TODO if not set */
448409a357bSJakub Kruzik       ((PC_Deflation*)(pcinner))->ksptype = def->ksptype;
449409a357bSJakub Kruzik       ((PC_Deflation*)(pcinner))->correct = def->correct;
450409a357bSJakub Kruzik       ((PC_Deflation*)(pcinner))->adaptiveconv = def->adaptiveconv;
451409a357bSJakub Kruzik       ((PC_Deflation*)(pcinner))->adaptiveconst = def->adaptiveconst;
452409a357bSJakub Kruzik       ierr = MatDestroy(&nextDef);CHKERRQ(ierr);
453409a357bSJakub Kruzik     } else { /* the last level */
454409a357bSJakub Kruzik       ierr = KSPSetType(innerksp,KSPPREONLY);CHKERRQ(ierr);
455409a357bSJakub Kruzik       /* TODO Cholesky if flgspd? */
456409a357bSJakub Kruzik       ierr = PCSetType(pc,PCLU);CHKERRQ(ierr);
457409a357bSJakub Kruzik       //TODO remove explicit matSolverPackage
458409a357bSJakub Kruzik       if (commsize == red) {
459409a357bSJakub Kruzik         ierr = PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU);CHKERRQ(ierr);
460409a357bSJakub Kruzik       } else {
461409a357bSJakub Kruzik         ierr = PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU_DIST);CHKERRQ(ierr);
462409a357bSJakub Kruzik       }
463409a357bSJakub Kruzik     }
464409a357bSJakub Kruzik     /* TODO use def_[lvl]_ if lvl > 0? */
465409a357bSJakub Kruzik     ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr);
466409a357bSJakub Kruzik     if (prefix) {
467409a357bSJakub Kruzik       ierr = KSPSetOptionsPrefix(innerksp,prefix);CHKERRQ(ierr);
468409a357bSJakub Kruzik       ierr = KSPAppendOptionsPrefix(innerksp,"def_");CHKERRQ(ierr);
469409a357bSJakub Kruzik     } else {
470409a357bSJakub Kruzik       ierr = KSPSetOptionsPrefix(innerksp,"def_");CHKERRQ(ierr);
471409a357bSJakub Kruzik     }
472409a357bSJakub Kruzik     /* TODO: check if WtAWinv is KSP and move following from this if */
473409a357bSJakub Kruzik     ierr = KSPSetFromOptions(def->WtAWinv);CHKERRQ(ierr);
474409a357bSJakub Kruzik     //if (def->adaptiveconv) {
475409a357bSJakub Kruzik     //  PetscReal *rnorm;
476409a357bSJakub Kruzik     //  PetscNew(&rnorm);
477409a357bSJakub Kruzik     //  ierr = KSPSetConvergenceTest(def->WtAWinv,KSPDCGConvergedAdaptive_DCG,rnorm,NULL);CHKERRQ(ierr);
478409a357bSJakub Kruzik     //}
479409a357bSJakub Kruzik     ierr = KSPSetUp(def->WtAWinv);CHKERRQ(ierr);
480409a357bSJakub Kruzik   }
481409a357bSJakub Kruzik 
482f8bf229cSJakub Kruzik   /* create work vecs */
483f8bf229cSJakub Kruzik   ierr = MatCreateVecs(Amat,NULL,&def->work);CHKERRQ(ierr);
484f8bf229cSJakub Kruzik   ierr = KSPCreateVecs(def->WtAWinv,2,&def->workcoarse,0,NULL);CHKERRQ(ierr);
48537eeb815SJakub Kruzik   PetscFunctionReturn(0);
48637eeb815SJakub Kruzik }
48737eeb815SJakub Kruzik /* -------------------------------------------------------------------------- */
48837eeb815SJakub Kruzik static PetscErrorCode PCReset_Deflation(PC pc)
48937eeb815SJakub Kruzik {
49037eeb815SJakub Kruzik   PC_Deflation      *jac = (PC_Deflation*)pc->data;
49137eeb815SJakub Kruzik   PetscErrorCode ierr;
49237eeb815SJakub Kruzik 
49337eeb815SJakub Kruzik   PetscFunctionBegin;
49437eeb815SJakub Kruzik   PetscFunctionReturn(0);
49537eeb815SJakub Kruzik }
49637eeb815SJakub Kruzik 
49737eeb815SJakub Kruzik /*
49837eeb815SJakub Kruzik    PCDestroy_Deflation - Destroys the private context for the Deflation preconditioner
49937eeb815SJakub Kruzik    that was created with PCCreate_Deflation().
50037eeb815SJakub Kruzik 
50137eeb815SJakub Kruzik    Input Parameter:
50237eeb815SJakub Kruzik .  pc - the preconditioner context
50337eeb815SJakub Kruzik 
50437eeb815SJakub Kruzik    Application Interface Routine: PCDestroy()
50537eeb815SJakub Kruzik */
50637eeb815SJakub Kruzik static PetscErrorCode PCDestroy_Deflation(PC pc)
50737eeb815SJakub Kruzik {
50837eeb815SJakub Kruzik   PetscErrorCode ierr;
50937eeb815SJakub Kruzik 
51037eeb815SJakub Kruzik   PetscFunctionBegin;
51137eeb815SJakub Kruzik   ierr = PCReset_Deflation(pc);CHKERRQ(ierr);
51237eeb815SJakub Kruzik 
51337eeb815SJakub Kruzik   /*
51437eeb815SJakub Kruzik       Free the private data structure that was hanging off the PC
51537eeb815SJakub Kruzik   */
51637eeb815SJakub Kruzik   ierr = PetscFree(pc->data);CHKERRQ(ierr);
51737eeb815SJakub Kruzik   PetscFunctionReturn(0);
51837eeb815SJakub Kruzik }
51937eeb815SJakub Kruzik 
52037eeb815SJakub Kruzik static PetscErrorCode PCSetFromOptions_Deflation(PetscOptionItems *PetscOptionsObject,PC pc)
52137eeb815SJakub Kruzik {
52237eeb815SJakub Kruzik   PC_Deflation      *jac = (PC_Deflation*)pc->data;
52337eeb815SJakub Kruzik   PetscErrorCode ierr;
52437eeb815SJakub Kruzik   PetscBool      flg;
52537eeb815SJakub Kruzik   PCDeflationType   deflt,type;
52637eeb815SJakub Kruzik 
52737eeb815SJakub Kruzik   PetscFunctionBegin;
52837eeb815SJakub Kruzik   ierr = PCDeflationGetType(pc,&deflt);CHKERRQ(ierr);
52937eeb815SJakub Kruzik   ierr = PetscOptionsHead(PetscOptionsObject,"Deflation options");CHKERRQ(ierr);
53037eeb815SJakub Kruzik   ierr = PetscOptionsEnum("-pc_jacobi_type","How to construct diagonal matrix","PCDeflationSetType",PCDeflationTypes,(PetscEnum)deflt,(PetscEnum*)&type,&flg);CHKERRQ(ierr);
53137eeb815SJakub Kruzik   if (flg) {
53237eeb815SJakub Kruzik     ierr = PCDeflationSetType(pc,type);CHKERRQ(ierr);
53337eeb815SJakub Kruzik   }
53437eeb815SJakub Kruzik   ierr = PetscOptionsTail();CHKERRQ(ierr);
53537eeb815SJakub Kruzik   PetscFunctionReturn(0);
53637eeb815SJakub Kruzik }
53737eeb815SJakub Kruzik 
53837eeb815SJakub Kruzik /*MC
539e361f8b5SJakub Kruzik      PCDEFLATION - Deflation preconditioner shifts part of the spectrum to zero (deflates)
540e361f8b5SJakub Kruzik      or to a predefined value
54137eeb815SJakub Kruzik 
54237eeb815SJakub Kruzik    Options Database Key:
543e361f8b5SJakub Kruzik +    -pc_deflation_type <init,pre,post> - selects approach to deflation (default: pre)
54437eeb815SJakub Kruzik -    -pc_jacobi_abs - use the absolute value of the diagonal entry
54537eeb815SJakub Kruzik 
54637eeb815SJakub Kruzik    Level: beginner
54737eeb815SJakub Kruzik 
54837eeb815SJakub Kruzik   Notes:
549e361f8b5SJakub Kruzik     todo
55037eeb815SJakub Kruzik 
55137eeb815SJakub Kruzik .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
552e662bc50SJakub Kruzik            PCDeflationSetType(), PCDeflationSetSpace()
55337eeb815SJakub Kruzik M*/
55437eeb815SJakub Kruzik 
55537eeb815SJakub Kruzik PETSC_EXTERN PetscErrorCode PCCreate_Deflation(PC pc)
55637eeb815SJakub Kruzik {
55737eeb815SJakub Kruzik   PC_Deflation   *def;
55837eeb815SJakub Kruzik   PetscErrorCode ierr;
55937eeb815SJakub Kruzik 
56037eeb815SJakub Kruzik   PetscFunctionBegin;
56137eeb815SJakub Kruzik   ierr     = PetscNewLog(pc,&def);CHKERRQ(ierr);
56237eeb815SJakub Kruzik   pc->data = (void*)def;
56337eeb815SJakub Kruzik 
564e662bc50SJakub Kruzik   def->init          = PETSC_FALSE;
565e662bc50SJakub Kruzik   def->pre           = PETSC_TRUE;
566e662bc50SJakub Kruzik   def->correct       = PETSC_FALSE;
567e662bc50SJakub Kruzik   def->truenorm      = PETSC_TRUE;
568e662bc50SJakub Kruzik   def->reductionfact = -1;
569e662bc50SJakub Kruzik   def->spacetype     = PC_DEFLATION_SPACE_HAAR;
570e662bc50SJakub Kruzik   def->spacesize     = 1;
571e662bc50SJakub Kruzik   def->extendsp      = PETSC_FALSE;
572e662bc50SJakub Kruzik   def->nestedlvl     = 0;
573e662bc50SJakub Kruzik   def->maxnestedlvl  = 0;
574e662bc50SJakub Kruzik   def->adaptiveconv  = PETSC_FALSE;
575e662bc50SJakub Kruzik   def->adaptiveconst = 1.0;
57637eeb815SJakub Kruzik 
57737eeb815SJakub Kruzik   /*
57837eeb815SJakub Kruzik       Set the pointers for the functions that are provided above.
57937eeb815SJakub Kruzik       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
58037eeb815SJakub Kruzik       are called, they will automatically call these functions.  Note we
58137eeb815SJakub Kruzik       choose not to provide a couple of these functions since they are
58237eeb815SJakub Kruzik       not needed.
58337eeb815SJakub Kruzik   */
58437eeb815SJakub Kruzik   pc->ops->apply               = PCApply_Deflation;
58537eeb815SJakub Kruzik   pc->ops->applytranspose      = PCApply_Deflation;
586*b27c8092SJakub Kruzik   pc->ops->presolve            = PCPreSolve_Deflation;
587*b27c8092SJakub Kruzik   pc->ops->postsolve           = 0;
58837eeb815SJakub Kruzik   pc->ops->setup               = PCSetUp_Deflation;
58937eeb815SJakub Kruzik   pc->ops->reset               = PCReset_Deflation;
59037eeb815SJakub Kruzik   pc->ops->destroy             = PCDestroy_Deflation;
59137eeb815SJakub Kruzik   pc->ops->setfromoptions      = PCSetFromOptions_Deflation;
59237eeb815SJakub Kruzik   pc->ops->view                = 0;
59337eeb815SJakub Kruzik   pc->ops->applyrichardson     = 0;
59437eeb815SJakub Kruzik   pc->ops->applysymmetricleft  = 0;
59537eeb815SJakub Kruzik   pc->ops->applysymmetricright = 0;
59637eeb815SJakub Kruzik 
59737eeb815SJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetType_C",PCDeflationSetType_Deflation);CHKERRQ(ierr);
59837eeb815SJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationGetType_C",PCDeflationGetType_Deflation);CHKERRQ(ierr);
599e662bc50SJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetSpace_C",PCDeflationSetSpace_Deflation);CHKERRQ(ierr);
6005c0c31c5SJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetLvl_C",PCDeflationSetLvl_Deflation);CHKERRQ(ierr);
60137eeb815SJakub Kruzik   PetscFunctionReturn(0);
60237eeb815SJakub Kruzik }
60337eeb815SJakub Kruzik 
604