xref: /petsc/src/ksp/pc/impls/deflation/deflation.c (revision 292e2e671b084eef5d973b8b5e75ade2ec8cff20)
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 
51*292e2e67SJakub Kruzik #include <../src/ksp/pc/impls/deflation/deflation.h> /*I "petscpc.h" I*/  /* includes for fortran wrappers */
5237eeb815SJakub Kruzik 
5337eeb815SJakub Kruzik const char *const PCDeflationTypes[]    = {"INIT","PRE","POST","PCDeflationType","PC_DEFLATION_",0};
5437eeb815SJakub Kruzik 
55e662bc50SJakub Kruzik static PetscErrorCode PCDeflationSetSpace_Deflation(PC pc,Mat W,PetscBool transpose)
56e662bc50SJakub Kruzik {
57e662bc50SJakub Kruzik   PC_Deflation   *def = (PC_Deflation*)pc->data;
58e662bc50SJakub Kruzik   PetscErrorCode ierr;
59e662bc50SJakub Kruzik 
60e662bc50SJakub Kruzik   PetscFunctionBegin;
61e662bc50SJakub Kruzik   if (transpose) {
62e662bc50SJakub Kruzik     def->Wt = W;
63e662bc50SJakub Kruzik     def->W = NULL;
64e662bc50SJakub Kruzik   } else {
65e662bc50SJakub Kruzik     def->W = W;
66e662bc50SJakub Kruzik   }
67e662bc50SJakub Kruzik   ierr = PetscObjectReference((PetscObject)W);CHKERRQ(ierr);
68e662bc50SJakub Kruzik   PetscFunctionReturn(0);
69e662bc50SJakub Kruzik }
70e662bc50SJakub Kruzik 
71e662bc50SJakub Kruzik /* TODO create PCDeflationSetSpaceTranspose? */
72e662bc50SJakub Kruzik /*@
73e662bc50SJakub Kruzik    PCDeflationSetSpace - Set deflation space matrix (or its transpose).
74e662bc50SJakub Kruzik 
75e662bc50SJakub Kruzik    Logically Collective on PC
76e662bc50SJakub Kruzik 
77e662bc50SJakub Kruzik    Input Parameters:
78e662bc50SJakub Kruzik +  pc - the preconditioner context
79e662bc50SJakub Kruzik .  W  - deflation matrix
80e662bc50SJakub Kruzik -  tranpose - indicates that W is an explicit transpose of the deflation matrix
81e662bc50SJakub Kruzik 
82e662bc50SJakub Kruzik    Level: intermediate
83e662bc50SJakub Kruzik 
84e662bc50SJakub Kruzik .seealso: PCDEFLATION
85e662bc50SJakub Kruzik @*/
86e662bc50SJakub Kruzik PetscErrorCode PCDeflationSetSpace(PC pc,Mat W,PetscBool transpose)
87e662bc50SJakub Kruzik {
88e662bc50SJakub Kruzik   PetscErrorCode ierr;
89e662bc50SJakub Kruzik 
90e662bc50SJakub Kruzik   PetscFunctionBegin;
91e662bc50SJakub Kruzik   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
92e662bc50SJakub Kruzik   PetscValidHeaderSpecific(W,MAT_CLASSID,2);
93e662bc50SJakub Kruzik   PetscValidLogicalCollectiveBool(pc,transpose,3);
94e662bc50SJakub Kruzik   ierr = PetscTryMethod(pc,"PCDeflationSetSpace_C",(PC,Mat,PetscBool),(pc,W,transpose));CHKERRQ(ierr);
95e662bc50SJakub Kruzik   PetscFunctionReturn(0);
96e662bc50SJakub Kruzik }
97e662bc50SJakub Kruzik 
985c0c31c5SJakub Kruzik static PetscErrorCode PCDeflationSetLvl_Deflation(PC pc,PetscInt current,PetscInt max)
995c0c31c5SJakub Kruzik {
1005c0c31c5SJakub Kruzik   PC_Deflation   *def = (PC_Deflation*)pc->data;
1015c0c31c5SJakub Kruzik 
1025c0c31c5SJakub Kruzik   PetscFunctionBegin;
1035c0c31c5SJakub Kruzik   def->nestedlvl = current;
1045c0c31c5SJakub Kruzik   def->maxnestedlvl = max;
1055c0c31c5SJakub Kruzik   PetscFunctionReturn(0);
1065c0c31c5SJakub Kruzik }
1075c0c31c5SJakub Kruzik 
1085c0c31c5SJakub Kruzik /*@
1095c0c31c5SJakub Kruzik    PCDeflationSetMaxLvl - Set maximum level of deflation.
1105c0c31c5SJakub Kruzik 
1115c0c31c5SJakub Kruzik    Logically Collective on PC
1125c0c31c5SJakub Kruzik 
1135c0c31c5SJakub Kruzik    Input Parameters:
1145c0c31c5SJakub Kruzik +  pc  - the preconditioner context
1155c0c31c5SJakub Kruzik .  max - maximum deflation level
1165c0c31c5SJakub Kruzik 
1175c0c31c5SJakub Kruzik    Level: intermediate
1185c0c31c5SJakub Kruzik 
1195c0c31c5SJakub Kruzik .seealso: PCDEFLATION
1205c0c31c5SJakub Kruzik @*/
1215c0c31c5SJakub Kruzik PetscErrorCode PCDeflationSetMaxLvl(PC pc,PetscInt max)
1225c0c31c5SJakub Kruzik {
1235c0c31c5SJakub Kruzik   PetscErrorCode ierr;
1245c0c31c5SJakub Kruzik 
1255c0c31c5SJakub Kruzik   PetscFunctionBegin;
1265c0c31c5SJakub Kruzik   PetscValidLogicalCollectiveInt(pc,max,2);
1275c0c31c5SJakub Kruzik   ierr = PetscTryMethod(pc,"PCDeflationSetLvl_C",(PC,PetscInt,PetscInt),(pc,0,max));CHKERRQ(ierr);
1285c0c31c5SJakub Kruzik   PetscFunctionReturn(0);
1295c0c31c5SJakub Kruzik }
130e662bc50SJakub Kruzik 
13137eeb815SJakub Kruzik /*
132b27c8092SJakub Kruzik   TODO CP corection?
133b27c8092SJakub Kruzik   x <- x + W*(W'*A*W)^{-1}*W'*r  = x + Q*r
134b27c8092SJakub Kruzik */
135b27c8092SJakub Kruzik static PetscErrorCode PCPreSolve_Deflation(PC pc,KSP ksp,Vec b, Vec x)
136b27c8092SJakub Kruzik {
137b27c8092SJakub Kruzik   PC_Deflation     *def = (PC_Deflation*)pc->data;
138b27c8092SJakub Kruzik   Mat              A;
139b27c8092SJakub Kruzik   Vec              r,w1,w2;
140b27c8092SJakub Kruzik   PetscBool        nonzero;
141b27c8092SJakub Kruzik   PetscErrorCode   ierr;
14237eeb815SJakub Kruzik 
143b27c8092SJakub Kruzik   PetscFunctionBegin;
144b27c8092SJakub Kruzik   w1 = def->workcoarse[0];
145b27c8092SJakub Kruzik   w2 = def->workcoarse[1];
146b27c8092SJakub Kruzik   r  = def->work;
147b27c8092SJakub Kruzik   ierr = PCGetOperators(pc,NULL,&A);CHKERRQ(ierr);
14837eeb815SJakub Kruzik 
149b27c8092SJakub Kruzik   ierr = KSPGetInitialGuessNonzero(ksp,&nonzero);CHKERRQ(ierr);
150b27c8092SJakub Kruzik   ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr);
151b27c8092SJakub Kruzik   if (nonzero) {
152b27c8092SJakub Kruzik     ierr = MatMult(A,x,r);CHKERRQ(ierr);               /*    r  <- b - Ax              */
153b27c8092SJakub Kruzik     ierr = VecAYPX(r,-1.0,b);CHKERRQ(ierr);
154b27c8092SJakub Kruzik   } else {
155b27c8092SJakub Kruzik     ierr = VecCopy(b,r);CHKERRQ(ierr);                 /*    r  <- b (x is 0)          */
156b27c8092SJakub Kruzik   }
157b27c8092SJakub Kruzik 
158b27c8092SJakub Kruzik   ierr = MatMultTranspose(def->W,r,w1);CHKERRQ(ierr);  /*    w1 <- W'*r                */
159b27c8092SJakub Kruzik   ierr = KSPSolve(def->WtAWinv,w1,w2);CHKERRQ(ierr);   /*    w2 <- (W'*A*W)^{-1}*w1    */
160b27c8092SJakub Kruzik   ierr = MatMult(def->W,w2,r);CHKERRQ(ierr);           /*    r  <- W*w2                */
161b27c8092SJakub Kruzik   ierr = VecAYPX(x,1.0,r);CHKERRQ(ierr);
162b27c8092SJakub Kruzik   PetscFunctionReturn(0);
163b27c8092SJakub Kruzik }
16437eeb815SJakub Kruzik 
165f8bf229cSJakub Kruzik /*
166f8bf229cSJakub Kruzik   if (def->correct) {
167f8bf229cSJakub Kruzik     z <- r - W*(W'*A*W)^{-1}*W'*(A*r -r) = (P-Q)*r
168f8bf229cSJakub Kruzik   } else {
169f8bf229cSJakub Kruzik     z <- r - W*(W'*A*W)^{-1}*W'*A*r = P*r
170f8bf229cSJakub Kruzik   }
17137eeb815SJakub Kruzik */
172f8bf229cSJakub Kruzik static PetscErrorCode PCApply_Deflation(PC pc,Vec r, Vec z)
173f8bf229cSJakub Kruzik {
174f8bf229cSJakub Kruzik   PC_Deflation     *def = (PC_Deflation*)pc->data;
175f8bf229cSJakub Kruzik   Mat              A;
176f8bf229cSJakub Kruzik   Vec              u,w1,w2;
177f8bf229cSJakub Kruzik   PetscErrorCode   ierr;
178f8bf229cSJakub Kruzik 
179f8bf229cSJakub Kruzik   PetscFunctionBegin;
180f8bf229cSJakub Kruzik   w1 = def->workcoarse[0];
181f8bf229cSJakub Kruzik   w2 = def->workcoarse[1];
182f8bf229cSJakub Kruzik   u  = def->work;
183f8bf229cSJakub Kruzik   ierr = PCGetOperators(pc,NULL,&A);CHKERRQ(ierr);
184f8bf229cSJakub Kruzik 
185f8bf229cSJakub Kruzik   if (!def->AW) {
186f8bf229cSJakub Kruzik     ierr = MatMult(A,r,u);CHKERRQ(ierr);                       /*    u  <- A*r                 */
187f8bf229cSJakub Kruzik     if (def->correct) ierr = VecAXPY(u,-1.0,r);CHKERRQ(ierr);  /*    u  <- A*r -r              */
188f8bf229cSJakub Kruzik     ierr = MatMultTranspose(def->W,u,w1);CHKERRQ(ierr);        /*    w1 <- W'*u                */
189f8bf229cSJakub Kruzik   } else {
190f8bf229cSJakub Kruzik     ierr = MatMultTranspose(def->AW,u,w1);CHKERRQ(ierr);       /*    u  <- A*r                 */
191f8bf229cSJakub Kruzik     if (def->correct) {
192f8bf229cSJakub Kruzik       ierr = MatMultTranspose(def->W,r,w2);CHKERRQ(ierr);      /*    w2 <- W'*u                */
193f8bf229cSJakub Kruzik       ierr = VecAXPY(w1,-1.0,w2);CHKERRQ(ierr);                /*    w1 <- w1 - w2             */
194f8bf229cSJakub Kruzik     }
195f8bf229cSJakub Kruzik   }
196f8bf229cSJakub Kruzik   ierr = KSPSolve(def->WtAWinv,w1,w2);CHKERRQ(ierr);           /*    w2 <- (W'*A*W)^{-1}*w1    */
197f8bf229cSJakub Kruzik   ierr = MatMult(def->W,w2,u);CHKERRQ(ierr);                   /*    u  <- W*w2                */
198f8bf229cSJakub Kruzik   ierr = VecWAXPY(z,-1.0,u,r);CHKERRQ(ierr);                   /*    z  <- r - u               */
199f8bf229cSJakub Kruzik   PetscFunctionReturn(0);
200f8bf229cSJakub Kruzik }
201f8bf229cSJakub Kruzik 
202b27c8092SJakub Kruzik static PetscErrorCode  PCDeflationSetType_Deflation(PC pc,PCDeflationType type)
203b27c8092SJakub Kruzik {
204b27c8092SJakub Kruzik   PC_Deflation *def = (PC_Deflation*)pc->data;
205b27c8092SJakub Kruzik 
206b27c8092SJakub Kruzik   PetscFunctionBegin;
207b27c8092SJakub Kruzik   def->init = PETSC_FALSE;
208b27c8092SJakub Kruzik   def->pre = PETSC_FALSE;
209b27c8092SJakub Kruzik   if (type == PC_DEFLATION_POST) {
210b27c8092SJakub Kruzik     //pc->ops->postsolve = PCPostSolve_Deflation;
211b27c8092SJakub Kruzik     pc->ops->presolve = 0;
212b27c8092SJakub Kruzik   } else {
213b27c8092SJakub Kruzik     pc->ops->presolve = PCPreSolve_Deflation;
214b27c8092SJakub Kruzik     pc->ops->postsolve = 0;
215b27c8092SJakub Kruzik     if (type == PC_DEFLATION_INIT) {
216b27c8092SJakub Kruzik       def->init = PETSC_TRUE;
217b27c8092SJakub Kruzik       pc->ops->apply = 0;
218b27c8092SJakub Kruzik     } else {
219b27c8092SJakub Kruzik       def->pre  = PETSC_TRUE;
220b27c8092SJakub Kruzik     }
221b27c8092SJakub Kruzik   }
222b27c8092SJakub Kruzik   PetscFunctionReturn(0);
223b27c8092SJakub Kruzik }
224b27c8092SJakub Kruzik 
225b27c8092SJakub Kruzik /*@
226b27c8092SJakub Kruzik    PCDeflationSetType - Causes the deflation preconditioner to use only a special
227b27c8092SJakub Kruzik     initial gues or pre/post solve solution update
228b27c8092SJakub Kruzik 
229b27c8092SJakub Kruzik    Logically Collective on PC
230b27c8092SJakub Kruzik 
231b27c8092SJakub Kruzik    Input Parameters:
232b27c8092SJakub Kruzik +  pc - the preconditioner context
233b27c8092SJakub Kruzik -  type - PC_DEFLATION_PRE, PC_DEFLATION_INIT, PC_DEFLATION_POST
234b27c8092SJakub Kruzik 
235b27c8092SJakub Kruzik    Options Database Key:
236b27c8092SJakub Kruzik .  -pc_deflation_type <pre,init,post>
237b27c8092SJakub Kruzik 
238b27c8092SJakub Kruzik    Level: intermediate
239b27c8092SJakub Kruzik 
240b27c8092SJakub Kruzik    Concepts: Deflation preconditioner
241b27c8092SJakub Kruzik 
242b27c8092SJakub Kruzik .seealso: PCDeflationGetType()
243b27c8092SJakub Kruzik @*/
244b27c8092SJakub Kruzik PetscErrorCode  PCDeflationSetType(PC pc,PCDeflationType type)
245b27c8092SJakub Kruzik {
246b27c8092SJakub Kruzik   PetscErrorCode ierr;
247b27c8092SJakub Kruzik 
248b27c8092SJakub Kruzik   PetscFunctionBegin;
249b27c8092SJakub Kruzik   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
250b27c8092SJakub Kruzik   ierr = PetscTryMethod(pc,"PCDeflationSetType_C",(PC,PCDeflationType),(pc,type));CHKERRQ(ierr);
251b27c8092SJakub Kruzik   PetscFunctionReturn(0);
252b27c8092SJakub Kruzik }
253b27c8092SJakub Kruzik 
254b27c8092SJakub Kruzik static PetscErrorCode  PCDeflationGetType_Deflation(PC pc,PCDeflationType *type)
255b27c8092SJakub Kruzik {
256b27c8092SJakub Kruzik   PC_Deflation *def = (PC_Deflation*)pc->data;
257b27c8092SJakub Kruzik 
258b27c8092SJakub Kruzik   PetscFunctionBegin;
259b27c8092SJakub Kruzik   if (def->init) {
260b27c8092SJakub Kruzik     *type = PC_DEFLATION_INIT;
261b27c8092SJakub Kruzik   } else if (def->pre) {
262b27c8092SJakub Kruzik     *type = PC_DEFLATION_PRE;
263b27c8092SJakub Kruzik   } else {
264b27c8092SJakub Kruzik     *type = PC_DEFLATION_POST;
265b27c8092SJakub Kruzik   }
266b27c8092SJakub Kruzik   PetscFunctionReturn(0);
267b27c8092SJakub Kruzik }
268b27c8092SJakub Kruzik 
269b27c8092SJakub Kruzik /*@
270b27c8092SJakub Kruzik    PCDeflationGetType - Gets how the diagonal matrix is produced for the preconditioner
271b27c8092SJakub Kruzik 
272b27c8092SJakub Kruzik    Not Collective on PC
273b27c8092SJakub Kruzik 
274b27c8092SJakub Kruzik    Input Parameter:
275b27c8092SJakub Kruzik .  pc - the preconditioner context
276b27c8092SJakub Kruzik 
277b27c8092SJakub Kruzik    Output Parameter:
278b27c8092SJakub Kruzik -  type - PC_DEFLATION_PRE, PC_DEFLATION_INIT, PC_DEFLATION_POST
279b27c8092SJakub Kruzik 
280b27c8092SJakub Kruzik    Level: intermediate
281b27c8092SJakub Kruzik 
282b27c8092SJakub Kruzik    Concepts: Deflation preconditioner
283b27c8092SJakub Kruzik 
284b27c8092SJakub Kruzik .seealso: PCDeflationSetType()
285b27c8092SJakub Kruzik @*/
286b27c8092SJakub Kruzik PetscErrorCode  PCDeflationGetType(PC pc,PCDeflationType *type)
287b27c8092SJakub Kruzik {
288b27c8092SJakub Kruzik   PetscErrorCode ierr;
289b27c8092SJakub Kruzik 
290b27c8092SJakub Kruzik   PetscFunctionBegin;
291b27c8092SJakub Kruzik   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
292b27c8092SJakub Kruzik   ierr = PetscUseMethod(pc,"PCDeflationGetType_C",(PC,PCDeflationType*),(pc,type));CHKERRQ(ierr);
293b27c8092SJakub Kruzik   PetscFunctionReturn(0);
294b27c8092SJakub Kruzik }
295b27c8092SJakub Kruzik 
29637eeb815SJakub Kruzik static PetscErrorCode PCSetUp_Deflation(PC pc)
29737eeb815SJakub Kruzik {
298409a357bSJakub Kruzik   PC_Deflation     *def = (PC_Deflation*)pc->data;
299409a357bSJakub Kruzik   KSP              innerksp;
300409a357bSJakub Kruzik   PC               pcinner;
301409a357bSJakub Kruzik   Mat              Amat,nextDef=NULL,*mats;
302409a357bSJakub Kruzik   PetscInt         i,m,red,size,commsize;
303409a357bSJakub Kruzik   PetscBool        match,flgspd,transp=PETSC_FALSE;
304409a357bSJakub Kruzik   MatCompositeType ctype;
305409a357bSJakub Kruzik   MPI_Comm         comm;
306409a357bSJakub Kruzik   const char       *prefix;
30737eeb815SJakub Kruzik   PetscErrorCode   ierr;
30837eeb815SJakub Kruzik 
30937eeb815SJakub Kruzik   PetscFunctionBegin;
310409a357bSJakub Kruzik   if (pc->setupcalled) PetscFunctionReturn(0);
311409a357bSJakub Kruzik   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
312f8bf229cSJakub Kruzik   ierr = PCGetOperators(pc,NULL,&Amat);CHKERRQ(ierr);
313f8bf229cSJakub Kruzik 
314f8bf229cSJakub Kruzik   /* compute a deflation space */
315409a357bSJakub Kruzik   if (def->W || def->Wt) {
316409a357bSJakub Kruzik     def->spacetype = PC_DEFLATION_SPACE_USER;
317409a357bSJakub Kruzik   } else {
318409a357bSJakub Kruzik     //ierr = KSPDCGComputeDeflationSpace(ksp);CHKERRQ(ierr);
31937eeb815SJakub Kruzik   }
32037eeb815SJakub Kruzik 
321409a357bSJakub Kruzik   /* nested deflation */
322409a357bSJakub Kruzik   if (def->W) {
323409a357bSJakub Kruzik     ierr = PetscObjectTypeCompare((PetscObject)def->W,MATCOMPOSITE,&match);CHKERRQ(ierr);
324409a357bSJakub Kruzik     if (match) {
325409a357bSJakub Kruzik       ierr = MatCompositeGetType(def->W,&ctype);CHKERRQ(ierr);
326409a357bSJakub Kruzik       ierr = MatCompositeGetNumberMat(def->W,&size);CHKERRQ(ierr);
32737eeb815SJakub Kruzik     }
328409a357bSJakub Kruzik   } else {
329409a357bSJakub Kruzik     ierr = MatCreateTranspose(def->Wt,&def->W);CHKERRQ(ierr);
330409a357bSJakub Kruzik     ierr = PetscObjectTypeCompare((PetscObject)def->Wt,MATCOMPOSITE,&match);CHKERRQ(ierr);
331409a357bSJakub Kruzik     if (match) {
332409a357bSJakub Kruzik       ierr = MatCompositeGetType(def->Wt,&ctype);CHKERRQ(ierr);
333409a357bSJakub Kruzik       ierr = MatCompositeGetNumberMat(def->Wt,&size);CHKERRQ(ierr);
334409a357bSJakub Kruzik     }
335409a357bSJakub Kruzik     transp = PETSC_TRUE;
336409a357bSJakub Kruzik   }
337409a357bSJakub Kruzik   if (match && ctype == MAT_COMPOSITE_MULTIPLICATIVE) {
338409a357bSJakub Kruzik     ierr = PetscMalloc1(size,&mats);CHKERRQ(ierr);
339409a357bSJakub Kruzik     if (!transp) {
340409a357bSJakub Kruzik       for (i=0; i<size; i++) {
341409a357bSJakub Kruzik         ierr = MatCompositeGetMat(def->W,i,&mats[i]);CHKERRQ(ierr);
342409a357bSJakub Kruzik         //ierr = PetscObjectReference((PetscObject)mats[i]);CHKERRQ(ierr);
343409a357bSJakub Kruzik       }
344409a357bSJakub Kruzik       if (def->nestedlvl < def->maxnestedlvl) {
345409a357bSJakub Kruzik         size -= 1;
346409a357bSJakub Kruzik         ierr = MatDestroy(&def->W);CHKERRQ(ierr);
347409a357bSJakub Kruzik         def->W = mats[size];
348409a357bSJakub Kruzik         ierr = PetscObjectReference((PetscObject)mats[size]);CHKERRQ(ierr);
349409a357bSJakub Kruzik         if (size > 1) {
350409a357bSJakub Kruzik           ierr = MatCreateComposite(comm,size,mats,&nextDef);CHKERRQ(ierr);
351409a357bSJakub Kruzik           ierr = MatCompositeSetType(nextDef,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr);
352409a357bSJakub Kruzik         } else {
353409a357bSJakub Kruzik           nextDef = mats[0];
354409a357bSJakub Kruzik           ierr = PetscObjectReference((PetscObject)mats[0]);CHKERRQ(ierr);
355409a357bSJakub Kruzik         }
356409a357bSJakub Kruzik       } else {
357409a357bSJakub Kruzik         /* ierr = MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT);CHKERRQ(ierr); */
358409a357bSJakub Kruzik         ierr = MatCompositeMerge(def->W);CHKERRQ(ierr);
359409a357bSJakub Kruzik       }
360409a357bSJakub Kruzik     }
361409a357bSJakub Kruzik     ierr = PetscFree(mats);CHKERRQ(ierr);
362409a357bSJakub Kruzik   }
363409a357bSJakub Kruzik 
364409a357bSJakub Kruzik   /* setup coarse problem */
365409a357bSJakub Kruzik   if (!def->WtAWinv) {
366409a357bSJakub Kruzik     ierr = MatGetSize(def->W,NULL,&m);CHKERRQ(ierr); /* TODO works for W MatTranspose? */
367409a357bSJakub Kruzik     if (!def->WtAW) {
368409a357bSJakub Kruzik       /* TODO add implicit product version ? */
369409a357bSJakub Kruzik       if (!def->AW) {
370409a357bSJakub Kruzik         ierr = MatPtAP(Amat,def->W,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtAW);CHKERRQ(ierr);
371409a357bSJakub Kruzik       } else {
372409a357bSJakub Kruzik         ierr = MatTransposeMatMult(def->W,def->AW,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtAW);CHKERRQ(ierr);
373409a357bSJakub Kruzik       }
374409a357bSJakub Kruzik       /* TODO create MatInheritOption(Mat,MatOption) */
375409a357bSJakub Kruzik       ierr = MatGetOption(Amat,MAT_SPD,&flgspd);CHKERRQ(ierr);
376409a357bSJakub Kruzik       ierr = MatSetOption(def->WtAW,MAT_SPD,flgspd);CHKERRQ(ierr);
377409a357bSJakub Kruzik #if defined(PETSC_USE_DEBUG)
378409a357bSJakub Kruzik       /* Check WtAW is not sigular */
379409a357bSJakub Kruzik       PetscReal *norms;
380409a357bSJakub Kruzik       ierr = PetscMalloc1(m,&norms);CHKERRQ(ierr);
381409a357bSJakub Kruzik       ierr = MatGetColumnNorms(def->WtAW,NORM_INFINITY,norms);CHKERRQ(ierr);
382409a357bSJakub Kruzik       for (i=0; i<m; i++) {
383409a357bSJakub Kruzik         if (norms[i] < 100*PETSC_MACHINE_EPSILON) {
384409a357bSJakub Kruzik           SETERRQ1(comm,PETSC_ERR_SUP,"Column %D of W is in kernel of A.",i);
385409a357bSJakub Kruzik         }
386409a357bSJakub Kruzik       }
387409a357bSJakub Kruzik       ierr = PetscFree(norms);CHKERRQ(ierr);
388409a357bSJakub Kruzik #endif
389409a357bSJakub Kruzik     } else {
390409a357bSJakub Kruzik       ierr = MatGetOption(def->WtAW,MAT_SPD,&flgspd);CHKERRQ(ierr);
391409a357bSJakub Kruzik     }
392409a357bSJakub Kruzik     /* TODO use MATINV */
393ac66f006SJakub Kruzik     ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr);
394409a357bSJakub Kruzik     ierr = KSPCreate(comm,&def->WtAWinv);CHKERRQ(ierr);
395409a357bSJakub Kruzik     ierr = KSPSetOperators(def->WtAWinv,def->WtAW,def->WtAW);CHKERRQ(ierr);
396409a357bSJakub Kruzik     ierr = KSPSetType(def->WtAWinv,KSPPREONLY);CHKERRQ(ierr);
397409a357bSJakub Kruzik     ierr = KSPGetPC(def->WtAWinv,&pcinner);CHKERRQ(ierr);
398409a357bSJakub Kruzik     ierr = PCSetType(pcinner,PCTELESCOPE);CHKERRQ(ierr);
399ac66f006SJakub Kruzik     /* ugly hack to not have overwritten PCTELESCOPE */
400ac66f006SJakub Kruzik     if (prefix) {
401ac66f006SJakub Kruzik       ierr = KSPSetOptionsPrefix(def->WtAWinv,prefix);CHKERRQ(ierr);
402ac66f006SJakub Kruzik       ierr = KSPAppendOptionsPrefix(def->WtAWinv,"tel_");CHKERRQ(ierr);
403ac66f006SJakub Kruzik     } else {
404ac66f006SJakub Kruzik       ierr = KSPSetOptionsPrefix(def->WtAWinv,"tel_");CHKERRQ(ierr);
405ac66f006SJakub Kruzik     }
406ac66f006SJakub Kruzik     ierr = PCSetFromOptions(pcinner);CHKERRQ(ierr);
407409a357bSJakub Kruzik     /* Reduction factor choice */
408409a357bSJakub Kruzik     red = def->reductionfact;
409409a357bSJakub Kruzik     if (red < 0) {
410409a357bSJakub Kruzik       ierr = MPI_Comm_size(comm,&commsize);CHKERRQ(ierr);
411409a357bSJakub Kruzik       red  = ceil((float)commsize/ceil((float)m/commsize));
412409a357bSJakub Kruzik       ierr = PetscObjectTypeCompareAny((PetscObject)(def->WtAW),&match,MATSEQDENSE,MATMPIDENSE,MATDENSE,"");CHKERRQ(ierr);
413409a357bSJakub Kruzik       if (match) red = commsize;
414409a357bSJakub Kruzik       ierr = PetscInfo1(pc,"Auto choosing reduction factor %D\n",red);CHKERRQ(ierr); /* TODO add level? */
415409a357bSJakub Kruzik     }
416409a357bSJakub Kruzik     ierr = PCTelescopeSetReductionFactor(pcinner,red);CHKERRQ(ierr);
417ac66f006SJakub Kruzik     ierr = PCSetUp(pcinner);CHKERRQ(ierr);
418409a357bSJakub Kruzik     ierr = PCTelescopeGetKSP(pcinner,&innerksp);CHKERRQ(ierr);
419ac66f006SJakub Kruzik     if (innerksp) {
420409a357bSJakub Kruzik       ierr = KSPGetPC(innerksp,&pcinner);CHKERRQ(ierr);
421409a357bSJakub Kruzik       /* Setup KSP and PC */
422409a357bSJakub Kruzik       if (nextDef) { /* next level for multilevel deflation */
423409a357bSJakub Kruzik         /* set default KSPtype */
424409a357bSJakub Kruzik         if (!def->ksptype) {
425409a357bSJakub Kruzik           def->ksptype = KSPFGMRES;
426409a357bSJakub Kruzik           if (flgspd) { /* SPD system */
427409a357bSJakub Kruzik             def->ksptype = KSPFCG;
428409a357bSJakub Kruzik           }
429409a357bSJakub Kruzik         }
430409a357bSJakub Kruzik         ierr = KSPSetType(innerksp,def->ksptype);CHKERRQ(ierr); /* TODO iherit from KSP */
431409a357bSJakub Kruzik         ierr = PCSetType(pcinner,PCDEFLATION);CHKERRQ(ierr); /* TODO create coarse preconditinoner M_c = WtMW ? */
432409a357bSJakub Kruzik         ierr = PCDeflationSetSpace(pcinner,nextDef,transp);CHKERRQ(ierr);
4335c0c31c5SJakub Kruzik         ierr = PCDeflationSetLvl_Deflation(pcinner,def->nestedlvl+1,def->maxnestedlvl);CHKERRQ(ierr);
434409a357bSJakub Kruzik         /* inherit options TODO if not set */
435409a357bSJakub Kruzik         ((PC_Deflation*)(pcinner))->ksptype = def->ksptype;
436409a357bSJakub Kruzik         ((PC_Deflation*)(pcinner))->correct = def->correct;
437409a357bSJakub Kruzik         ((PC_Deflation*)(pcinner))->adaptiveconv = def->adaptiveconv;
438409a357bSJakub Kruzik         ((PC_Deflation*)(pcinner))->adaptiveconst = def->adaptiveconst;
439409a357bSJakub Kruzik         ierr = MatDestroy(&nextDef);CHKERRQ(ierr);
440409a357bSJakub Kruzik       } else { /* the last level */
441ac66f006SJakub Kruzik         ierr = KSPSetType(innerksp,KSPGMRES);CHKERRQ(ierr);
442409a357bSJakub Kruzik         /* TODO Cholesky if flgspd? */
443ac66f006SJakub Kruzik         ierr = PCSetType(pcinner,PCLU);CHKERRQ(ierr);
444409a357bSJakub Kruzik         //TODO remove explicit matSolverPackage
445409a357bSJakub Kruzik         if (commsize == red) {
446ac66f006SJakub Kruzik           ierr = PCFactorSetMatSolverType(pcinner,MATSOLVERSUPERLU);CHKERRQ(ierr);
447409a357bSJakub Kruzik         } else {
448ac66f006SJakub Kruzik           ierr = PCFactorSetMatSolverType(pcinner,MATSOLVERSUPERLU_DIST);CHKERRQ(ierr);
449409a357bSJakub Kruzik         }
450409a357bSJakub Kruzik       }
451409a357bSJakub Kruzik       /* TODO use def_[lvl]_ if lvl > 0? */
452409a357bSJakub Kruzik       if (prefix) {
453409a357bSJakub Kruzik         ierr = KSPSetOptionsPrefix(innerksp,prefix);CHKERRQ(ierr);
454409a357bSJakub Kruzik         ierr = KSPAppendOptionsPrefix(innerksp,"def_");CHKERRQ(ierr);
455409a357bSJakub Kruzik       } else {
456409a357bSJakub Kruzik         ierr = KSPSetOptionsPrefix(innerksp,"def_");CHKERRQ(ierr);
457409a357bSJakub Kruzik       }
458ac66f006SJakub Kruzik     }
459409a357bSJakub Kruzik     /* TODO: check if WtAWinv is KSP and move following from this if */
460409a357bSJakub Kruzik     ierr = KSPSetFromOptions(def->WtAWinv);CHKERRQ(ierr);
461409a357bSJakub Kruzik     //if (def->adaptiveconv) {
462409a357bSJakub Kruzik     //  PetscReal *rnorm;
463409a357bSJakub Kruzik     //  PetscNew(&rnorm);
464409a357bSJakub Kruzik     //  ierr = KSPSetConvergenceTest(def->WtAWinv,KSPDCGConvergedAdaptive_DCG,rnorm,NULL);CHKERRQ(ierr);
465409a357bSJakub Kruzik     //}
466409a357bSJakub Kruzik     ierr = KSPSetUp(def->WtAWinv);CHKERRQ(ierr);
467409a357bSJakub Kruzik   }
468409a357bSJakub Kruzik 
469f8bf229cSJakub Kruzik   /* create work vecs */
470f8bf229cSJakub Kruzik   ierr = MatCreateVecs(Amat,NULL,&def->work);CHKERRQ(ierr);
471f8bf229cSJakub Kruzik   ierr = KSPCreateVecs(def->WtAWinv,2,&def->workcoarse,0,NULL);CHKERRQ(ierr);
47237eeb815SJakub Kruzik   PetscFunctionReturn(0);
47337eeb815SJakub Kruzik }
47437eeb815SJakub Kruzik /* -------------------------------------------------------------------------- */
47537eeb815SJakub Kruzik static PetscErrorCode PCReset_Deflation(PC pc)
47637eeb815SJakub Kruzik {
47737eeb815SJakub Kruzik   PC_Deflation      *jac = (PC_Deflation*)pc->data;
47837eeb815SJakub Kruzik   PetscErrorCode ierr;
47937eeb815SJakub Kruzik 
48037eeb815SJakub Kruzik   PetscFunctionBegin;
48137eeb815SJakub Kruzik   PetscFunctionReturn(0);
48237eeb815SJakub Kruzik }
48337eeb815SJakub Kruzik 
48437eeb815SJakub Kruzik /*
48537eeb815SJakub Kruzik    PCDestroy_Deflation - Destroys the private context for the Deflation preconditioner
48637eeb815SJakub Kruzik    that was created with PCCreate_Deflation().
48737eeb815SJakub Kruzik 
48837eeb815SJakub Kruzik    Input Parameter:
48937eeb815SJakub Kruzik .  pc - the preconditioner context
49037eeb815SJakub Kruzik 
49137eeb815SJakub Kruzik    Application Interface Routine: PCDestroy()
49237eeb815SJakub Kruzik */
49337eeb815SJakub Kruzik static PetscErrorCode PCDestroy_Deflation(PC pc)
49437eeb815SJakub Kruzik {
49537eeb815SJakub Kruzik   PetscErrorCode ierr;
49637eeb815SJakub Kruzik 
49737eeb815SJakub Kruzik   PetscFunctionBegin;
49837eeb815SJakub Kruzik   ierr = PCReset_Deflation(pc);CHKERRQ(ierr);
49937eeb815SJakub Kruzik 
50037eeb815SJakub Kruzik   /*
50137eeb815SJakub Kruzik       Free the private data structure that was hanging off the PC
50237eeb815SJakub Kruzik   */
50337eeb815SJakub Kruzik   ierr = PetscFree(pc->data);CHKERRQ(ierr);
50437eeb815SJakub Kruzik   PetscFunctionReturn(0);
50537eeb815SJakub Kruzik }
50637eeb815SJakub Kruzik 
50737eeb815SJakub Kruzik static PetscErrorCode PCSetFromOptions_Deflation(PetscOptionItems *PetscOptionsObject,PC pc)
50837eeb815SJakub Kruzik {
50937eeb815SJakub Kruzik   PC_Deflation      *jac = (PC_Deflation*)pc->data;
51037eeb815SJakub Kruzik   PetscErrorCode ierr;
51137eeb815SJakub Kruzik   PetscBool      flg;
51237eeb815SJakub Kruzik   PCDeflationType   deflt,type;
51337eeb815SJakub Kruzik 
51437eeb815SJakub Kruzik   PetscFunctionBegin;
51537eeb815SJakub Kruzik   ierr = PCDeflationGetType(pc,&deflt);CHKERRQ(ierr);
51637eeb815SJakub Kruzik   ierr = PetscOptionsHead(PetscOptionsObject,"Deflation options");CHKERRQ(ierr);
51737eeb815SJakub Kruzik   ierr = PetscOptionsEnum("-pc_jacobi_type","How to construct diagonal matrix","PCDeflationSetType",PCDeflationTypes,(PetscEnum)deflt,(PetscEnum*)&type,&flg);CHKERRQ(ierr);
51837eeb815SJakub Kruzik   if (flg) {
51937eeb815SJakub Kruzik     ierr = PCDeflationSetType(pc,type);CHKERRQ(ierr);
52037eeb815SJakub Kruzik   }
52137eeb815SJakub Kruzik   ierr = PetscOptionsTail();CHKERRQ(ierr);
52237eeb815SJakub Kruzik   PetscFunctionReturn(0);
52337eeb815SJakub Kruzik }
52437eeb815SJakub Kruzik 
52537eeb815SJakub Kruzik /*MC
526e361f8b5SJakub Kruzik      PCDEFLATION - Deflation preconditioner shifts part of the spectrum to zero (deflates)
527e361f8b5SJakub Kruzik      or to a predefined value
52837eeb815SJakub Kruzik 
52937eeb815SJakub Kruzik    Options Database Key:
530e361f8b5SJakub Kruzik +    -pc_deflation_type <init,pre,post> - selects approach to deflation (default: pre)
53137eeb815SJakub Kruzik -    -pc_jacobi_abs - use the absolute value of the diagonal entry
53237eeb815SJakub Kruzik 
53337eeb815SJakub Kruzik    Level: beginner
53437eeb815SJakub Kruzik 
53537eeb815SJakub Kruzik   Notes:
536e361f8b5SJakub Kruzik     todo
53737eeb815SJakub Kruzik 
53837eeb815SJakub Kruzik .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
539e662bc50SJakub Kruzik            PCDeflationSetType(), PCDeflationSetSpace()
54037eeb815SJakub Kruzik M*/
54137eeb815SJakub Kruzik 
54237eeb815SJakub Kruzik PETSC_EXTERN PetscErrorCode PCCreate_Deflation(PC pc)
54337eeb815SJakub Kruzik {
54437eeb815SJakub Kruzik   PC_Deflation   *def;
54537eeb815SJakub Kruzik   PetscErrorCode ierr;
54637eeb815SJakub Kruzik 
54737eeb815SJakub Kruzik   PetscFunctionBegin;
54837eeb815SJakub Kruzik   ierr     = PetscNewLog(pc,&def);CHKERRQ(ierr);
54937eeb815SJakub Kruzik   pc->data = (void*)def;
55037eeb815SJakub Kruzik 
551e662bc50SJakub Kruzik   def->init          = PETSC_FALSE;
552e662bc50SJakub Kruzik   def->pre           = PETSC_TRUE;
553e662bc50SJakub Kruzik   def->correct       = PETSC_FALSE;
554e662bc50SJakub Kruzik   def->truenorm      = PETSC_TRUE;
555e662bc50SJakub Kruzik   def->reductionfact = -1;
556e662bc50SJakub Kruzik   def->spacetype     = PC_DEFLATION_SPACE_HAAR;
557e662bc50SJakub Kruzik   def->spacesize     = 1;
558e662bc50SJakub Kruzik   def->extendsp      = PETSC_FALSE;
559e662bc50SJakub Kruzik   def->nestedlvl     = 0;
560e662bc50SJakub Kruzik   def->maxnestedlvl  = 0;
561e662bc50SJakub Kruzik   def->adaptiveconv  = PETSC_FALSE;
562e662bc50SJakub Kruzik   def->adaptiveconst = 1.0;
56337eeb815SJakub Kruzik 
56437eeb815SJakub Kruzik   /*
56537eeb815SJakub Kruzik       Set the pointers for the functions that are provided above.
56637eeb815SJakub Kruzik       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
56737eeb815SJakub Kruzik       are called, they will automatically call these functions.  Note we
56837eeb815SJakub Kruzik       choose not to provide a couple of these functions since they are
56937eeb815SJakub Kruzik       not needed.
57037eeb815SJakub Kruzik   */
57137eeb815SJakub Kruzik   pc->ops->apply               = PCApply_Deflation;
57237eeb815SJakub Kruzik   pc->ops->applytranspose      = PCApply_Deflation;
573b27c8092SJakub Kruzik   pc->ops->presolve            = PCPreSolve_Deflation;
574b27c8092SJakub Kruzik   pc->ops->postsolve           = 0;
57537eeb815SJakub Kruzik   pc->ops->setup               = PCSetUp_Deflation;
57637eeb815SJakub Kruzik   pc->ops->reset               = PCReset_Deflation;
57737eeb815SJakub Kruzik   pc->ops->destroy             = PCDestroy_Deflation;
57837eeb815SJakub Kruzik   pc->ops->setfromoptions      = PCSetFromOptions_Deflation;
57937eeb815SJakub Kruzik   pc->ops->view                = 0;
58037eeb815SJakub Kruzik   pc->ops->applyrichardson     = 0;
58137eeb815SJakub Kruzik   pc->ops->applysymmetricleft  = 0;
58237eeb815SJakub Kruzik   pc->ops->applysymmetricright = 0;
58337eeb815SJakub Kruzik 
58437eeb815SJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetType_C",PCDeflationSetType_Deflation);CHKERRQ(ierr);
58537eeb815SJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationGetType_C",PCDeflationGetType_Deflation);CHKERRQ(ierr);
586e662bc50SJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetSpace_C",PCDeflationSetSpace_Deflation);CHKERRQ(ierr);
5875c0c31c5SJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetLvl_C",PCDeflationSetLvl_Deflation);CHKERRQ(ierr);
58837eeb815SJakub Kruzik   PetscFunctionReturn(0);
58937eeb815SJakub Kruzik }
59037eeb815SJakub Kruzik 
591