xref: /petsc/src/ksp/pc/impls/deflation/deflation.c (revision f8bf229ca243c5c3d537374664be453d1ec7953a)
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;
69*f8bf229cSJakub Kruzik   Vec       work;
70*f8bf229cSJakub 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 
7937eeb815SJakub Kruzik static PetscErrorCode  PCDeflationSetType_Deflation(PC pc,PCDeflationType type)
8037eeb815SJakub Kruzik {
8137eeb815SJakub Kruzik   PC_Deflation *def = (PC_Deflation*)pc->data;
8237eeb815SJakub Kruzik 
8337eeb815SJakub Kruzik   PetscFunctionBegin;
8437eeb815SJakub Kruzik   def->init = PETSC_FALSE;
8537eeb815SJakub Kruzik   def->pre = PETSC_FALSE;
8637eeb815SJakub Kruzik   if (type == PC_DEFLATION_INIT) {
8737eeb815SJakub Kruzik     def->init = PETSC_TRUE;
8837eeb815SJakub Kruzik     def->pre  = PETSC_TRUE;
8937eeb815SJakub Kruzik   } else if (type == PC_DEFLATION_PRE) {
9037eeb815SJakub Kruzik     def->pre  = PETSC_TRUE;
9137eeb815SJakub Kruzik   }
9237eeb815SJakub Kruzik   PetscFunctionReturn(0);
9337eeb815SJakub Kruzik }
9437eeb815SJakub Kruzik 
95bb9edbfaSJakub Kruzik /*@
96bb9edbfaSJakub Kruzik    PCDeflationSetType - Causes the deflation preconditioner to use only a special
97bb9edbfaSJakub Kruzik     initial gues or pre/post solve solution update
98bb9edbfaSJakub Kruzik 
99bb9edbfaSJakub Kruzik    Logically Collective on PC
100bb9edbfaSJakub Kruzik 
101bb9edbfaSJakub Kruzik    Input Parameters:
102bb9edbfaSJakub Kruzik +  pc - the preconditioner context
103bb9edbfaSJakub Kruzik -  type - PC_DEFLATION_PRE, PC_DEFLATION_INIT, PC_DEFLATION_POST
104bb9edbfaSJakub Kruzik 
105bb9edbfaSJakub Kruzik    Options Database Key:
106bb9edbfaSJakub Kruzik .  -pc_deflation_type <pre,init,post>
107bb9edbfaSJakub Kruzik 
108bb9edbfaSJakub Kruzik    Level: intermediate
109bb9edbfaSJakub Kruzik 
110bb9edbfaSJakub Kruzik    Concepts: Deflation preconditioner
111bb9edbfaSJakub Kruzik 
112bb9edbfaSJakub Kruzik .seealso: PCDeflationGetType()
113bb9edbfaSJakub Kruzik @*/
114bb9edbfaSJakub Kruzik PetscErrorCode  PCDeflationSetType(PC pc,PCDeflationType type)
115bb9edbfaSJakub Kruzik {
116bb9edbfaSJakub Kruzik   PetscErrorCode ierr;
117bb9edbfaSJakub Kruzik 
118bb9edbfaSJakub Kruzik   PetscFunctionBegin;
119bb9edbfaSJakub Kruzik   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
120bb9edbfaSJakub Kruzik   ierr = PetscTryMethod(pc,"PCDeflationSetType_C",(PC,PCDeflationType),(pc,type));CHKERRQ(ierr);
121bb9edbfaSJakub Kruzik   PetscFunctionReturn(0);
122bb9edbfaSJakub Kruzik }
123bb9edbfaSJakub Kruzik 
12437eeb815SJakub Kruzik static PetscErrorCode  PCDeflationGetType_Deflation(PC pc,PCDeflationType *type)
12537eeb815SJakub Kruzik {
12637eeb815SJakub Kruzik   PC_Deflation *def = (PC_Deflation*)pc->data;
12737eeb815SJakub Kruzik 
12837eeb815SJakub Kruzik   PetscFunctionBegin;
12937eeb815SJakub Kruzik   if (def->init) {
13037eeb815SJakub Kruzik     *type = PC_DEFLATION_INIT;
13137eeb815SJakub Kruzik   } else if (def->pre) {
13237eeb815SJakub Kruzik     *type = PC_DEFLATION_PRE;
13337eeb815SJakub Kruzik   } else {
13437eeb815SJakub Kruzik     *type = PC_DEFLATION_POST;
13537eeb815SJakub Kruzik   }
13637eeb815SJakub Kruzik   PetscFunctionReturn(0);
13737eeb815SJakub Kruzik }
13837eeb815SJakub Kruzik 
139bb9edbfaSJakub Kruzik /*@
140bb9edbfaSJakub Kruzik    PCDeflationGetType - Gets how the diagonal matrix is produced for the preconditioner
141bb9edbfaSJakub Kruzik 
142bb9edbfaSJakub Kruzik    Not Collective on PC
143bb9edbfaSJakub Kruzik 
144bb9edbfaSJakub Kruzik    Input Parameter:
145bb9edbfaSJakub Kruzik .  pc - the preconditioner context
146bb9edbfaSJakub Kruzik 
147bb9edbfaSJakub Kruzik    Output Parameter:
148bb9edbfaSJakub Kruzik -  type - PC_DEFLATION_PRE, PC_DEFLATION_INIT, PC_DEFLATION_POST
149bb9edbfaSJakub Kruzik 
150bb9edbfaSJakub Kruzik    Level: intermediate
151bb9edbfaSJakub Kruzik 
152bb9edbfaSJakub Kruzik    Concepts: Deflation preconditioner
153bb9edbfaSJakub Kruzik 
154bb9edbfaSJakub Kruzik .seealso: PCDeflationSetType()
155bb9edbfaSJakub Kruzik @*/
156bb9edbfaSJakub Kruzik PetscErrorCode  PCDeflationGetType(PC pc,PCDeflationType *type)
157bb9edbfaSJakub Kruzik {
158bb9edbfaSJakub Kruzik   PetscErrorCode ierr;
159bb9edbfaSJakub Kruzik 
160bb9edbfaSJakub Kruzik   PetscFunctionBegin;
161bb9edbfaSJakub Kruzik   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
162bb9edbfaSJakub Kruzik   ierr = PetscUseMethod(pc,"PCDeflationGetType_C",(PC,PCDeflationType*),(pc,type));CHKERRQ(ierr);
163bb9edbfaSJakub Kruzik   PetscFunctionReturn(0);
164bb9edbfaSJakub Kruzik }
165bb9edbfaSJakub Kruzik 
166e662bc50SJakub Kruzik static PetscErrorCode PCDeflationSetSpace_Deflation(PC pc,Mat W,PetscBool transpose)
167e662bc50SJakub Kruzik {
168e662bc50SJakub Kruzik   PC_Deflation   *def = (PC_Deflation*)pc->data;
169e662bc50SJakub Kruzik   PetscErrorCode ierr;
170e662bc50SJakub Kruzik 
171e662bc50SJakub Kruzik   PetscFunctionBegin;
172e662bc50SJakub Kruzik   if (transpose) {
173e662bc50SJakub Kruzik     def->Wt = W;
174e662bc50SJakub Kruzik     def->W = NULL;
175e662bc50SJakub Kruzik   } else {
176e662bc50SJakub Kruzik     def->W = W;
177e662bc50SJakub Kruzik   }
178e662bc50SJakub Kruzik   ierr = PetscObjectReference((PetscObject)W);CHKERRQ(ierr);
179e662bc50SJakub Kruzik   PetscFunctionReturn(0);
180e662bc50SJakub Kruzik }
181e662bc50SJakub Kruzik 
182e662bc50SJakub Kruzik /* TODO create PCDeflationSetSpaceTranspose? */
183e662bc50SJakub Kruzik /*@
184e662bc50SJakub Kruzik    PCDeflationSetSpace - Set deflation space matrix (or its transpose).
185e662bc50SJakub Kruzik 
186e662bc50SJakub Kruzik    Logically Collective on PC
187e662bc50SJakub Kruzik 
188e662bc50SJakub Kruzik    Input Parameters:
189e662bc50SJakub Kruzik +  pc - the preconditioner context
190e662bc50SJakub Kruzik .  W  - deflation matrix
191e662bc50SJakub Kruzik -  tranpose - indicates that W is an explicit transpose of the deflation matrix
192e662bc50SJakub Kruzik 
193e662bc50SJakub Kruzik    Level: intermediate
194e662bc50SJakub Kruzik 
195e662bc50SJakub Kruzik .seealso: PCDEFLATION
196e662bc50SJakub Kruzik @*/
197e662bc50SJakub Kruzik PetscErrorCode PCDeflationSetSpace(PC pc,Mat W,PetscBool transpose)
198e662bc50SJakub Kruzik {
199e662bc50SJakub Kruzik   PetscErrorCode ierr;
200e662bc50SJakub Kruzik 
201e662bc50SJakub Kruzik   PetscFunctionBegin;
202e662bc50SJakub Kruzik   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
203e662bc50SJakub Kruzik   PetscValidHeaderSpecific(W,MAT_CLASSID,2);
204e662bc50SJakub Kruzik   PetscValidLogicalCollectiveBool(pc,transpose,3);
205e662bc50SJakub Kruzik   ierr = PetscTryMethod(pc,"PCDeflationSetSpace_C",(PC,Mat,PetscBool),(pc,W,transpose));CHKERRQ(ierr);
206e662bc50SJakub Kruzik   PetscFunctionReturn(0);
207e662bc50SJakub Kruzik }
208e662bc50SJakub Kruzik 
2095c0c31c5SJakub Kruzik static PetscErrorCode PCDeflationSetLvl_Deflation(PC pc,PetscInt current,PetscInt max)
2105c0c31c5SJakub Kruzik {
2115c0c31c5SJakub Kruzik   PC_Deflation   *def = (PC_Deflation*)pc->data;
2125c0c31c5SJakub Kruzik 
2135c0c31c5SJakub Kruzik   PetscFunctionBegin;
2145c0c31c5SJakub Kruzik   def->nestedlvl = current;
2155c0c31c5SJakub Kruzik   def->maxnestedlvl = max;
2165c0c31c5SJakub Kruzik   PetscFunctionReturn(0);
2175c0c31c5SJakub Kruzik }
2185c0c31c5SJakub Kruzik 
2195c0c31c5SJakub Kruzik /*@
2205c0c31c5SJakub Kruzik    PCDeflationSetMaxLvl - Set maximum level of deflation.
2215c0c31c5SJakub Kruzik 
2225c0c31c5SJakub Kruzik    Logically Collective on PC
2235c0c31c5SJakub Kruzik 
2245c0c31c5SJakub Kruzik    Input Parameters:
2255c0c31c5SJakub Kruzik +  pc  - the preconditioner context
2265c0c31c5SJakub Kruzik .  max - maximum deflation level
2275c0c31c5SJakub Kruzik 
2285c0c31c5SJakub Kruzik    Level: intermediate
2295c0c31c5SJakub Kruzik 
2305c0c31c5SJakub Kruzik .seealso: PCDEFLATION
2315c0c31c5SJakub Kruzik @*/
2325c0c31c5SJakub Kruzik PetscErrorCode PCDeflationSetMaxLvl(PC pc,PetscInt max)
2335c0c31c5SJakub Kruzik {
2345c0c31c5SJakub Kruzik   PetscErrorCode ierr;
2355c0c31c5SJakub Kruzik 
2365c0c31c5SJakub Kruzik   PetscFunctionBegin;
2375c0c31c5SJakub Kruzik   PetscValidLogicalCollectiveInt(pc,max,2);
2385c0c31c5SJakub Kruzik   ierr = PetscTryMethod(pc,"PCDeflationSetLvl_C",(PC,PetscInt,PetscInt),(pc,0,max));CHKERRQ(ierr);
2395c0c31c5SJakub Kruzik   PetscFunctionReturn(0);
2405c0c31c5SJakub Kruzik }
241e662bc50SJakub Kruzik 
24237eeb815SJakub Kruzik /*
24337eeb815SJakub Kruzik    PCSetUp_Deflation - Prepares for the use of the Deflation preconditioner
24437eeb815SJakub Kruzik                     by setting data structures and options.
24537eeb815SJakub Kruzik 
24637eeb815SJakub Kruzik    Input Parameter:
24737eeb815SJakub Kruzik .  pc - the preconditioner context
24837eeb815SJakub Kruzik 
24937eeb815SJakub Kruzik    Application Interface Routine: PCSetUp()
25037eeb815SJakub Kruzik 
251*f8bf229cSJakub Kruzik /*
252*f8bf229cSJakub Kruzik   if (def->correct) {
253*f8bf229cSJakub Kruzik     z <- r - W*(W'*A*W)^{-1}*W'*(A*r -r) = (P-Q)*r
254*f8bf229cSJakub Kruzik   } else {
255*f8bf229cSJakub Kruzik     z <- r - W*(W'*A*W)^{-1}*W'*A*r = P*r
256*f8bf229cSJakub Kruzik   }
25737eeb815SJakub Kruzik */
258*f8bf229cSJakub Kruzik static PetscErrorCode PCApply_Deflation(PC pc,Vec r, Vec z)
259*f8bf229cSJakub Kruzik {
260*f8bf229cSJakub Kruzik   PC_Deflation     *def = (PC_Deflation*)pc->data;
261*f8bf229cSJakub Kruzik   Mat              A;
262*f8bf229cSJakub Kruzik   Vec              u,w1,w2;
263*f8bf229cSJakub Kruzik   PetscErrorCode   ierr;
264*f8bf229cSJakub Kruzik 
265*f8bf229cSJakub Kruzik   PetscFunctionBegin;
266*f8bf229cSJakub Kruzik   w1 = def->workcoarse[0];
267*f8bf229cSJakub Kruzik   w2 = def->workcoarse[1];
268*f8bf229cSJakub Kruzik   u  = def->work;
269*f8bf229cSJakub Kruzik   ierr = PCGetOperators(pc,NULL,&A);CHKERRQ(ierr);
270*f8bf229cSJakub Kruzik 
271*f8bf229cSJakub Kruzik   if (!def->AW) {
272*f8bf229cSJakub Kruzik     ierr = MatMult(A,r,u);CHKERRQ(ierr);                       /*    u  <- A*r                 */
273*f8bf229cSJakub Kruzik     if (def->correct) ierr = VecAXPY(u,-1.0,r);CHKERRQ(ierr);  /*    u  <- A*r -r              */
274*f8bf229cSJakub Kruzik     ierr = MatMultTranspose(def->W,u,w1);CHKERRQ(ierr);        /*    w1 <- W'*u                */
275*f8bf229cSJakub Kruzik   } else {
276*f8bf229cSJakub Kruzik     ierr = MatMultTranspose(def->AW,u,w1);CHKERRQ(ierr);       /*    u  <- A*r                 */
277*f8bf229cSJakub Kruzik     if (def->correct) {
278*f8bf229cSJakub Kruzik       ierr = MatMultTranspose(def->W,r,w2);CHKERRQ(ierr);      /*    w2 <- W'*u                */
279*f8bf229cSJakub Kruzik       ierr = VecAXPY(w1,-1.0,w2);CHKERRQ(ierr);                /*    w1 <- w1 - w2             */
280*f8bf229cSJakub Kruzik     }
281*f8bf229cSJakub Kruzik   }
282*f8bf229cSJakub Kruzik   ierr = KSPSolve(def->WtAWinv,w1,w2);CHKERRQ(ierr);           /*    w2 <- (W'*A*W)^{-1}*w1    */
283*f8bf229cSJakub Kruzik   ierr = MatMult(def->W,w2,u);CHKERRQ(ierr);                   /*    u  <- W*w2                */
284*f8bf229cSJakub Kruzik   ierr = VecWAXPY(z,-1.0,u,r);CHKERRQ(ierr);                   /*    z  <- r - u               */
285*f8bf229cSJakub Kruzik   PetscFunctionReturn(0);
286*f8bf229cSJakub Kruzik }
287*f8bf229cSJakub Kruzik 
28837eeb815SJakub Kruzik static PetscErrorCode PCSetUp_Deflation(PC pc)
28937eeb815SJakub Kruzik {
290409a357bSJakub Kruzik   PC_Deflation     *def = (PC_Deflation*)pc->data;
291409a357bSJakub Kruzik   KSP              innerksp;
292409a357bSJakub Kruzik   PC               pcinner;
293409a357bSJakub Kruzik   Mat              Amat,nextDef=NULL,*mats;
294409a357bSJakub Kruzik   PetscInt         i,m,red,size,commsize;
295409a357bSJakub Kruzik   PetscBool        match,flgspd,transp=PETSC_FALSE;
296409a357bSJakub Kruzik   MatCompositeType ctype;
297409a357bSJakub Kruzik   MPI_Comm         comm;
298409a357bSJakub Kruzik   const char       *prefix;
29937eeb815SJakub Kruzik   PetscErrorCode   ierr;
30037eeb815SJakub Kruzik 
30137eeb815SJakub Kruzik   PetscFunctionBegin;
302409a357bSJakub Kruzik   if (pc->setupcalled) PetscFunctionReturn(0);
303409a357bSJakub Kruzik   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
304*f8bf229cSJakub Kruzik   ierr = PCGetOperators(pc,NULL,&Amat);CHKERRQ(ierr);
305*f8bf229cSJakub Kruzik 
306*f8bf229cSJakub Kruzik   /* compute a deflation space */
307409a357bSJakub Kruzik   if (def->W || def->Wt) {
308409a357bSJakub Kruzik     def->spacetype = PC_DEFLATION_SPACE_USER;
309409a357bSJakub Kruzik   } else {
310409a357bSJakub Kruzik     //ierr = KSPDCGComputeDeflationSpace(ksp);CHKERRQ(ierr);
31137eeb815SJakub Kruzik   }
31237eeb815SJakub Kruzik 
313409a357bSJakub Kruzik   /* nested deflation */
314409a357bSJakub Kruzik   if (def->W) {
315409a357bSJakub Kruzik     ierr = PetscObjectTypeCompare((PetscObject)def->W,MATCOMPOSITE,&match);CHKERRQ(ierr);
316409a357bSJakub Kruzik     if (match) {
317409a357bSJakub Kruzik       ierr = MatCompositeGetType(def->W,&ctype);CHKERRQ(ierr);
318409a357bSJakub Kruzik       ierr = MatCompositeGetNumberMat(def->W,&size);CHKERRQ(ierr);
31937eeb815SJakub Kruzik     }
320409a357bSJakub Kruzik   } else {
321409a357bSJakub Kruzik     ierr = MatCreateTranspose(def->Wt,&def->W);CHKERRQ(ierr);
322409a357bSJakub Kruzik     ierr = PetscObjectTypeCompare((PetscObject)def->Wt,MATCOMPOSITE,&match);CHKERRQ(ierr);
323409a357bSJakub Kruzik     if (match) {
324409a357bSJakub Kruzik       ierr = MatCompositeGetType(def->Wt,&ctype);CHKERRQ(ierr);
325409a357bSJakub Kruzik       ierr = MatCompositeGetNumberMat(def->Wt,&size);CHKERRQ(ierr);
326409a357bSJakub Kruzik     }
327409a357bSJakub Kruzik     transp = PETSC_TRUE;
328409a357bSJakub Kruzik   }
329409a357bSJakub Kruzik   if (match && ctype == MAT_COMPOSITE_MULTIPLICATIVE) {
330409a357bSJakub Kruzik     ierr = PetscMalloc1(size,&mats);CHKERRQ(ierr);
331409a357bSJakub Kruzik     if (!transp) {
332409a357bSJakub Kruzik       for (i=0; i<size; i++) {
333409a357bSJakub Kruzik         ierr = MatCompositeGetMat(def->W,i,&mats[i]);CHKERRQ(ierr);
334409a357bSJakub Kruzik         //ierr = PetscObjectReference((PetscObject)mats[i]);CHKERRQ(ierr);
335409a357bSJakub Kruzik       }
336409a357bSJakub Kruzik       if (def->nestedlvl < def->maxnestedlvl) {
337409a357bSJakub Kruzik         size -= 1;
338409a357bSJakub Kruzik         ierr = MatDestroy(&def->W);CHKERRQ(ierr);
339409a357bSJakub Kruzik         def->W = mats[size];
340409a357bSJakub Kruzik         ierr = PetscObjectReference((PetscObject)mats[size]);CHKERRQ(ierr);
341409a357bSJakub Kruzik         if (size > 1) {
342409a357bSJakub Kruzik           ierr = MatCreateComposite(comm,size,mats,&nextDef);CHKERRQ(ierr);
343409a357bSJakub Kruzik           ierr = MatCompositeSetType(nextDef,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr);
344409a357bSJakub Kruzik         } else {
345409a357bSJakub Kruzik           nextDef = mats[0];
346409a357bSJakub Kruzik           ierr = PetscObjectReference((PetscObject)mats[0]);CHKERRQ(ierr);
347409a357bSJakub Kruzik         }
348409a357bSJakub Kruzik       } else {
349409a357bSJakub Kruzik         /* ierr = MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT);CHKERRQ(ierr); */
350409a357bSJakub Kruzik         ierr = MatCompositeMerge(def->W);CHKERRQ(ierr);
351409a357bSJakub Kruzik       }
352409a357bSJakub Kruzik     }
353409a357bSJakub Kruzik     ierr = PetscFree(mats);CHKERRQ(ierr);
354409a357bSJakub Kruzik   }
355409a357bSJakub Kruzik 
356409a357bSJakub Kruzik   /* setup coarse problem */
357409a357bSJakub Kruzik   if (!def->WtAWinv) {
358409a357bSJakub Kruzik     ierr = MatGetSize(def->W,NULL,&m);CHKERRQ(ierr); /* TODO works for W MatTranspose? */
359409a357bSJakub Kruzik     if (!def->WtAW) {
360409a357bSJakub Kruzik       /* TODO add implicit product version ? */
361409a357bSJakub Kruzik       if (!def->AW) {
362409a357bSJakub Kruzik         ierr = MatPtAP(Amat,def->W,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtAW);CHKERRQ(ierr);
363409a357bSJakub Kruzik       } else {
364409a357bSJakub Kruzik         ierr = MatTransposeMatMult(def->W,def->AW,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtAW);CHKERRQ(ierr);
365409a357bSJakub Kruzik       }
366409a357bSJakub Kruzik       /* TODO create MatInheritOption(Mat,MatOption) */
367409a357bSJakub Kruzik       ierr = MatGetOption(Amat,MAT_SPD,&flgspd);CHKERRQ(ierr);
368409a357bSJakub Kruzik       ierr = MatSetOption(def->WtAW,MAT_SPD,flgspd);CHKERRQ(ierr);
369409a357bSJakub Kruzik #if defined(PETSC_USE_DEBUG)
370409a357bSJakub Kruzik       /* Check WtAW is not sigular */
371409a357bSJakub Kruzik       PetscReal *norms;
372409a357bSJakub Kruzik       ierr = PetscMalloc1(m,&norms);CHKERRQ(ierr);
373409a357bSJakub Kruzik       ierr = MatGetColumnNorms(def->WtAW,NORM_INFINITY,norms);CHKERRQ(ierr);
374409a357bSJakub Kruzik       for (i=0; i<m; i++) {
375409a357bSJakub Kruzik         if (norms[i] < 100*PETSC_MACHINE_EPSILON) {
376409a357bSJakub Kruzik           SETERRQ1(comm,PETSC_ERR_SUP,"Column %D of W is in kernel of A.",i);
377409a357bSJakub Kruzik         }
378409a357bSJakub Kruzik       }
379409a357bSJakub Kruzik       ierr = PetscFree(norms);CHKERRQ(ierr);
380409a357bSJakub Kruzik #endif
381409a357bSJakub Kruzik     } else {
382409a357bSJakub Kruzik       ierr = MatGetOption(def->WtAW,MAT_SPD,&flgspd);CHKERRQ(ierr);
383409a357bSJakub Kruzik     }
384409a357bSJakub Kruzik     /* TODO use MATINV */
385409a357bSJakub Kruzik     ierr = KSPCreate(comm,&def->WtAWinv);CHKERRQ(ierr);
386409a357bSJakub Kruzik     ierr = KSPSetOperators(def->WtAWinv,def->WtAW,def->WtAW);CHKERRQ(ierr);
387409a357bSJakub Kruzik     ierr = KSPSetType(def->WtAWinv,KSPPREONLY);CHKERRQ(ierr);
388409a357bSJakub Kruzik     ierr = KSPGetPC(def->WtAWinv,&pcinner);CHKERRQ(ierr);
389409a357bSJakub Kruzik     ierr = PCSetType(pcinner,PCTELESCOPE);CHKERRQ(ierr);
390409a357bSJakub Kruzik     /* Reduction factor choice */
391409a357bSJakub Kruzik     red = def->reductionfact;
392409a357bSJakub Kruzik     if (red < 0) {
393409a357bSJakub Kruzik       ierr = MPI_Comm_size(comm,&commsize);CHKERRQ(ierr);
394409a357bSJakub Kruzik       red  = ceil((float)commsize/ceil((float)m/commsize));
395409a357bSJakub Kruzik       ierr = PetscObjectTypeCompareAny((PetscObject)(def->WtAW),&match,MATSEQDENSE,MATMPIDENSE,MATDENSE,"");CHKERRQ(ierr);
396409a357bSJakub Kruzik       if (match) red = commsize;
397409a357bSJakub Kruzik       ierr = PetscInfo1(pc,"Auto choosing reduction factor %D\n",red);CHKERRQ(ierr); /* TODO add level? */
398409a357bSJakub Kruzik     }
399409a357bSJakub Kruzik     ierr = PCTelescopeSetReductionFactor(pcinner,red);CHKERRQ(ierr);
400409a357bSJakub Kruzik     ierr = PCTelescopeGetKSP(pcinner,&innerksp);CHKERRQ(ierr);
401409a357bSJakub Kruzik     ierr = KSPGetPC(innerksp,&pcinner);CHKERRQ(ierr);
402409a357bSJakub Kruzik     /* Setup KSP and PC */
403409a357bSJakub Kruzik     if (nextDef) { /* next level for multilevel deflation */
404409a357bSJakub Kruzik       /* set default KSPtype */
405409a357bSJakub Kruzik       if (!def->ksptype) {
406409a357bSJakub Kruzik         def->ksptype = KSPFGMRES;
407409a357bSJakub Kruzik         if (flgspd) { /* SPD system */
408409a357bSJakub Kruzik           def->ksptype = KSPFCG;
409409a357bSJakub Kruzik         }
410409a357bSJakub Kruzik       }
411409a357bSJakub Kruzik       ierr = KSPSetType(innerksp,def->ksptype);CHKERRQ(ierr); /* TODO iherit from KSP */
412409a357bSJakub Kruzik       ierr = PCSetType(pcinner,PCDEFLATION);CHKERRQ(ierr); /* TODO create coarse preconditinoner M_c = WtMW ? */
413409a357bSJakub Kruzik       ierr = PCDeflationSetSpace(pcinner,nextDef,transp);CHKERRQ(ierr);
4145c0c31c5SJakub Kruzik       ierr = PCDeflationSetLvl_Deflation(pcinner,def->nestedlvl+1,def->maxnestedlvl);CHKERRQ(ierr);
415409a357bSJakub Kruzik       /* inherit options TODO if not set */
416409a357bSJakub Kruzik       ((PC_Deflation*)(pcinner))->ksptype = def->ksptype;
417409a357bSJakub Kruzik       ((PC_Deflation*)(pcinner))->correct = def->correct;
418409a357bSJakub Kruzik       ((PC_Deflation*)(pcinner))->adaptiveconv = def->adaptiveconv;
419409a357bSJakub Kruzik       ((PC_Deflation*)(pcinner))->adaptiveconst = def->adaptiveconst;
420409a357bSJakub Kruzik       ierr = MatDestroy(&nextDef);CHKERRQ(ierr);
421409a357bSJakub Kruzik     } else { /* the last level */
422409a357bSJakub Kruzik       ierr = KSPSetType(innerksp,KSPPREONLY);CHKERRQ(ierr);
423409a357bSJakub Kruzik       /* TODO Cholesky if flgspd? */
424409a357bSJakub Kruzik       ierr = PCSetType(pc,PCLU);CHKERRQ(ierr);
425409a357bSJakub Kruzik       //TODO remove explicit matSolverPackage
426409a357bSJakub Kruzik       if (commsize == red) {
427409a357bSJakub Kruzik         ierr = PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU);CHKERRQ(ierr);
428409a357bSJakub Kruzik       } else {
429409a357bSJakub Kruzik         ierr = PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU_DIST);CHKERRQ(ierr);
430409a357bSJakub Kruzik       }
431409a357bSJakub Kruzik     }
432409a357bSJakub Kruzik     /* TODO use def_[lvl]_ if lvl > 0? */
433409a357bSJakub Kruzik     ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr);
434409a357bSJakub Kruzik     if (prefix) {
435409a357bSJakub Kruzik       ierr = KSPSetOptionsPrefix(innerksp,prefix);CHKERRQ(ierr);
436409a357bSJakub Kruzik       ierr = KSPAppendOptionsPrefix(innerksp,"def_");CHKERRQ(ierr);
437409a357bSJakub Kruzik     } else {
438409a357bSJakub Kruzik       ierr = KSPSetOptionsPrefix(innerksp,"def_");CHKERRQ(ierr);
439409a357bSJakub Kruzik     }
440409a357bSJakub Kruzik     /* TODO: check if WtAWinv is KSP and move following from this if */
441409a357bSJakub Kruzik     ierr = KSPSetFromOptions(def->WtAWinv);CHKERRQ(ierr);
442409a357bSJakub Kruzik     //if (def->adaptiveconv) {
443409a357bSJakub Kruzik     //  PetscReal *rnorm;
444409a357bSJakub Kruzik     //  PetscNew(&rnorm);
445409a357bSJakub Kruzik     //  ierr = KSPSetConvergenceTest(def->WtAWinv,KSPDCGConvergedAdaptive_DCG,rnorm,NULL);CHKERRQ(ierr);
446409a357bSJakub Kruzik     //}
447409a357bSJakub Kruzik     ierr = KSPSetUp(def->WtAWinv);CHKERRQ(ierr);
448409a357bSJakub Kruzik   }
449409a357bSJakub Kruzik 
450*f8bf229cSJakub Kruzik   /* create work vecs */
451*f8bf229cSJakub Kruzik   ierr = MatCreateVecs(Amat,NULL,&def->work);CHKERRQ(ierr);
452*f8bf229cSJakub Kruzik   ierr = KSPCreateVecs(def->WtAWinv,2,&def->workcoarse,0,NULL);CHKERRQ(ierr);
45337eeb815SJakub Kruzik   PetscFunctionReturn(0);
45437eeb815SJakub Kruzik }
45537eeb815SJakub Kruzik /* -------------------------------------------------------------------------- */
45637eeb815SJakub Kruzik static PetscErrorCode PCReset_Deflation(PC pc)
45737eeb815SJakub Kruzik {
45837eeb815SJakub Kruzik   PC_Deflation      *jac = (PC_Deflation*)pc->data;
45937eeb815SJakub Kruzik   PetscErrorCode ierr;
46037eeb815SJakub Kruzik 
46137eeb815SJakub Kruzik   PetscFunctionBegin;
46237eeb815SJakub Kruzik   PetscFunctionReturn(0);
46337eeb815SJakub Kruzik }
46437eeb815SJakub Kruzik 
46537eeb815SJakub Kruzik /*
46637eeb815SJakub Kruzik    PCDestroy_Deflation - Destroys the private context for the Deflation preconditioner
46737eeb815SJakub Kruzik    that was created with PCCreate_Deflation().
46837eeb815SJakub Kruzik 
46937eeb815SJakub Kruzik    Input Parameter:
47037eeb815SJakub Kruzik .  pc - the preconditioner context
47137eeb815SJakub Kruzik 
47237eeb815SJakub Kruzik    Application Interface Routine: PCDestroy()
47337eeb815SJakub Kruzik */
47437eeb815SJakub Kruzik static PetscErrorCode PCDestroy_Deflation(PC pc)
47537eeb815SJakub Kruzik {
47637eeb815SJakub Kruzik   PetscErrorCode ierr;
47737eeb815SJakub Kruzik 
47837eeb815SJakub Kruzik   PetscFunctionBegin;
47937eeb815SJakub Kruzik   ierr = PCReset_Deflation(pc);CHKERRQ(ierr);
48037eeb815SJakub Kruzik 
48137eeb815SJakub Kruzik   /*
48237eeb815SJakub Kruzik       Free the private data structure that was hanging off the PC
48337eeb815SJakub Kruzik   */
48437eeb815SJakub Kruzik   ierr = PetscFree(pc->data);CHKERRQ(ierr);
48537eeb815SJakub Kruzik   PetscFunctionReturn(0);
48637eeb815SJakub Kruzik }
48737eeb815SJakub Kruzik 
48837eeb815SJakub Kruzik static PetscErrorCode PCSetFromOptions_Deflation(PetscOptionItems *PetscOptionsObject,PC pc)
48937eeb815SJakub Kruzik {
49037eeb815SJakub Kruzik   PC_Deflation      *jac = (PC_Deflation*)pc->data;
49137eeb815SJakub Kruzik   PetscErrorCode ierr;
49237eeb815SJakub Kruzik   PetscBool      flg;
49337eeb815SJakub Kruzik   PCDeflationType   deflt,type;
49437eeb815SJakub Kruzik 
49537eeb815SJakub Kruzik   PetscFunctionBegin;
49637eeb815SJakub Kruzik   ierr = PCDeflationGetType(pc,&deflt);CHKERRQ(ierr);
49737eeb815SJakub Kruzik   ierr = PetscOptionsHead(PetscOptionsObject,"Deflation options");CHKERRQ(ierr);
49837eeb815SJakub Kruzik   ierr = PetscOptionsEnum("-pc_jacobi_type","How to construct diagonal matrix","PCDeflationSetType",PCDeflationTypes,(PetscEnum)deflt,(PetscEnum*)&type,&flg);CHKERRQ(ierr);
49937eeb815SJakub Kruzik   if (flg) {
50037eeb815SJakub Kruzik     ierr = PCDeflationSetType(pc,type);CHKERRQ(ierr);
50137eeb815SJakub Kruzik   }
50237eeb815SJakub Kruzik   ierr = PetscOptionsTail();CHKERRQ(ierr);
50337eeb815SJakub Kruzik   PetscFunctionReturn(0);
50437eeb815SJakub Kruzik }
50537eeb815SJakub Kruzik 
50637eeb815SJakub Kruzik /*MC
507e361f8b5SJakub Kruzik      PCDEFLATION - Deflation preconditioner shifts part of the spectrum to zero (deflates)
508e361f8b5SJakub Kruzik      or to a predefined value
50937eeb815SJakub Kruzik 
51037eeb815SJakub Kruzik    Options Database Key:
511e361f8b5SJakub Kruzik +    -pc_deflation_type <init,pre,post> - selects approach to deflation (default: pre)
51237eeb815SJakub Kruzik -    -pc_jacobi_abs - use the absolute value of the diagonal entry
51337eeb815SJakub Kruzik 
51437eeb815SJakub Kruzik    Level: beginner
51537eeb815SJakub Kruzik 
51637eeb815SJakub Kruzik   Notes:
517e361f8b5SJakub Kruzik     todo
51837eeb815SJakub Kruzik 
51937eeb815SJakub Kruzik .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
520e662bc50SJakub Kruzik            PCDeflationSetType(), PCDeflationSetSpace()
52137eeb815SJakub Kruzik M*/
52237eeb815SJakub Kruzik 
52337eeb815SJakub Kruzik PETSC_EXTERN PetscErrorCode PCCreate_Deflation(PC pc)
52437eeb815SJakub Kruzik {
52537eeb815SJakub Kruzik   PC_Deflation   *def;
52637eeb815SJakub Kruzik   PetscErrorCode ierr;
52737eeb815SJakub Kruzik 
52837eeb815SJakub Kruzik   PetscFunctionBegin;
52937eeb815SJakub Kruzik   ierr     = PetscNewLog(pc,&def);CHKERRQ(ierr);
53037eeb815SJakub Kruzik   pc->data = (void*)def;
53137eeb815SJakub Kruzik 
532e662bc50SJakub Kruzik   def->init          = PETSC_FALSE;
533e662bc50SJakub Kruzik   def->pre           = PETSC_TRUE;
534e662bc50SJakub Kruzik   def->correct       = PETSC_FALSE;
535e662bc50SJakub Kruzik   def->truenorm      = PETSC_TRUE;
536e662bc50SJakub Kruzik   def->reductionfact = -1;
537e662bc50SJakub Kruzik   def->spacetype     = PC_DEFLATION_SPACE_HAAR;
538e662bc50SJakub Kruzik   def->spacesize     = 1;
539e662bc50SJakub Kruzik   def->extendsp      = PETSC_FALSE;
540e662bc50SJakub Kruzik   def->nestedlvl     = 0;
541e662bc50SJakub Kruzik   def->maxnestedlvl  = 0;
542e662bc50SJakub Kruzik   def->adaptiveconv  = PETSC_FALSE;
543e662bc50SJakub Kruzik   def->adaptiveconst = 1.0;
54437eeb815SJakub Kruzik 
54537eeb815SJakub Kruzik   /*
54637eeb815SJakub Kruzik       Set the pointers for the functions that are provided above.
54737eeb815SJakub Kruzik       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
54837eeb815SJakub Kruzik       are called, they will automatically call these functions.  Note we
54937eeb815SJakub Kruzik       choose not to provide a couple of these functions since they are
55037eeb815SJakub Kruzik       not needed.
55137eeb815SJakub Kruzik   */
55237eeb815SJakub Kruzik   pc->ops->apply               = PCApply_Deflation;
55337eeb815SJakub Kruzik   pc->ops->applytranspose      = PCApply_Deflation;
55437eeb815SJakub Kruzik   pc->ops->setup               = PCSetUp_Deflation;
55537eeb815SJakub Kruzik   pc->ops->reset               = PCReset_Deflation;
55637eeb815SJakub Kruzik   pc->ops->destroy             = PCDestroy_Deflation;
55737eeb815SJakub Kruzik   pc->ops->setfromoptions      = PCSetFromOptions_Deflation;
55837eeb815SJakub Kruzik   pc->ops->view                = 0;
55937eeb815SJakub Kruzik   pc->ops->applyrichardson     = 0;
56037eeb815SJakub Kruzik   pc->ops->applysymmetricleft  = 0;
56137eeb815SJakub Kruzik   pc->ops->applysymmetricright = 0;
56237eeb815SJakub Kruzik 
56337eeb815SJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetType_C",PCDeflationSetType_Deflation);CHKERRQ(ierr);
56437eeb815SJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationGetType_C",PCDeflationGetType_Deflation);CHKERRQ(ierr);
565e662bc50SJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetSpace_C",PCDeflationSetSpace_Deflation);CHKERRQ(ierr);
5665c0c31c5SJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetLvl_C",PCDeflationSetLvl_Deflation);CHKERRQ(ierr);
56737eeb815SJakub Kruzik   PetscFunctionReturn(0);
56837eeb815SJakub Kruzik }
56937eeb815SJakub Kruzik 
570