xref: /petsc/src/ksp/pc/impls/deflation/deflation.c (revision 9e56ec8ad35c371f54e154e3c50c5316551eb8e2)
1292e2e67SJakub Kruzik #include <../src/ksp/pc/impls/deflation/deflation.h> /*I "petscpc.h" I*/  /* includes for fortran wrappers */
237eeb815SJakub Kruzik 
38513960bSJakub Kruzik const char *const PCDeflationSpaceTypes[] = {
48513960bSJakub Kruzik   "haar",
58513960bSJakub Kruzik   "db2",
68513960bSJakub Kruzik   "db4",
78513960bSJakub Kruzik   "db8",
88513960bSJakub Kruzik   "db16",
98513960bSJakub Kruzik   "biorth22",
108513960bSJakub Kruzik   "meyer",
118513960bSJakub Kruzik   "aggregation",
128513960bSJakub Kruzik   "slepc",
138513960bSJakub Kruzik   "slepc-cheap",
148513960bSJakub Kruzik   "user",
150a78af22SJakub Kruzik   "PCDeflationSpaceType",
160a78af22SJakub Kruzik   "PC_DEFLATION_SPACE_",
178513960bSJakub Kruzik   0
188513960bSJakub Kruzik };
198513960bSJakub Kruzik 
20a122ebaeSJakub Kruzik static PetscErrorCode PCDeflationSetInitOnly_Deflation(PC pc,PetscBool flg)
21a122ebaeSJakub Kruzik {
22a122ebaeSJakub Kruzik   PC_Deflation   *def = (PC_Deflation*)pc->data;
23a122ebaeSJakub Kruzik 
24a122ebaeSJakub Kruzik   PetscFunctionBegin;
25a122ebaeSJakub Kruzik   def->init = flg;
26a122ebaeSJakub Kruzik   PetscFunctionReturn(0);
27a122ebaeSJakub Kruzik }
28a122ebaeSJakub Kruzik 
29a122ebaeSJakub Kruzik /*@
30a122ebaeSJakub Kruzik    PCDeflationSetInitOnly - Do only initialization step.
31e26ad46dSJakub Kruzik     Sets initial guess to the solution on the deflation space but does not apply
32e26ad46dSJakub Kruzik     the deflation preconditioner. The additional preconditioner is still applied.
33a122ebaeSJakub Kruzik 
34e26ad46dSJakub Kruzik    Logically Collective
35a122ebaeSJakub Kruzik 
36a122ebaeSJakub Kruzik    Input Parameters:
37a122ebaeSJakub Kruzik +  pc  - the preconditioner context
38a122ebaeSJakub Kruzik -  flg - default PETSC_FALSE
39a122ebaeSJakub Kruzik 
40e26ad46dSJakub Kruzik    Options Database Keys:
41e26ad46dSJakub Kruzik .    -pc_deflation_init_only <false> - if true computes only the special guess
42e26ad46dSJakub Kruzik 
43a122ebaeSJakub Kruzik    Level: intermediate
44a122ebaeSJakub Kruzik 
45a122ebaeSJakub Kruzik .seealso: PCDEFLATION
46a122ebaeSJakub Kruzik @*/
47a122ebaeSJakub Kruzik PetscErrorCode PCDeflationSetInitOnly(PC pc,PetscBool flg)
48a122ebaeSJakub Kruzik {
49a122ebaeSJakub Kruzik   PetscErrorCode ierr;
50a122ebaeSJakub Kruzik 
51a122ebaeSJakub Kruzik   PetscFunctionBegin;
52a122ebaeSJakub Kruzik   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
53a122ebaeSJakub Kruzik   PetscValidLogicalCollectiveBool(pc,flg,2);
54a122ebaeSJakub Kruzik   ierr = PetscTryMethod(pc,"PCDeflationSetInitOnly_C",(PC,PetscBool),(pc,flg));CHKERRQ(ierr);
55a122ebaeSJakub Kruzik   PetscFunctionReturn(0);
56a122ebaeSJakub Kruzik }
57a122ebaeSJakub Kruzik 
58a122ebaeSJakub Kruzik 
5998d6e3deSJakub Kruzik static PetscErrorCode PCDeflationSetLvl_Deflation(PC pc,PetscInt current,PetscInt max)
6098d6e3deSJakub Kruzik {
6198d6e3deSJakub Kruzik   PC_Deflation   *def = (PC_Deflation*)pc->data;
6298d6e3deSJakub Kruzik 
6398d6e3deSJakub Kruzik   PetscFunctionBegin;
646c93e71cSJakub Kruzik   if (current) def->lvl = current;
656c93e71cSJakub Kruzik   def->maxlvl = max;
6698d6e3deSJakub Kruzik   PetscFunctionReturn(0);
6798d6e3deSJakub Kruzik }
6898d6e3deSJakub Kruzik 
6998d6e3deSJakub Kruzik /*@
70e26ad46dSJakub Kruzik    PCDeflationSetMaxLvl - Set the maximum level of deflation nesting.
7198d6e3deSJakub Kruzik 
72e26ad46dSJakub Kruzik    Logically Collective
7398d6e3deSJakub Kruzik 
7498d6e3deSJakub Kruzik    Input Parameters:
7598d6e3deSJakub Kruzik +  pc  - the preconditioner context
7698d6e3deSJakub Kruzik -  max - maximum deflation level
7798d6e3deSJakub Kruzik 
78e26ad46dSJakub Kruzik    Options Database Keys:
79e26ad46dSJakub Kruzik .    -pc_deflation_max_lvl <0> - maximum number of levels for multilevel deflation
80e26ad46dSJakub Kruzik 
8198d6e3deSJakub Kruzik    Level: intermediate
8298d6e3deSJakub Kruzik 
83e26ad46dSJakub Kruzik .seealso: PCDeflationSetSpaceToCompute(), PCDeflationSetSpace(), PCDEFLATION
8498d6e3deSJakub Kruzik @*/
8598d6e3deSJakub Kruzik PetscErrorCode PCDeflationSetMaxLvl(PC pc,PetscInt max)
8698d6e3deSJakub Kruzik {
8798d6e3deSJakub Kruzik   PetscErrorCode ierr;
8898d6e3deSJakub Kruzik 
8998d6e3deSJakub Kruzik   PetscFunctionBegin;
9098d6e3deSJakub Kruzik   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
9198d6e3deSJakub Kruzik   PetscValidLogicalCollectiveInt(pc,max,2);
9298d6e3deSJakub Kruzik   ierr = PetscTryMethod(pc,"PCDeflationSetLvl_C",(PC,PetscInt,PetscInt),(pc,0,max));CHKERRQ(ierr);
9398d6e3deSJakub Kruzik   PetscFunctionReturn(0);
9498d6e3deSJakub Kruzik }
9598d6e3deSJakub Kruzik 
96859a873cSJakub Kruzik static PetscErrorCode PCDeflationSetReductionFactor_Deflation(PC pc,PetscInt red)
97859a873cSJakub Kruzik {
98859a873cSJakub Kruzik   PC_Deflation   *def = (PC_Deflation*)pc->data;
99859a873cSJakub Kruzik 
100859a873cSJakub Kruzik   PetscFunctionBegin;
101859a873cSJakub Kruzik   def->reductionfact = red;
102859a873cSJakub Kruzik   PetscFunctionReturn(0);
103859a873cSJakub Kruzik }
104859a873cSJakub Kruzik 
105859a873cSJakub Kruzik /*@
106e26ad46dSJakub Kruzik    PCDeflationSetReductionFactor - Set reduction factor for the bottom PCTELESCOPE coarse problem solver.
107859a873cSJakub Kruzik 
108e26ad46dSJakub Kruzik    Logically Collective
109859a873cSJakub Kruzik 
110859a873cSJakub Kruzik    Input Parameters:
111859a873cSJakub Kruzik +  pc  - the preconditioner context
112859a873cSJakub Kruzik -  red - reduction factor (or PETSC_DETERMINE)
113859a873cSJakub Kruzik 
114e26ad46dSJakub Kruzik    Options Database Keys:
115e26ad46dSJakub Kruzik .    -pc_deflation_reduction_factor <-1> - reduction factor on bottom level coarse problem for PCTELESCOPE
116e26ad46dSJakub Kruzik 
117e26ad46dSJakub Kruzik    Notes:
118e26ad46dSJakub Kruzik      Default is computed based on the size of the coarse problem.
119e26ad46dSJakub Kruzik 
120859a873cSJakub Kruzik    Level: intermediate
121859a873cSJakub Kruzik 
122e26ad46dSJakub Kruzik .seealso: PCTELESCOPE, PCDEFLATION
123859a873cSJakub Kruzik @*/
124859a873cSJakub Kruzik PetscErrorCode PCDeflationSetReductionFactor(PC pc,PetscInt red)
125859a873cSJakub Kruzik {
126859a873cSJakub Kruzik   PetscErrorCode ierr;
127859a873cSJakub Kruzik 
128859a873cSJakub Kruzik   PetscFunctionBegin;
129859a873cSJakub Kruzik   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
130859a873cSJakub Kruzik   PetscValidLogicalCollectiveInt(pc,red,2);
131859a873cSJakub Kruzik   ierr = PetscTryMethod(pc,"PCDeflationSetReductionFactor_C",(PC,PetscInt),(pc,red));CHKERRQ(ierr);
132859a873cSJakub Kruzik   PetscFunctionReturn(0);
133859a873cSJakub Kruzik }
134859a873cSJakub Kruzik 
1358a71cb68SJakub Kruzik static PetscErrorCode PCDeflationSetCorrectionFactor_Deflation(PC pc,PetscScalar fact)
1368a71cb68SJakub Kruzik {
1378a71cb68SJakub Kruzik   PC_Deflation   *def = (PC_Deflation*)pc->data;
1388a71cb68SJakub Kruzik 
1398a71cb68SJakub Kruzik   PetscFunctionBegin;
1408a71cb68SJakub Kruzik   /* TODO PETSC_DETERMINE -> compute max eigenvalue with power method */
1418a71cb68SJakub Kruzik   def->correct = PETSC_TRUE;
1428a71cb68SJakub Kruzik   def->correctfact = fact;
143e26ad46dSJakub Kruzik   if (def->correct == 0.0) {
1448a71cb68SJakub Kruzik     def->correct = PETSC_FALSE;
1458a71cb68SJakub Kruzik   }
1468a71cb68SJakub Kruzik   PetscFunctionReturn(0);
1478a71cb68SJakub Kruzik }
1488a71cb68SJakub Kruzik 
1498a71cb68SJakub Kruzik /*@
1508a71cb68SJakub Kruzik    PCDeflationSetCorrectionFactor - Set coarse problem correction factor.
151e26ad46dSJakub Kruzik     The Preconditioner becomes P*M^{-1} + fact*Q.
1528a71cb68SJakub Kruzik 
153e26ad46dSJakub Kruzik    Logically Collective
1548a71cb68SJakub Kruzik 
1558a71cb68SJakub Kruzik    Input Parameters:
1568a71cb68SJakub Kruzik +  pc   - the preconditioner context
157e26ad46dSJakub Kruzik -  fact - correction factor
158e26ad46dSJakub Kruzik 
159e26ad46dSJakub Kruzik    Options Database Keys:
160e26ad46dSJakub Kruzik +    -pc_deflation_correction        <false> - if true apply coarse problem correction
161e26ad46dSJakub Kruzik -    -pc_deflation_correction_factor <1.0>   - sets coarse problem correction factor
162e26ad46dSJakub Kruzik 
163e26ad46dSJakub Kruzik    Notes:
164e26ad46dSJakub Kruzik     Any non-zero fact enables the coarse problem correction.
1658a71cb68SJakub Kruzik 
1668a71cb68SJakub Kruzik    Level: intermediate
1678a71cb68SJakub Kruzik 
1688a71cb68SJakub Kruzik .seealso: PCDEFLATION
1698a71cb68SJakub Kruzik @*/
1708a71cb68SJakub Kruzik PetscErrorCode PCDeflationSetCorrectionFactor(PC pc,PetscScalar fact)
1718a71cb68SJakub Kruzik {
1728a71cb68SJakub Kruzik   PetscErrorCode ierr;
1738a71cb68SJakub Kruzik 
1748a71cb68SJakub Kruzik   PetscFunctionBegin;
1758a71cb68SJakub Kruzik   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1768a71cb68SJakub Kruzik   PetscValidLogicalCollectiveScalar(pc,fact,2);
1778a71cb68SJakub Kruzik   ierr = PetscTryMethod(pc,"PCDeflationSetCorrectionFactor_C",(PC,PetscScalar),(pc,fact));CHKERRQ(ierr);
1788a71cb68SJakub Kruzik   PetscFunctionReturn(0);
1798a71cb68SJakub Kruzik }
1808a71cb68SJakub Kruzik 
18139417d7eSJakub Kruzik static PetscErrorCode PCDeflationSetSpaceToCompute_Deflation(PC pc,PCDeflationSpaceType type,PetscInt size)
18239417d7eSJakub Kruzik {
18339417d7eSJakub Kruzik   PC_Deflation   *def = (PC_Deflation*)pc->data;
18439417d7eSJakub Kruzik 
18539417d7eSJakub Kruzik   PetscFunctionBegin;
18639417d7eSJakub Kruzik   if (type) def->spacetype = type;
18739417d7eSJakub Kruzik   if (size > 0) def->spacesize = size;
18839417d7eSJakub Kruzik   PetscFunctionReturn(0);
18939417d7eSJakub Kruzik }
19039417d7eSJakub Kruzik 
19139417d7eSJakub Kruzik /*@
19239417d7eSJakub Kruzik    PCDeflationSetSpaceToCompute - Set deflation space type and size to compute.
19339417d7eSJakub Kruzik 
194e26ad46dSJakub Kruzik    Logically Collective
19539417d7eSJakub Kruzik 
19639417d7eSJakub Kruzik    Input Parameters:
19739417d7eSJakub Kruzik +  pc   - the preconditioner context
19839417d7eSJakub Kruzik .  type - deflation space type to compute (or PETSC_IGNORE)
19939417d7eSJakub Kruzik -  size - size of the space to compute (or PETSC_DEFAULT)
20039417d7eSJakub Kruzik 
201e26ad46dSJakub Kruzik    Options Database Keys:
202e26ad46dSJakub Kruzik +    -pc_deflation_compute_space      <haar> - compute PCDeflationSpaceType deflation space
203e26ad46dSJakub Kruzik -    -pc_deflation_compute_space_size <1>    - size of the deflation space
204e26ad46dSJakub Kruzik 
20539417d7eSJakub Kruzik    Notes:
20639417d7eSJakub Kruzik     For wavelet-based deflation, size represents number of levels.
207e26ad46dSJakub Kruzik 
20839417d7eSJakub Kruzik     The deflation space is computed in PCSetUP().
20939417d7eSJakub Kruzik 
21039417d7eSJakub Kruzik    Level: intermediate
21139417d7eSJakub Kruzik 
212e26ad46dSJakub Kruzik .seealso: PCDeflationSetMaxLvl(), PCDEFLATION
21339417d7eSJakub Kruzik @*/
21439417d7eSJakub Kruzik PetscErrorCode PCDeflationSetSpaceToCompute(PC pc,PCDeflationSpaceType type,PetscInt size)
21539417d7eSJakub Kruzik {
21639417d7eSJakub Kruzik   PetscErrorCode ierr;
21739417d7eSJakub Kruzik 
21839417d7eSJakub Kruzik   PetscFunctionBegin;
21939417d7eSJakub Kruzik   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
22039417d7eSJakub Kruzik   if (type) PetscValidLogicalCollectiveEnum(pc,type,2);
22139417d7eSJakub Kruzik   if (size > 0) PetscValidLogicalCollectiveInt(pc,size,3);
22239417d7eSJakub Kruzik   ierr = PetscTryMethod(pc,"PCDeflationSetSpaceToCompute_C",(PC,PCDeflationSpaceType,PetscInt),(pc,type,size));CHKERRQ(ierr);
22339417d7eSJakub Kruzik   PetscFunctionReturn(0);
22439417d7eSJakub Kruzik }
2258513960bSJakub Kruzik 
226e662bc50SJakub Kruzik static PetscErrorCode PCDeflationSetSpace_Deflation(PC pc,Mat W,PetscBool transpose)
227e662bc50SJakub Kruzik {
228e662bc50SJakub Kruzik   PC_Deflation   *def = (PC_Deflation*)pc->data;
229e662bc50SJakub Kruzik   PetscErrorCode ierr;
230e662bc50SJakub Kruzik 
231e662bc50SJakub Kruzik   PetscFunctionBegin;
232e662bc50SJakub Kruzik   if (transpose) {
233e662bc50SJakub Kruzik     def->Wt = W;
234e662bc50SJakub Kruzik     def->W = NULL;
235e662bc50SJakub Kruzik   } else {
236e662bc50SJakub Kruzik     def->W = W;
237e662bc50SJakub Kruzik   }
238e662bc50SJakub Kruzik   ierr = PetscObjectReference((PetscObject)W);CHKERRQ(ierr);
239e662bc50SJakub Kruzik   PetscFunctionReturn(0);
240e662bc50SJakub Kruzik }
241e662bc50SJakub Kruzik 
242e662bc50SJakub Kruzik /*@
243e26ad46dSJakub Kruzik    PCDeflationSetSpace - Set the deflation space matrix (or its (hermitian) transpose).
244e662bc50SJakub Kruzik 
245e26ad46dSJakub Kruzik    Logically Collective
246e662bc50SJakub Kruzik 
247e662bc50SJakub Kruzik    Input Parameters:
248e662bc50SJakub Kruzik +  pc        - the preconditioner context
249e662bc50SJakub Kruzik .  W         - deflation matrix
250e26ad46dSJakub Kruzik -  transpose - indicates that W is an explicit transpose of the deflation matrix
251e26ad46dSJakub Kruzik 
252e26ad46dSJakub Kruzik    Notes:
253e26ad46dSJakub Kruzik     Setting W as a multipliplicative MATCOMPOSITE enables use of the multilevel
254e26ad46dSJakub Kruzik     deflation. If W = W0*W1*W2*...*Wn, W0 is taken as the first deflation space and
255e26ad46dSJakub Kruzik     the coarse problem (W0'*A*W0)^{-1} is again preconditioned by deflation with
256e26ad46dSJakub Kruzik     W1 as the deflation matrix. This repeats until the maximum level set by
257e26ad46dSJakub Kruzik     PCDeflationSetMaxLvl is reached or there are no more matrices available.
258e26ad46dSJakub Kruzik     If there are matrices left after reaching the maximum level,
259e26ad46dSJakub Kruzik     they are merged into a deflation matrix ...*W{n-1}*Wn.
260e662bc50SJakub Kruzik 
261e662bc50SJakub Kruzik    Level: intermediate
262e662bc50SJakub Kruzik 
263e26ad46dSJakub Kruzik .seealso: PCDeflationSetMaxLvl(), PCDEFLATION
264e662bc50SJakub Kruzik @*/
265e662bc50SJakub Kruzik PetscErrorCode PCDeflationSetSpace(PC pc,Mat W,PetscBool transpose)
266e662bc50SJakub Kruzik {
267e662bc50SJakub Kruzik   PetscErrorCode ierr;
268e662bc50SJakub Kruzik 
269e662bc50SJakub Kruzik   PetscFunctionBegin;
270e662bc50SJakub Kruzik   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
271e662bc50SJakub Kruzik   PetscValidHeaderSpecific(W,MAT_CLASSID,2);
272e662bc50SJakub Kruzik   PetscValidLogicalCollectiveBool(pc,transpose,3);
273e662bc50SJakub Kruzik   ierr = PetscTryMethod(pc,"PCDeflationSetSpace_C",(PC,Mat,PetscBool),(pc,W,transpose));CHKERRQ(ierr);
274e662bc50SJakub Kruzik   PetscFunctionReturn(0);
275e662bc50SJakub Kruzik }
276e662bc50SJakub Kruzik 
2773cf3a049SJakub Kruzik static PetscErrorCode PCDeflationSetProjectionNullSpaceMat_Deflation(PC pc,Mat mat)
2783cf3a049SJakub Kruzik {
2793cf3a049SJakub Kruzik   PC_Deflation     *def = (PC_Deflation*)pc->data;
2803cf3a049SJakub Kruzik   PetscErrorCode   ierr;
2813cf3a049SJakub Kruzik 
2823cf3a049SJakub Kruzik   PetscFunctionBegin;
2833cf3a049SJakub Kruzik   ierr = MatDestroy(&def->WtA);CHKERRQ(ierr);
2843cf3a049SJakub Kruzik   def->WtA = mat;
2853cf3a049SJakub Kruzik   ierr = PetscObjectReference((PetscObject)mat);CHKERRQ(ierr);
2863cf3a049SJakub Kruzik   ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)def->WtA);CHKERRQ(ierr);
2873cf3a049SJakub Kruzik   PetscFunctionReturn(0);
2883cf3a049SJakub Kruzik }
2893cf3a049SJakub Kruzik 
2903cf3a049SJakub Kruzik /*@
291e26ad46dSJakub Kruzik    PCDeflationSetProjectionNullSpaceMat - Set the projection null space matrix (W'*A).
2923cf3a049SJakub Kruzik 
293e26ad46dSJakub Kruzik    Collective
2943cf3a049SJakub Kruzik 
2953cf3a049SJakub Kruzik    Input Parameters:
2963cf3a049SJakub Kruzik +  pc  - preconditioner context
2973cf3a049SJakub Kruzik -  mat - projection null space matrix
2983cf3a049SJakub Kruzik 
2993cf3a049SJakub Kruzik    Level: developer
3003cf3a049SJakub Kruzik 
3013cf3a049SJakub Kruzik .seealso: PCDEFLATION
3023cf3a049SJakub Kruzik @*/
3033cf3a049SJakub Kruzik PetscErrorCode  PCDeflationSetProjectionNullSpaceMat(PC pc,Mat mat)
3043cf3a049SJakub Kruzik {
3053cf3a049SJakub Kruzik   PetscErrorCode ierr;
3063cf3a049SJakub Kruzik 
3073cf3a049SJakub Kruzik   PetscFunctionBegin;
3083cf3a049SJakub Kruzik   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
3093cf3a049SJakub Kruzik   PetscValidHeaderSpecific(mat,MAT_CLASSID,2);
3103cf3a049SJakub Kruzik   ierr = PetscTryMethod(pc,"PCDeflationSetProjectionNullSpaceMat_C",(PC,Mat),(pc,mat));CHKERRQ(ierr);
3113cf3a049SJakub Kruzik   PetscFunctionReturn(0);
3123cf3a049SJakub Kruzik }
3133cf3a049SJakub Kruzik 
314e924b002SJakub Kruzik static PetscErrorCode PCDeflationSetCoarseMat_Deflation(PC pc,Mat mat)
315e924b002SJakub Kruzik {
316e924b002SJakub Kruzik   PC_Deflation     *def = (PC_Deflation*)pc->data;
317e924b002SJakub Kruzik   PetscErrorCode   ierr;
318e924b002SJakub Kruzik 
319e924b002SJakub Kruzik   PetscFunctionBegin;
320e924b002SJakub Kruzik   ierr = MatDestroy(&def->WtAW);CHKERRQ(ierr);
321e924b002SJakub Kruzik   def->WtAW = mat;
322e924b002SJakub Kruzik   ierr = PetscObjectReference((PetscObject)mat);CHKERRQ(ierr);
323e924b002SJakub Kruzik   ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)def->WtAW);CHKERRQ(ierr);
324e924b002SJakub Kruzik   PetscFunctionReturn(0);
325e924b002SJakub Kruzik }
326e924b002SJakub Kruzik 
327e924b002SJakub Kruzik /*@
328e26ad46dSJakub Kruzik    PCDeflationSetCoarseMat - Set the coarse problem Mat.
329e924b002SJakub Kruzik 
330e26ad46dSJakub Kruzik    Collective
331e924b002SJakub Kruzik 
332e924b002SJakub Kruzik    Input Parameters:
333e924b002SJakub Kruzik +  pc  - preconditioner context
334e924b002SJakub Kruzik -  mat - coarse problem mat
335e924b002SJakub Kruzik 
336e924b002SJakub Kruzik    Level: developer
337e924b002SJakub Kruzik 
338e924b002SJakub Kruzik .seealso: PCDEFLATION
339e924b002SJakub Kruzik @*/
340e924b002SJakub Kruzik PetscErrorCode  PCDeflationSetCoarseMat(PC pc,Mat mat)
341e924b002SJakub Kruzik {
342e924b002SJakub Kruzik   PetscErrorCode ierr;
343e924b002SJakub Kruzik 
344e924b002SJakub Kruzik   PetscFunctionBegin;
345e924b002SJakub Kruzik   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
346e924b002SJakub Kruzik   PetscValidHeaderSpecific(mat,MAT_CLASSID,2);
347e924b002SJakub Kruzik   ierr = PetscTryMethod(pc,"PCDeflationSetCoarseMat_C",(PC,Mat),(pc,mat));CHKERRQ(ierr);
348e924b002SJakub Kruzik   PetscFunctionReturn(0);
349e924b002SJakub Kruzik }
350e924b002SJakub Kruzik 
35198d6e3deSJakub Kruzik static PetscErrorCode PCDeflationGetCoarseKSP_Deflation(PC pc,KSP *ksp)
3525c0c31c5SJakub Kruzik {
3535c0c31c5SJakub Kruzik   PC_Deflation     *def = (PC_Deflation*)pc->data;
3545c0c31c5SJakub Kruzik 
3555c0c31c5SJakub Kruzik   PetscFunctionBegin;
35698d6e3deSJakub Kruzik   *ksp = def->WtAWinv;
3575c0c31c5SJakub Kruzik   PetscFunctionReturn(0);
3585c0c31c5SJakub Kruzik }
3595c0c31c5SJakub Kruzik 
3605c0c31c5SJakub Kruzik /*@
36198d6e3deSJakub Kruzik    PCDeflationGetCoarseKSP - Returns a pointer to the coarse problem KSP.
3625c0c31c5SJakub Kruzik 
36398d6e3deSJakub Kruzik    Not Collective
3645c0c31c5SJakub Kruzik 
3655c0c31c5SJakub Kruzik    Input Parameters:
36698d6e3deSJakub Kruzik .  pc - preconditioner context
3675c0c31c5SJakub Kruzik 
368e26ad46dSJakub Kruzik    Output Parameters:
36998d6e3deSJakub Kruzik .  ksp - coarse problem KSP context
3705c0c31c5SJakub Kruzik 
371e26ad46dSJakub Kruzik    Level: advanced
37298d6e3deSJakub Kruzik 
37398d6e3deSJakub Kruzik .seealso: PCDeflationSetCoarseKSP(), PCDEFLATION
3745c0c31c5SJakub Kruzik @*/
37598d6e3deSJakub Kruzik PetscErrorCode  PCDeflationGetCoarseKSP(PC pc,KSP *ksp)
3765c0c31c5SJakub Kruzik {
3775c0c31c5SJakub Kruzik   PetscErrorCode ierr;
3785c0c31c5SJakub Kruzik 
3795c0c31c5SJakub Kruzik   PetscFunctionBegin;
38022b0793eSJakub Kruzik   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
38198d6e3deSJakub Kruzik   PetscValidPointer(ksp,2);
38298d6e3deSJakub Kruzik   ierr = PetscTryMethod(pc,"PCDeflationGetCoarseKSP_C",(PC,KSP*),(pc,ksp));CHKERRQ(ierr);
38398d6e3deSJakub Kruzik   PetscFunctionReturn(0);
38498d6e3deSJakub Kruzik }
38598d6e3deSJakub Kruzik 
38698d6e3deSJakub Kruzik static PetscErrorCode PCDeflationSetCoarseKSP_Deflation(PC pc,KSP ksp)
38798d6e3deSJakub Kruzik {
38898d6e3deSJakub Kruzik   PC_Deflation     *def = (PC_Deflation*)pc->data;
38998d6e3deSJakub Kruzik   PetscErrorCode   ierr;
39098d6e3deSJakub Kruzik 
39198d6e3deSJakub Kruzik   PetscFunctionBegin;
39298d6e3deSJakub Kruzik   ierr = KSPDestroy(&def->WtAWinv);CHKERRQ(ierr);
39398d6e3deSJakub Kruzik   def->WtAWinv = ksp;
39498d6e3deSJakub Kruzik   ierr = PetscObjectReference((PetscObject)ksp);CHKERRQ(ierr);
39598d6e3deSJakub Kruzik   ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)def->WtAWinv);CHKERRQ(ierr);
39698d6e3deSJakub Kruzik   PetscFunctionReturn(0);
39798d6e3deSJakub Kruzik }
39898d6e3deSJakub Kruzik 
39998d6e3deSJakub Kruzik /*@
400e26ad46dSJakub Kruzik    PCDeflationSetCoarseKSP - Set the coarse problem KSP.
40198d6e3deSJakub Kruzik 
402e26ad46dSJakub Kruzik    Collective
40398d6e3deSJakub Kruzik 
40498d6e3deSJakub Kruzik    Input Parameters:
40598d6e3deSJakub Kruzik +  pc  - preconditioner context
40698d6e3deSJakub Kruzik -  ksp - coarse problem KSP context
40798d6e3deSJakub Kruzik 
40898d6e3deSJakub Kruzik    Level: developer
40998d6e3deSJakub Kruzik 
41098d6e3deSJakub Kruzik .seealso: PCDeflationGetCoarseKSP(), PCDEFLATION
41198d6e3deSJakub Kruzik @*/
41298d6e3deSJakub Kruzik PetscErrorCode  PCDeflationSetCoarseKSP(PC pc,KSP ksp)
41398d6e3deSJakub Kruzik {
41498d6e3deSJakub Kruzik   PetscErrorCode ierr;
41598d6e3deSJakub Kruzik 
41698d6e3deSJakub Kruzik   PetscFunctionBegin;
41798d6e3deSJakub Kruzik   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
41898d6e3deSJakub Kruzik   PetscValidHeaderSpecific(ksp,KSP_CLASSID,2);
41998d6e3deSJakub Kruzik   PetscCheckSameComm(pc,1,ksp,2);
42098d6e3deSJakub Kruzik   ierr = PetscTryMethod(pc,"PCDeflationSetCoarseKSP_C",(PC,KSP),(pc,ksp));CHKERRQ(ierr);
4215c0c31c5SJakub Kruzik   PetscFunctionReturn(0);
4225c0c31c5SJakub Kruzik }
423e662bc50SJakub Kruzik 
424268c6673SJakub Kruzik static PetscErrorCode PCDeflationGetPC_Deflation(PC pc,PC *apc)
425268c6673SJakub Kruzik {
426268c6673SJakub Kruzik   PC_Deflation   *def = (PC_Deflation*)pc->data;
427268c6673SJakub Kruzik 
428268c6673SJakub Kruzik   PetscFunctionBegin;
429268c6673SJakub Kruzik   *apc = def->pc;
430268c6673SJakub Kruzik   PetscFunctionReturn(0);
431268c6673SJakub Kruzik }
432268c6673SJakub Kruzik 
433268c6673SJakub Kruzik /*@
434e26ad46dSJakub Kruzik    PCDeflationGetPC - Returns a pointer to the additional preconditioner.
435268c6673SJakub Kruzik 
436268c6673SJakub Kruzik    Not Collective
437268c6673SJakub Kruzik 
438268c6673SJakub Kruzik    Input Parameters:
439268c6673SJakub Kruzik .  pc  - the preconditioner context
440268c6673SJakub Kruzik 
441e26ad46dSJakub Kruzik    Output Parameters:
442268c6673SJakub Kruzik .  apc - additional preconditioner
443268c6673SJakub Kruzik 
444268c6673SJakub Kruzik    Level: advanced
445268c6673SJakub Kruzik 
446268c6673SJakub Kruzik .seealso: PCDeflationSetPC(), PCDEFLATION
447268c6673SJakub Kruzik @*/
448268c6673SJakub Kruzik PetscErrorCode PCDeflationGetPC(PC pc,PC *apc)
449268c6673SJakub Kruzik {
450268c6673SJakub Kruzik   PetscErrorCode ierr;
451268c6673SJakub Kruzik 
452268c6673SJakub Kruzik   PetscFunctionBegin;
453268c6673SJakub Kruzik   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
454268c6673SJakub Kruzik   PetscValidPointer(pc,2);
455268c6673SJakub Kruzik   ierr = PetscTryMethod(pc,"PCDeflationGetPC_C",(PC,PC*),(pc,apc));CHKERRQ(ierr);
456268c6673SJakub Kruzik   PetscFunctionReturn(0);
457268c6673SJakub Kruzik }
458268c6673SJakub Kruzik 
45922b0793eSJakub Kruzik static PetscErrorCode PCDeflationSetPC_Deflation(PC pc,PC apc)
46022b0793eSJakub Kruzik {
46122b0793eSJakub Kruzik   PC_Deflation   *def = (PC_Deflation*)pc->data;
46222b0793eSJakub Kruzik   PetscErrorCode ierr;
46322b0793eSJakub Kruzik 
46422b0793eSJakub Kruzik   PetscFunctionBegin;
46522b0793eSJakub Kruzik   ierr = PCDestroy(&def->pc);CHKERRQ(ierr);
46622b0793eSJakub Kruzik   def->pc = apc;
46722b0793eSJakub Kruzik   ierr = PetscObjectReference((PetscObject)apc);CHKERRQ(ierr);
46822b0793eSJakub Kruzik   ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)def->pc);CHKERRQ(ierr);
46922b0793eSJakub Kruzik   PetscFunctionReturn(0);
47022b0793eSJakub Kruzik }
47122b0793eSJakub Kruzik 
47222b0793eSJakub Kruzik /*@
473e26ad46dSJakub Kruzik    PCDeflationSetPC - Set the additional preconditioner.
47422b0793eSJakub Kruzik 
475e26ad46dSJakub Kruzik    Collective
47622b0793eSJakub Kruzik 
47722b0793eSJakub Kruzik    Input Parameters:
47822b0793eSJakub Kruzik +  pc  - the preconditioner context
47922b0793eSJakub Kruzik -  apc - additional preconditioner
48022b0793eSJakub Kruzik 
481268c6673SJakub Kruzik    Level: developer
48222b0793eSJakub Kruzik 
483268c6673SJakub Kruzik .seealso: PCDeflationGetPC(), PCDEFLATION
48422b0793eSJakub Kruzik @*/
48522b0793eSJakub Kruzik PetscErrorCode PCDeflationSetPC(PC pc,PC apc)
48622b0793eSJakub Kruzik {
48722b0793eSJakub Kruzik   PetscErrorCode ierr;
48822b0793eSJakub Kruzik 
48922b0793eSJakub Kruzik   PetscFunctionBegin;
49022b0793eSJakub Kruzik   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
49122b0793eSJakub Kruzik   PetscValidHeaderSpecific(apc,PC_CLASSID,2);
49222b0793eSJakub Kruzik   PetscCheckSameComm(pc,1,apc,2);
49322b0793eSJakub Kruzik   ierr = PetscTryMethod(pc,"PCDeflationSetPC_C",(PC,PC),(pc,apc));CHKERRQ(ierr);
49422b0793eSJakub Kruzik   PetscFunctionReturn(0);
49522b0793eSJakub Kruzik }
49622b0793eSJakub Kruzik 
49737eeb815SJakub Kruzik /*
498b27c8092SJakub Kruzik   x <- x + W*(W'*A*W)^{-1}*W'*r  = x + Q*r
499b27c8092SJakub Kruzik */
500b27c8092SJakub Kruzik static PetscErrorCode PCPreSolve_Deflation(PC pc,KSP ksp,Vec b, Vec x)
501b27c8092SJakub Kruzik {
502b27c8092SJakub Kruzik   PC_Deflation     *def = (PC_Deflation*)pc->data;
503b27c8092SJakub Kruzik   Mat              A;
504b27c8092SJakub Kruzik   Vec              r,w1,w2;
505b27c8092SJakub Kruzik   PetscBool        nonzero;
506b27c8092SJakub Kruzik   PetscErrorCode   ierr;
50737eeb815SJakub Kruzik 
508b27c8092SJakub Kruzik   PetscFunctionBegin;
509b27c8092SJakub Kruzik   w1 = def->workcoarse[0];
510b27c8092SJakub Kruzik   w2 = def->workcoarse[1];
511b27c8092SJakub Kruzik   r  = def->work;
512b27c8092SJakub Kruzik   ierr = PCGetOperators(pc,NULL,&A);CHKERRQ(ierr);
51337eeb815SJakub Kruzik 
514b27c8092SJakub Kruzik   ierr = KSPGetInitialGuessNonzero(ksp,&nonzero);CHKERRQ(ierr);
515b27c8092SJakub Kruzik   ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr);
516b27c8092SJakub Kruzik   if (nonzero) {
517b27c8092SJakub Kruzik     ierr = MatMult(A,x,r);CHKERRQ(ierr);                          /*    r  <- b - Ax              */
518b27c8092SJakub Kruzik     ierr = VecAYPX(r,-1.0,b);CHKERRQ(ierr);
519b27c8092SJakub Kruzik   } else {
520b27c8092SJakub Kruzik     ierr = VecCopy(b,r);CHKERRQ(ierr);                            /*    r  <- b (x is 0)          */
521b27c8092SJakub Kruzik   }
522b27c8092SJakub Kruzik 
523a3177931SJakub Kruzik   if (def->Wt) {
524a3177931SJakub Kruzik     ierr = MatMult(def->Wt,r,w1);CHKERRQ(ierr);                   /*    w1 <- W'*r                */
525a3177931SJakub Kruzik   } else {
526a3177931SJakub Kruzik     ierr = MatMultHermitianTranspose(def->W,r,w1);CHKERRQ(ierr);  /*    w1 <- W'*r                */
527a3177931SJakub Kruzik   }
528b27c8092SJakub Kruzik   ierr = KSPSolve(def->WtAWinv,w1,w2);CHKERRQ(ierr);              /*    w2 <- (W'*A*W)^{-1}*w1    */
529b27c8092SJakub Kruzik   ierr = MatMult(def->W,w2,r);CHKERRQ(ierr);                      /*    r  <- W*w2                */
530b27c8092SJakub Kruzik   ierr = VecAYPX(x,1.0,r);CHKERRQ(ierr);
531b27c8092SJakub Kruzik   PetscFunctionReturn(0);
532b27c8092SJakub Kruzik }
53337eeb815SJakub Kruzik 
534f8bf229cSJakub Kruzik /*
535f8bf229cSJakub Kruzik   if (def->correct) {
536ae029463SJakub Kruzik     z <- M^{-1}r - W*(W'*A*W)^{-1}*(W'*A*M^{-1}r - l*W'*r) = (P*M^{-1} + l*Q)*r
537f8bf229cSJakub Kruzik   } else {
538ae029463SJakub Kruzik     z <- M^{-1}*r - W*(W'*A*W)^{-1}*W'*A*M{-1}*r = P*M^{-1}*r
539f8bf229cSJakub Kruzik   }
54037eeb815SJakub Kruzik */
541f8bf229cSJakub Kruzik static PetscErrorCode PCApply_Deflation(PC pc,Vec r,Vec z)
542f8bf229cSJakub Kruzik {
543f8bf229cSJakub Kruzik   PC_Deflation     *def = (PC_Deflation*)pc->data;
544f8bf229cSJakub Kruzik   Mat              A;
545f8bf229cSJakub Kruzik   Vec              u,w1,w2;
546f8bf229cSJakub Kruzik   PetscErrorCode   ierr;
547f8bf229cSJakub Kruzik 
548f8bf229cSJakub Kruzik   PetscFunctionBegin;
549f8bf229cSJakub Kruzik   w1 = def->workcoarse[0];
550f8bf229cSJakub Kruzik   w2 = def->workcoarse[1];
551f8bf229cSJakub Kruzik   u  = def->work;
552f8bf229cSJakub Kruzik   ierr = PCGetOperators(pc,NULL,&A);CHKERRQ(ierr);
553f8bf229cSJakub Kruzik 
554ae029463SJakub Kruzik   ierr = PCApply(def->pc,r,z);CHKERRQ(ierr);                          /*    z <- M^{-1}*r             */
55522b0793eSJakub Kruzik   if (!def->init) {
556ae029463SJakub Kruzik     ierr = MatMult(def->WtA,z,w1);CHKERRQ(ierr);                      /*    w1 <- W'*A*z              */
557f8bf229cSJakub Kruzik     if (def->correct) {
558ae029463SJakub Kruzik       if (def->Wt) {
559ae029463SJakub Kruzik         ierr = MatMult(def->Wt,r,w2);CHKERRQ(ierr);                   /*    w2 <- W'*r                */
560ae029463SJakub Kruzik       } else {
561a3177931SJakub Kruzik         ierr = MatMultHermitianTranspose(def->W,r,w2);CHKERRQ(ierr);  /*    w2 <- W'*r                */
562f8bf229cSJakub Kruzik       }
563ae029463SJakub Kruzik       ierr = VecAXPY(w1,-1.0*def->correctfact,w2);CHKERRQ(ierr);      /*    w1 <- w1 - l*w2           */
564f8bf229cSJakub Kruzik     }
565f8bf229cSJakub Kruzik     ierr = KSPSolve(def->WtAWinv,w1,w2);CHKERRQ(ierr);                /*    w2 <- (W'*A*W)^{-1}*w1    */
566f8bf229cSJakub Kruzik     ierr = MatMult(def->W,w2,u);CHKERRQ(ierr);                        /*    u  <- W*w2                */
56722b0793eSJakub Kruzik     ierr = VecAXPY(z,-1.0,u);CHKERRQ(ierr);                           /*    z  <- z - u               */
56822b0793eSJakub Kruzik   }
569f8bf229cSJakub Kruzik   PetscFunctionReturn(0);
570f8bf229cSJakub Kruzik }
571f8bf229cSJakub Kruzik 
57237eeb815SJakub Kruzik static PetscErrorCode PCSetUp_Deflation(PC pc)
57337eeb815SJakub Kruzik {
574409a357bSJakub Kruzik   PC_Deflation     *def = (PC_Deflation*)pc->data;
575409a357bSJakub Kruzik   KSP              innerksp;
576409a357bSJakub Kruzik   PC               pcinner;
577409a357bSJakub Kruzik   Mat              Amat,nextDef=NULL,*mats;
578409a357bSJakub Kruzik   PetscInt         i,m,red,size,commsize;
579409a357bSJakub Kruzik   PetscBool        match,flgspd,transp=PETSC_FALSE;
580409a357bSJakub Kruzik   MatCompositeType ctype;
581409a357bSJakub Kruzik   MPI_Comm         comm;
5826c93e71cSJakub Kruzik   char             prefix[128]="";
58337eeb815SJakub Kruzik   PetscErrorCode   ierr;
58437eeb815SJakub Kruzik 
58537eeb815SJakub Kruzik   PetscFunctionBegin;
586409a357bSJakub Kruzik   if (pc->setupcalled) PetscFunctionReturn(0);
587409a357bSJakub Kruzik   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
588f8bf229cSJakub Kruzik   ierr = PCGetOperators(pc,NULL,&Amat);CHKERRQ(ierr);
5896c93e71cSJakub Kruzik   if (!def->lvl && !def->prefix) {
5906c93e71cSJakub Kruzik     ierr = PCGetOptionsPrefix(pc,&def->prefix);CHKERRQ(ierr);
5916c93e71cSJakub Kruzik   }
5926c93e71cSJakub Kruzik   if (def->lvl) {
5936c93e71cSJakub Kruzik     sprintf(prefix,"%d_",(int)def->lvl);
5946c93e71cSJakub Kruzik   }
595f8bf229cSJakub Kruzik 
596f8bf229cSJakub Kruzik   /* compute a deflation space */
597409a357bSJakub Kruzik   if (def->W || def->Wt) {
598409a357bSJakub Kruzik     def->spacetype = PC_DEFLATION_SPACE_USER;
599409a357bSJakub Kruzik   } else {
600e53e0a0dSJakub Kruzik     ierr = PCDeflationComputeSpace(pc);CHKERRQ(ierr);
60137eeb815SJakub Kruzik   }
60237eeb815SJakub Kruzik 
603409a357bSJakub Kruzik   /* nested deflation */
604409a357bSJakub Kruzik   if (def->W) {
605409a357bSJakub Kruzik     ierr = PetscObjectTypeCompare((PetscObject)def->W,MATCOMPOSITE,&match);CHKERRQ(ierr);
606409a357bSJakub Kruzik     if (match) {
607409a357bSJakub Kruzik       ierr = MatCompositeGetType(def->W,&ctype);CHKERRQ(ierr);
608409a357bSJakub Kruzik       ierr = MatCompositeGetNumberMat(def->W,&size);CHKERRQ(ierr);
60937eeb815SJakub Kruzik     }
610409a357bSJakub Kruzik   } else {
611a3177931SJakub Kruzik     ierr = MatCreateHermitianTranspose(def->Wt,&def->W);CHKERRQ(ierr);
612409a357bSJakub Kruzik     ierr = PetscObjectTypeCompare((PetscObject)def->Wt,MATCOMPOSITE,&match);CHKERRQ(ierr);
613409a357bSJakub Kruzik     if (match) {
614409a357bSJakub Kruzik       ierr = MatCompositeGetType(def->Wt,&ctype);CHKERRQ(ierr);
615409a357bSJakub Kruzik       ierr = MatCompositeGetNumberMat(def->Wt,&size);CHKERRQ(ierr);
616409a357bSJakub Kruzik     }
617409a357bSJakub Kruzik     transp = PETSC_TRUE;
618409a357bSJakub Kruzik   }
619409a357bSJakub Kruzik   if (match && ctype == MAT_COMPOSITE_MULTIPLICATIVE) {
620409a357bSJakub Kruzik     if (!transp) {
6216c93e71cSJakub Kruzik       if (def->lvl < def->maxlvl) {
62228da0170SJakub Kruzik         ierr = PetscMalloc1(size,&mats);CHKERRQ(ierr);
623409a357bSJakub Kruzik         for (i=0; i<size; i++) {
624409a357bSJakub Kruzik           ierr = MatCompositeGetMat(def->W,i,&mats[i]);CHKERRQ(ierr);
625409a357bSJakub Kruzik         }
626409a357bSJakub Kruzik         size -= 1;
627409a357bSJakub Kruzik         ierr = MatDestroy(&def->W);CHKERRQ(ierr);
628409a357bSJakub Kruzik         def->W = mats[size];
629409a357bSJakub Kruzik         ierr = PetscObjectReference((PetscObject)mats[size]);CHKERRQ(ierr);
630409a357bSJakub Kruzik         if (size > 1) {
631409a357bSJakub Kruzik           ierr = MatCreateComposite(comm,size,mats,&nextDef);CHKERRQ(ierr);
632409a357bSJakub Kruzik           ierr = MatCompositeSetType(nextDef,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr);
633409a357bSJakub Kruzik         } else {
634409a357bSJakub Kruzik           nextDef = mats[0];
635409a357bSJakub Kruzik           ierr = PetscObjectReference((PetscObject)mats[0]);CHKERRQ(ierr);
636409a357bSJakub Kruzik         }
63728da0170SJakub Kruzik         ierr = PetscFree(mats);CHKERRQ(ierr);
638409a357bSJakub Kruzik       } else {
639409a357bSJakub Kruzik         /* ierr = MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT);CHKERRQ(ierr); */
640409a357bSJakub Kruzik         ierr = MatCompositeMerge(def->W);CHKERRQ(ierr);
641409a357bSJakub Kruzik       }
64228da0170SJakub Kruzik     } else {
6436c93e71cSJakub Kruzik       if (def->lvl < def->maxlvl) {
64428da0170SJakub Kruzik         ierr = PetscMalloc1(size,&mats);CHKERRQ(ierr);
64528da0170SJakub Kruzik         for (i=0; i<size; i++) {
64628da0170SJakub Kruzik           ierr = MatCompositeGetMat(def->Wt,i,&mats[i]);CHKERRQ(ierr);
64728da0170SJakub Kruzik         }
64828da0170SJakub Kruzik         size -= 1;
64928da0170SJakub Kruzik         ierr = MatDestroy(&def->Wt);CHKERRQ(ierr);
65028da0170SJakub Kruzik         def->Wt = mats[0];
65128da0170SJakub Kruzik         ierr = PetscObjectReference((PetscObject)mats[0]);CHKERRQ(ierr);
65228da0170SJakub Kruzik         if (size > 1) {
65328da0170SJakub Kruzik           ierr = MatCreateComposite(comm,size,&mats[1],&nextDef);CHKERRQ(ierr);
65428da0170SJakub Kruzik           ierr = MatCompositeSetType(nextDef,MAT_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr);
65528da0170SJakub Kruzik         } else {
65628da0170SJakub Kruzik           nextDef = mats[1];
65728da0170SJakub Kruzik           ierr = PetscObjectReference((PetscObject)mats[1]);CHKERRQ(ierr);
658409a357bSJakub Kruzik         }
659409a357bSJakub Kruzik         ierr = PetscFree(mats);CHKERRQ(ierr);
66028da0170SJakub Kruzik       } else {
66128da0170SJakub Kruzik         /* ierr = MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT);CHKERRQ(ierr); */
66228da0170SJakub Kruzik         ierr = MatCompositeMerge(def->Wt);CHKERRQ(ierr);
66328da0170SJakub Kruzik       }
66428da0170SJakub Kruzik     }
66528da0170SJakub Kruzik   }
66628da0170SJakub Kruzik 
66728da0170SJakub Kruzik   if (transp) {
66828da0170SJakub Kruzik     ierr = MatDestroy(&def->W);CHKERRQ(ierr);
669a3177931SJakub Kruzik     ierr = MatHermitianTranspose(def->Wt,MAT_INITIAL_MATRIX,&def->W);CHKERRQ(ierr);
670409a357bSJakub Kruzik   }
671409a357bSJakub Kruzik 
672ae029463SJakub Kruzik   /* assemble WtA */
673ae029463SJakub Kruzik   if (!def->WtA) {
674ae029463SJakub Kruzik     if (def->Wt) {
675ae029463SJakub Kruzik       ierr = MatMatMult(def->Wt,Amat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtA);CHKERRQ(ierr);
676ae029463SJakub Kruzik     } else {
677a3177931SJakub Kruzik #if defined(PETSC_USE_COMPLEX)
678a3177931SJakub Kruzik       ierr = MatHermitianTranspose(def->W,MAT_INITIAL_MATRIX,&def->Wt);CHKERRQ(ierr);
679*9e56ec8aSJakub Kruzik       ierr = MatMatMult(def->Wt,Amat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtA);CHKERRQ(ierr);
680a3177931SJakub Kruzik #else
681ae029463SJakub Kruzik       ierr = MatTransposeMatMult(def->W,Amat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtA);CHKERRQ(ierr);
682a3177931SJakub Kruzik #endif
683ae029463SJakub Kruzik     }
684ae029463SJakub Kruzik   }
685409a357bSJakub Kruzik   /* setup coarse problem */
686409a357bSJakub Kruzik   if (!def->WtAWinv) {
68728da0170SJakub Kruzik     ierr = MatGetSize(def->W,NULL,&m);CHKERRQ(ierr);
688409a357bSJakub Kruzik     if (!def->WtAW) {
689ae029463SJakub Kruzik       ierr = MatMatMult(def->WtA,def->W,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtAW);CHKERRQ(ierr);
690409a357bSJakub Kruzik       /* TODO create MatInheritOption(Mat,MatOption) */
691409a357bSJakub Kruzik       ierr = MatGetOption(Amat,MAT_SPD,&flgspd);CHKERRQ(ierr);
692409a357bSJakub Kruzik       ierr = MatSetOption(def->WtAW,MAT_SPD,flgspd);CHKERRQ(ierr);
693409a357bSJakub Kruzik #if defined(PETSC_USE_DEBUG)
694ae029463SJakub Kruzik       /* Check columns of W are not in kernel of A */
695409a357bSJakub Kruzik       PetscReal *norms;
696409a357bSJakub Kruzik       ierr = PetscMalloc1(m,&norms);CHKERRQ(ierr);
697409a357bSJakub Kruzik       ierr = MatGetColumnNorms(def->WtAW,NORM_INFINITY,norms);CHKERRQ(ierr);
698409a357bSJakub Kruzik       for (i=0; i<m; i++) {
699409a357bSJakub Kruzik         if (norms[i] < 100*PETSC_MACHINE_EPSILON) {
700409a357bSJakub Kruzik           SETERRQ1(comm,PETSC_ERR_SUP,"Column %D of W is in kernel of A.",i);
701409a357bSJakub Kruzik         }
702409a357bSJakub Kruzik       }
703409a357bSJakub Kruzik       ierr = PetscFree(norms);CHKERRQ(ierr);
704409a357bSJakub Kruzik #endif
705409a357bSJakub Kruzik     } else {
706409a357bSJakub Kruzik       ierr = MatGetOption(def->WtAW,MAT_SPD,&flgspd);CHKERRQ(ierr);
707409a357bSJakub Kruzik     }
708409a357bSJakub Kruzik     /* TODO use MATINV */
709409a357bSJakub Kruzik     ierr = KSPCreate(comm,&def->WtAWinv);CHKERRQ(ierr);
710409a357bSJakub Kruzik     ierr = KSPSetOperators(def->WtAWinv,def->WtAW,def->WtAW);CHKERRQ(ierr);
711409a357bSJakub Kruzik     ierr = KSPGetPC(def->WtAWinv,&pcinner);CHKERRQ(ierr);
712557a2f7dSJakub Kruzik     /* Setup KSP and PC */
713557a2f7dSJakub Kruzik     if (nextDef) { /* next level for multilevel deflation */
714557a2f7dSJakub Kruzik       innerksp = def->WtAWinv;
715557a2f7dSJakub Kruzik       /* set default KSPtype */
716557a2f7dSJakub Kruzik       if (!def->ksptype) {
717557a2f7dSJakub Kruzik         def->ksptype = KSPFGMRES;
718557a2f7dSJakub Kruzik         if (flgspd) { /* SPD system */
719557a2f7dSJakub Kruzik           def->ksptype = KSPFCG;
720557a2f7dSJakub Kruzik         }
721557a2f7dSJakub Kruzik       }
722ae029463SJakub Kruzik       ierr = KSPSetType(innerksp,def->ksptype);CHKERRQ(ierr); /* TODO iherit from KSP + tolerances */
723557a2f7dSJakub Kruzik       ierr = PCSetType(pcinner,PCDEFLATION);CHKERRQ(ierr); /* TODO create coarse preconditinoner M_c = WtMW ? */
724557a2f7dSJakub Kruzik       ierr = PCDeflationSetSpace(pcinner,nextDef,transp);CHKERRQ(ierr);
7256c93e71cSJakub Kruzik       ierr = PCDeflationSetLvl_Deflation(pcinner,def->lvl+1,def->maxlvl);CHKERRQ(ierr);
726ae029463SJakub Kruzik       /* inherit options */
7276c93e71cSJakub Kruzik       if (def->prefix) ((PC_Deflation*)(pcinner->data))->prefix = def->prefix;
7289f604af8SJakub Kruzik       ((PC_Deflation*)(pcinner->data))->init          = def->init;
729557a2f7dSJakub Kruzik       ((PC_Deflation*)(pcinner->data))->ksptype       = def->ksptype;
730557a2f7dSJakub Kruzik       ((PC_Deflation*)(pcinner->data))->correct       = def->correct;
7319f604af8SJakub Kruzik       ((PC_Deflation*)(pcinner->data))->correctfact   = def->correctfact;
7329f604af8SJakub Kruzik       ((PC_Deflation*)(pcinner->data))->reductionfact = def->reductionfact;
733557a2f7dSJakub Kruzik       ierr = MatDestroy(&nextDef);CHKERRQ(ierr);
734557a2f7dSJakub Kruzik     } else { /* the last level */
735557a2f7dSJakub Kruzik       ierr = KSPSetType(def->WtAWinv,KSPPREONLY);CHKERRQ(ierr);
736409a357bSJakub Kruzik       ierr = PCSetType(pcinner,PCTELESCOPE);CHKERRQ(ierr);
7376c93e71cSJakub Kruzik       /* do not overwrite PCTELESCOPE */
7386c93e71cSJakub Kruzik       if (def->prefix) {
7396c93e71cSJakub Kruzik         ierr = KSPSetOptionsPrefix(def->WtAWinv,def->prefix);CHKERRQ(ierr);
740ac66f006SJakub Kruzik       }
7416c93e71cSJakub Kruzik       ierr = KSPAppendOptionsPrefix(def->WtAWinv,"def_tel_");CHKERRQ(ierr);
742ac66f006SJakub Kruzik       ierr = PCSetFromOptions(pcinner);CHKERRQ(ierr);
743409a357bSJakub Kruzik       /* Reduction factor choice */
744409a357bSJakub Kruzik       red = def->reductionfact;
745409a357bSJakub Kruzik       if (red < 0) {
746409a357bSJakub Kruzik         ierr = MPI_Comm_size(comm,&commsize);CHKERRQ(ierr);
747409a357bSJakub Kruzik         red  = ceil((float)commsize/ceil((float)m/commsize));
748409a357bSJakub Kruzik         ierr = PetscObjectTypeCompareAny((PetscObject)(def->WtAW),&match,MATSEQDENSE,MATMPIDENSE,MATDENSE,"");CHKERRQ(ierr);
749409a357bSJakub Kruzik         if (match) red = commsize;
7506c93e71cSJakub Kruzik         ierr = PetscInfo1(pc,"Auto choosing reduction factor %D\n",red);CHKERRQ(ierr);
751409a357bSJakub Kruzik       }
752409a357bSJakub Kruzik       ierr = PCTelescopeSetReductionFactor(pcinner,red);CHKERRQ(ierr);
753ac66f006SJakub Kruzik       ierr = PCSetUp(pcinner);CHKERRQ(ierr);
754409a357bSJakub Kruzik       ierr = PCTelescopeGetKSP(pcinner,&innerksp);CHKERRQ(ierr);
755ac66f006SJakub Kruzik       if (innerksp) {
756409a357bSJakub Kruzik         ierr = KSPGetPC(innerksp,&pcinner);CHKERRQ(ierr);
757ac66f006SJakub Kruzik         ierr = PCSetType(pcinner,PCLU);CHKERRQ(ierr);
758481b7641SJakub Kruzik #if defined(PETSC_HAVE_SUPERLU)
759481b7641SJakub Kruzik         ierr = MatGetFactorAvailable(def->WtAW,MATSOLVERSUPERLU,MAT_FACTOR_LU,&match);CHKERRQ(ierr);
760481b7641SJakub Kruzik         if (match) {
761ac66f006SJakub Kruzik           ierr = PCFactorSetMatSolverType(pcinner,MATSOLVERSUPERLU);CHKERRQ(ierr);
762481b7641SJakub Kruzik         }
763481b7641SJakub Kruzik #endif
764481b7641SJakub Kruzik #if defined(PETSC_HAVE_SUPERLU_DIST)
765481b7641SJakub Kruzik         ierr = MatGetFactorAvailable(def->WtAW,MATSOLVERSUPERLU_DIST,MAT_FACTOR_LU,&match);CHKERRQ(ierr);
766481b7641SJakub Kruzik         if (match) {
767ac66f006SJakub Kruzik           ierr = PCFactorSetMatSolverType(pcinner,MATSOLVERSUPERLU_DIST);CHKERRQ(ierr);
768409a357bSJakub Kruzik         }
769481b7641SJakub Kruzik #endif
770409a357bSJakub Kruzik       }
771557a2f7dSJakub Kruzik     }
772557a2f7dSJakub Kruzik 
773557a2f7dSJakub Kruzik     if (innerksp) {
7746c93e71cSJakub Kruzik       if (def->prefix) {
7756c93e71cSJakub Kruzik         ierr = KSPSetOptionsPrefix(innerksp,def->prefix);CHKERRQ(ierr);
77622b0793eSJakub Kruzik         ierr = KSPAppendOptionsPrefix(innerksp,"def_");CHKERRQ(ierr);
7776c93e71cSJakub Kruzik       } else {
7786c93e71cSJakub Kruzik         ierr = KSPSetOptionsPrefix(innerksp,"def_");CHKERRQ(ierr);
7796c93e71cSJakub Kruzik       }
7806c93e71cSJakub Kruzik       ierr = KSPAppendOptionsPrefix(innerksp,prefix);CHKERRQ(ierr);
781557a2f7dSJakub Kruzik       ierr = KSPSetFromOptions(innerksp);CHKERRQ(ierr);
782557a2f7dSJakub Kruzik       ierr = KSPSetUp(innerksp);CHKERRQ(ierr);
783ac66f006SJakub Kruzik     }
784409a357bSJakub Kruzik   }
785557a2f7dSJakub Kruzik   ierr = KSPSetFromOptions(def->WtAWinv);CHKERRQ(ierr);
786557a2f7dSJakub Kruzik   ierr = KSPSetUp(def->WtAWinv);CHKERRQ(ierr);
787409a357bSJakub Kruzik 
78822b0793eSJakub Kruzik   if (!def->pc) {
78922b0793eSJakub Kruzik     ierr = PCCreate(comm,&def->pc);CHKERRQ(ierr);
79022b0793eSJakub Kruzik     ierr = PCSetOperators(def->pc,Amat,Amat);CHKERRQ(ierr);
79122b0793eSJakub Kruzik     ierr = PCSetType(def->pc,PCNONE);CHKERRQ(ierr);
7926c93e71cSJakub Kruzik     if (def->prefix) {
7936c93e71cSJakub Kruzik       ierr = PCSetOptionsPrefix(def->pc,def->prefix);CHKERRQ(ierr);
79422b0793eSJakub Kruzik     }
7956c93e71cSJakub Kruzik     ierr = PCAppendOptionsPrefix(def->pc,"def_");CHKERRQ(ierr);
7966c93e71cSJakub Kruzik     ierr = PCAppendOptionsPrefix(def->pc,prefix);CHKERRQ(ierr);
7976c93e71cSJakub Kruzik     ierr = PCAppendOptionsPrefix(def->pc,"pc_");CHKERRQ(ierr);
79822b0793eSJakub Kruzik     ierr = PCSetFromOptions(def->pc);CHKERRQ(ierr);
79922b0793eSJakub Kruzik     ierr = PCSetUp(def->pc);CHKERRQ(ierr);
80022b0793eSJakub Kruzik   }
80122b0793eSJakub Kruzik 
802f8bf229cSJakub Kruzik   /* create work vecs */
803f8bf229cSJakub Kruzik   ierr = MatCreateVecs(Amat,NULL,&def->work);CHKERRQ(ierr);
804f8bf229cSJakub Kruzik   ierr = KSPCreateVecs(def->WtAWinv,2,&def->workcoarse,0,NULL);CHKERRQ(ierr);
80537eeb815SJakub Kruzik   PetscFunctionReturn(0);
80637eeb815SJakub Kruzik }
807b30d4775SJakub Kruzik 
80837eeb815SJakub Kruzik static PetscErrorCode PCReset_Deflation(PC pc)
80937eeb815SJakub Kruzik {
810b30d4775SJakub Kruzik   PC_Deflation      *def = (PC_Deflation*)pc->data;
81137eeb815SJakub Kruzik   PetscErrorCode    ierr;
81237eeb815SJakub Kruzik 
81337eeb815SJakub Kruzik   PetscFunctionBegin;
814b30d4775SJakub Kruzik   ierr = VecDestroy(&def->work);CHKERRQ(ierr);
815b30d4775SJakub Kruzik   ierr = VecDestroyVecs(2,&def->workcoarse);CHKERRQ(ierr);
816b30d4775SJakub Kruzik   ierr = MatDestroy(&def->W);CHKERRQ(ierr);
817b30d4775SJakub Kruzik   ierr = MatDestroy(&def->Wt);CHKERRQ(ierr);
818ae029463SJakub Kruzik   ierr = MatDestroy(&def->WtA);CHKERRQ(ierr);
819b30d4775SJakub Kruzik   ierr = MatDestroy(&def->WtAW);CHKERRQ(ierr);
820b30d4775SJakub Kruzik   ierr = KSPDestroy(&def->WtAWinv);CHKERRQ(ierr);
82122b0793eSJakub Kruzik   ierr = PCDestroy(&def->pc);CHKERRQ(ierr);
82237eeb815SJakub Kruzik   PetscFunctionReturn(0);
82337eeb815SJakub Kruzik }
82437eeb815SJakub Kruzik 
82537eeb815SJakub Kruzik static PetscErrorCode PCDestroy_Deflation(PC pc)
82637eeb815SJakub Kruzik {
82737eeb815SJakub Kruzik   PetscErrorCode ierr;
82837eeb815SJakub Kruzik 
82937eeb815SJakub Kruzik   PetscFunctionBegin;
83037eeb815SJakub Kruzik   ierr = PCReset_Deflation(pc);CHKERRQ(ierr);
83137eeb815SJakub Kruzik   ierr = PetscFree(pc->data);CHKERRQ(ierr);
83237eeb815SJakub Kruzik   PetscFunctionReturn(0);
83337eeb815SJakub Kruzik }
83437eeb815SJakub Kruzik 
8358513960bSJakub Kruzik static PetscErrorCode PCView_Deflation(PC pc,PetscViewer viewer)
8368513960bSJakub Kruzik {
8378513960bSJakub Kruzik   PC_Deflation      *def = (PC_Deflation*)pc->data;
8381ac6250aSJakub Kruzik   PetscInt          its;
8398513960bSJakub Kruzik   PetscBool         iascii;
8408513960bSJakub Kruzik   PetscErrorCode    ierr;
8418513960bSJakub Kruzik 
8428513960bSJakub Kruzik   PetscFunctionBegin;
8438513960bSJakub Kruzik   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
8448513960bSJakub Kruzik   if (iascii) {
8458513960bSJakub Kruzik     if (def->correct) {
8461ac6250aSJakub Kruzik       ierr = PetscViewerASCIIPrintf(viewer,"using CP correction, factor = %g+%gi\n",
8471ac6250aSJakub Kruzik                                     (double)PetscRealPart(def->correctfact),
8481ac6250aSJakub Kruzik                                     (double)PetscImaginaryPart(def->correctfact));CHKERRQ(ierr);
8498513960bSJakub Kruzik     }
8506c93e71cSJakub Kruzik     if (!def->lvl) {
8511ac6250aSJakub Kruzik       ierr = PetscViewerASCIIPrintf(viewer,"deflation space type: %s\n",PCDeflationSpaceTypes[def->spacetype]);CHKERRQ(ierr);
8528513960bSJakub Kruzik     }
8538513960bSJakub Kruzik 
8541ac6250aSJakub Kruzik     ierr = PetscViewerASCIIPrintf(viewer,"--- Additional PC:\n");CHKERRQ(ierr);
85522b0793eSJakub Kruzik     ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
85622b0793eSJakub Kruzik     ierr = PCView(def->pc,viewer);CHKERRQ(ierr);
85722b0793eSJakub Kruzik     ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
85822b0793eSJakub Kruzik 
8591ac6250aSJakub Kruzik     ierr = PetscViewerASCIIPrintf(viewer,"--- Coarse problem solver:\n");CHKERRQ(ierr);
8608513960bSJakub Kruzik     ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
8611ac6250aSJakub Kruzik     ierr = KSPGetTotalIterations(def->WtAWinv,&its);CHKERRQ(ierr);
8621ac6250aSJakub Kruzik     ierr = PetscViewerASCIIPrintf(viewer,"total number of iterations: %D\n",its);CHKERRQ(ierr);
8638513960bSJakub Kruzik     ierr = KSPView(def->WtAWinv,viewer);CHKERRQ(ierr);
8648513960bSJakub Kruzik     ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
8658513960bSJakub Kruzik   }
8668513960bSJakub Kruzik   PetscFunctionReturn(0);
8678513960bSJakub Kruzik }
8688513960bSJakub Kruzik 
86937eeb815SJakub Kruzik static PetscErrorCode PCSetFromOptions_Deflation(PetscOptionItems *PetscOptionsObject,PC pc)
87037eeb815SJakub Kruzik {
871880d05ceSJakub Kruzik   PC_Deflation      *def = (PC_Deflation*)pc->data;
87237eeb815SJakub Kruzik   PetscErrorCode    ierr;
87337eeb815SJakub Kruzik 
87437eeb815SJakub Kruzik   PetscFunctionBegin;
87537eeb815SJakub Kruzik   ierr = PetscOptionsHead(PetscOptionsObject,"Deflation options");CHKERRQ(ierr);
876a122ebaeSJakub Kruzik   ierr = PetscOptionsBool("-pc_deflation_init_only","Use only initialization step - Initdef","PCDeflationSetInitOnly",def->init,&def->init,NULL);CHKERRQ(ierr);
8776c93e71cSJakub Kruzik   ierr = PetscOptionsInt("-pc_deflation_max_lvl","Maximum of deflation levels","PCDeflationSetMaxLvl",def->maxlvl,&def->maxlvl,NULL);CHKERRQ(ierr);
878859a873cSJakub Kruzik   ierr = PetscOptionsInt("-pc_deflation_reduction_factor","Reduction factor for coarse problem solution using PCTELESCOPE","PCDeflationSetReductionFactor",def->reductionfact,&def->reductionfact,NULL);CHKERRQ(ierr);
8798a71cb68SJakub Kruzik   ierr = PetscOptionsBool("-pc_deflation_correction","Add coarse problem correction Q to P","PCDeflationSetCorrectionFactor",def->correct,&def->correct,NULL);CHKERRQ(ierr);
8808a71cb68SJakub Kruzik   ierr = PetscOptionsScalar("-pc_deflation_correction_factor","Set multiple of Q to use as coarse problem correction","PCDeflationSetCorrectionFactor",def->correctfact,&def->correctfact,NULL);CHKERRQ(ierr);
881880d05ceSJakub Kruzik   ierr = PetscOptionsEnum("-pc_deflation_compute_space","Compute deflation space","PCDeflationSetSpace",PCDeflationSpaceTypes,(PetscEnum)def->spacetype,(PetscEnum*)&def->spacetype,NULL);CHKERRQ(ierr);
882880d05ceSJakub Kruzik   ierr = PetscOptionsInt("-pc_deflation_compute_space_size","Set size of the deflation space to compute","PCDeflationSetSpace",def->spacesize,&def->spacesize,NULL);CHKERRQ(ierr);
883880d05ceSJakub Kruzik   ierr = PetscOptionsBool("-pc_deflation_space_extend","Extend deflation space instead of truncating (wavelets)","PCDeflation",def->extendsp,&def->extendsp,NULL);CHKERRQ(ierr);
884880d05ceSJakub Kruzik //TODO add set function and fix manpages
88537eeb815SJakub Kruzik   ierr = PetscOptionsTail();CHKERRQ(ierr);
88637eeb815SJakub Kruzik   PetscFunctionReturn(0);
88737eeb815SJakub Kruzik }
88837eeb815SJakub Kruzik 
88937eeb815SJakub Kruzik /*MC
890e26ad46dSJakub Kruzik    PCDEFLATION - Deflation preconditioner shifts (deflates) part of the spectrum to zero or to a predefined value.
89137eeb815SJakub Kruzik 
892e26ad46dSJakub Kruzik    Options Database Keys:
893e26ad46dSJakub Kruzik +    -pc_deflation_init_only          <false> - if true computes only the special guess
894e26ad46dSJakub Kruzik .    -pc_deflation_max_lvl            <0>     - maximum number of levels for multilevel deflation
895e26ad46dSJakub Kruzik .    -pc_deflation_reduction_factor   <-1>    - reduction factor on bottom level coarse problem for PCTELESCOPE (default based on the size of the coarse problem)
896e26ad46dSJakub Kruzik .    -pc_deflation_correction         <false> - if true apply coarse problem correction
897e26ad46dSJakub Kruzik .    -pc_deflation_correction_factor  <1.0>   - sets coarse problem correction factor
898e26ad46dSJakub Kruzik .    -pc_deflation_compute_space      <haar>  - compute PCDeflationSpaceType deflation space
899e26ad46dSJakub Kruzik -    -pc_deflation_compute_space_size <1>     - size of the deflation space (corresponds to number of levels for wavelet-based deflation)
90037eeb815SJakub Kruzik 
90137eeb815SJakub Kruzik    Notes:
902e26ad46dSJakub Kruzik     Given a (complex - transpose is always hermitian) full rank deflation matrix W, the deflation (introduced in [1,2])
903e26ad46dSJakub Kruzik     preconditioner uses projections Q = W*(W'*A*W)^{-1}*W' and P = I - Q*A, where A is pmat.
904e26ad46dSJakub Kruzik 
905e26ad46dSJakub Kruzik     The deflation computes initial guess x0 = x_{-1} - Q*r_{-1}, which is the solution on the deflation space.
906e26ad46dSJakub Kruzik     If PCDeflationSetInitOnly() or -pc_deflation_init_only is set to PETSC_TRUE (InitDef scheme), the application of the
907e26ad46dSJakub Kruzik     preconditioner consists only of application of the additional preconditioner M^{-1}. Otherwise, the preconditioner
908e26ad46dSJakub Kruzik     application consists of P*M{-1} + factor*Q. The first part of the preconditioner (PM^{-1}) shifts some eigenvalues
909e26ad46dSJakub Kruzik     to zero while the addition of the coarse problem correction (factor*Q) makes the preconditioner to shift some
910e26ad46dSJakub Kruzik     eigenvalues to the given factor. The InitDef scheme is recommended for deflation using high accuracy estimates
911e26ad46dSJakub Kruzik     of eigenvectors of A when it exhibits similar convergence to the full deflation but is cheaper.
912e26ad46dSJakub Kruzik 
913e26ad46dSJakub Kruzik     The deflation matrix is by default automatically computed. The type of deflation matrix and its size to compute can
914e26ad46dSJakub Kruzik     be controlled by PCDeflationSetSpaceToCompute() or -pc_deflation_compute_space and -pc_deflation_compute_space_size.
915e26ad46dSJakub Kruzik     User can set an arbitrary deflation space matrix with PCDeflationSetSpace(). If the deflation matrix
916e26ad46dSJakub Kruzik     is a multiplicative MATCOMPOSITE, a multilevel deflation [3] is used. The first matrix in the composite is used as the
917e26ad46dSJakub Kruzik     deflation matrix, and the coarse problem (W'*A*W)^{-1} is solved by FCG (if A is MAT_SPD) or FGMRES preconditioned
918e26ad46dSJakub Kruzik     by deflation with deflation matrix being the next matrix in the MATCOMPOSITE. This scheme repeats until the maximum
919e26ad46dSJakub Kruzik     level is reached or there are no more matrices. If the maximal level is reached, the remaining matrices are merged
920e26ad46dSJakub Kruzik     (multiplied) to create the last deflation matrix. The maximal level defaults to 0 and can be set by
921e26ad46dSJakub Kruzik     PCDeflationSetMaxLvl() or by -pc_deflation_max_lvl.
922e26ad46dSJakub Kruzik 
9236c93e71cSJakub Kruzik     The coarse problem KSP can be controlled from the command line with prefix -def_ for the first level and -def_[lvl-1]
9246c93e71cSJakub Kruzik     from the second level onward. You can also use
925e26ad46dSJakub Kruzik     PCDeflationGetCoarseKSP() or PCDeflationSetCoarseKSP() to control it from code. The bottom level KSP defaults to
926481b7641SJakub Kruzik     KSPREONLY with PCLU direct solver (MATSOLVERSUPERLU/MATSOLVERSUPERLU_DIST if available) wrapped into PCTELESCOPE.
927481b7641SJakub Kruzik     For convenience, the reduction factor can be set by PCDeflationSetReductionFactor()
928481b7641SJakub Kruzik     or -pc_deflation_recduction_factor. The default is chosen heuristically based on the coarse problem size.
929e26ad46dSJakub Kruzik 
9306c93e71cSJakub Kruzik     The additional preconditioner can be controlled from command line with prefix -def_[lvl]_pc (same rules used for
9316c93e71cSJakub Kruzik     coarse problem KSP apply for [lvl]_ part of prefix), e.g., -def_1_pc_pc_type bjacobi. You can also use
9326c93e71cSJakub Kruzik     PCDeflationGetPC() or PCDeflationSetPC() to control the additional preconditioner from code. It defaults to PCNONE.
933e26ad46dSJakub Kruzik 
934e26ad46dSJakub Kruzik     The coarse problem correction term (factor*Q) can be turned on by -pc_deflation_correction and the factor value can
935e26ad46dSJakub Kruzik     be set by pc_deflation_correction_factor or by PCDeflationSetCorrectionFactor(). The coarse problem can
936e26ad46dSJakub Kruzik     significantly improve convergence when the deflation coarse problem is not solved with high enough accuracy. We
937e26ad46dSJakub Kruzik     recommend setting factor to some eigenvalue, e.g., the largest eigenvalue so that the preconditioner does not create
938e26ad46dSJakub Kruzik     an isolated eigenvalue.
939e26ad46dSJakub Kruzik 
9409f604af8SJakub Kruzik     The options are automatically inherited from previous deflation level.
9419f604af8SJakub Kruzik 
942e26ad46dSJakub Kruzik     The preconditioner supports KSPMonitorDynamicTolerance(). This is useful for the multilevel scheme for which we also
943e26ad46dSJakub Kruzik     recommend limiting the number of iterations for the coarse problem.
944e26ad46dSJakub Kruzik 
945e26ad46dSJakub Kruzik     See section 3 of [4] for additional references and decription of the algorithm when used for conjugate gradients.
946e26ad46dSJakub Kruzik     Section 4 describes some possible choices for the deflation space.
947e26ad46dSJakub Kruzik 
948e26ad46dSJakub Kruzik    Developer Notes:
949e26ad46dSJakub Kruzik      Contributed by Jakub Kruzik (PERMON), Institute of Geonics of the Czech
950e26ad46dSJakub Kruzik      Academy of Sciences and VSB - TU Ostrava.
951e26ad46dSJakub Kruzik 
952e26ad46dSJakub Kruzik      Developed from PERMON code used in [4] while on a research stay with
953e26ad46dSJakub Kruzik      Prof. Reinhard Nabben at the Institute of Mathematics, TU Berlin.
954e26ad46dSJakub Kruzik 
955e26ad46dSJakub Kruzik    References:
956e26ad46dSJakub Kruzik +    [1] - A. Nicolaides. “Deflation of conjugate gradients with applications to boundary valueproblems”, SIAM J. Numer. Anal. 24.2, 1987.
957e26ad46dSJakub Kruzik .    [2] - Z. Dostal. "Conjugate gradient method with preconditioning by projector", Int J. Comput. Math. 23.3-4, 1988.
958e26ad46dSJakub Kruzik .    [3] - Y. A. Erlangga and R. Nabben. "Multilevel Projection-Based Nested Krylov Iteration for Boundary Value Problems", SIAM J. Sci. Comput. 30.3, 2008.
959e26ad46dSJakub Kruzik .    [4] - J. Kruzik "Implementation of the Deflated Variants of the Conjugate Gradient Method", Master's thesis, VSB-TUO, 2018 - http://dspace5.vsb.cz/bitstream/handle/10084/130303/KRU0097_USP_N2658_2612T078_2018.pdf
960e26ad46dSJakub Kruzik 
961e26ad46dSJakub Kruzik    Level: intermediate
96237eeb815SJakub Kruzik 
96337eeb815SJakub Kruzik .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
964e26ad46dSJakub Kruzik            PCDeflationSetInitOnly(), PCDeflationSetMaxLvl(), PCDeflationSetReductionFactor(),
965e26ad46dSJakub Kruzik            PCDeflationSetCorrectionFactor(), PCDeflationSetSpaceToCompoute(),
966e26ad46dSJakub Kruzik            PCDeflationSetSpace(), PCDeflationSpaceType, PCDeflationSetProjectionNullSpaceMat(),
967e26ad46dSJakub Kruzik            PCDeflationSetCoarseMat(), PCDeflationSetCoarseKSP(), PCDeflationGetCoarseKSP(),
968e26ad46dSJakub Kruzik            PCDeflationGetPC(), PCDeflationSetPC()
96937eeb815SJakub Kruzik M*/
97037eeb815SJakub Kruzik 
97137eeb815SJakub Kruzik PETSC_EXTERN PetscErrorCode PCCreate_Deflation(PC pc)
97237eeb815SJakub Kruzik {
97337eeb815SJakub Kruzik   PC_Deflation   *def;
97437eeb815SJakub Kruzik   PetscErrorCode ierr;
97537eeb815SJakub Kruzik 
97637eeb815SJakub Kruzik   PetscFunctionBegin;
97737eeb815SJakub Kruzik   ierr     = PetscNewLog(pc,&def);CHKERRQ(ierr);
97837eeb815SJakub Kruzik   pc->data = (void*)def;
97937eeb815SJakub Kruzik 
980e662bc50SJakub Kruzik   def->init          = PETSC_FALSE;
981e662bc50SJakub Kruzik   def->correct       = PETSC_FALSE;
982fcb31d99SJakub Kruzik   def->correctfact   = 1.0;
983e662bc50SJakub Kruzik   def->reductionfact = -1;
984e662bc50SJakub Kruzik   def->spacetype     = PC_DEFLATION_SPACE_HAAR;
985e662bc50SJakub Kruzik   def->spacesize     = 1;
986e662bc50SJakub Kruzik   def->extendsp      = PETSC_FALSE;
9876c93e71cSJakub Kruzik   def->lvl           = 0;
9886c93e71cSJakub Kruzik   def->maxlvl        = 0;
98937eeb815SJakub Kruzik 
99037eeb815SJakub Kruzik   pc->ops->apply          = PCApply_Deflation;
991b27c8092SJakub Kruzik   pc->ops->presolve       = PCPreSolve_Deflation;
99237eeb815SJakub Kruzik   pc->ops->setup          = PCSetUp_Deflation;
99337eeb815SJakub Kruzik   pc->ops->reset          = PCReset_Deflation;
99437eeb815SJakub Kruzik   pc->ops->destroy        = PCDestroy_Deflation;
99537eeb815SJakub Kruzik   pc->ops->setfromoptions = PCSetFromOptions_Deflation;
9968513960bSJakub Kruzik   pc->ops->view           = PCView_Deflation;
99737eeb815SJakub Kruzik 
998a122ebaeSJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetInitOnly_C",PCDeflationSetInitOnly_Deflation);CHKERRQ(ierr);
99998d6e3deSJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetLvl_C",PCDeflationSetLvl_Deflation);CHKERRQ(ierr);
1000859a873cSJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetReductionFactor_C",PCDeflationSetReductionFactor_Deflation);CHKERRQ(ierr);
10018a71cb68SJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetCorrectionFactor_C",PCDeflationSetCorrectionFactor_Deflation);CHKERRQ(ierr);
100239417d7eSJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetSpaceToCompute_C",PCDeflationSetSpaceToCompute_Deflation);CHKERRQ(ierr);
1003e662bc50SJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetSpace_C",PCDeflationSetSpace_Deflation);CHKERRQ(ierr);
10043cf3a049SJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetProjectionNullSpaceMat_C",PCDeflationSetProjectionNullSpaceMat_Deflation);CHKERRQ(ierr);
1005e924b002SJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetCoarseMat_C",PCDeflationSetCoarseMat_Deflation);CHKERRQ(ierr);
10064a99276eSJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationGetCoarseKSP_C",PCDeflationGetCoarseKSP_Deflation);CHKERRQ(ierr);
1007f3eaa4f0SJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetCoarseKSP_C",PCDeflationSetCoarseKSP_Deflation);CHKERRQ(ierr);
100898d6e3deSJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationGetPC_C",PCDeflationGetPC_Deflation);CHKERRQ(ierr);
100998d6e3deSJakub Kruzik   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetPC_C",PCDeflationSetPC_Deflation);CHKERRQ(ierr);
101037eeb815SJakub Kruzik   PetscFunctionReturn(0);
101137eeb815SJakub Kruzik }
101237eeb815SJakub Kruzik 
1013