xref: /petsc/src/ksp/ksp/impls/fetidp/fetidp.c (revision 547c9a8eef4c37d1d4925374a9f807749d52415c) 
1 #include <petsc/private/kspimpl.h> /*I <petscksp.h> I*/
2 #include <../src/ksp/pc/impls/bddc/bddc.h>
3 
4 /*
5     This file implements the FETI-DP method in PETSc as part of KSP.
6 */
7 typedef struct {
8   KSP parentksp;
9 } KSP_FETIDPMon;
10 
11 typedef struct {
12   KSP           innerksp;         /* the KSP for the Lagrange multipliers */
13   PC            innerbddc;        /* the inner BDDC object */
14   PetscBool     fully_redundant;  /* true for using a fully redundant set of multipliers */
15   PetscBool     userbddc;         /* true if the user provided the PCBDDC object */
16   PetscBool     saddlepoint;      /* support for saddle point problems */
17   KSP_FETIDPMon *monctx;          /* monitor context, used to pass user defined monitors
18                                      in the physical space */
19 } KSP_FETIDP;
20 
21 #undef __FUNCT__
22 #define __FUNCT__ "KSPFETIDPSetPressureOperators_FETIDP"
23 static PetscErrorCode KSPFETIDPSetPressureOperators_FETIDP(KSP ksp, Mat A, Mat P)
24 {
25   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
26   PetscErrorCode ierr;
27 
28   PetscFunctionBegin;
29   if (A || P) fetidp->saddlepoint = PETSC_TRUE;
30   ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_AAmat",(PetscObject)A);CHKERRQ(ierr);
31   ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PAmat",(PetscObject)P);CHKERRQ(ierr);
32   PetscFunctionReturn(0);
33 }
34 
35 #undef __FUNCT__
36 #define __FUNCT__ "KSPFETIDPSetPressureOperators"
37 /*@
38  KSPFETIDPSetPressureOperators - Sets the operators used to setup the pressure preconditioner for saddle point FETI-DP.
39 
40    Collective on KSP
41 
42    Input Parameters:
43 +  ksp - the FETI-DP Krylov solver
44 .  A - the linear operator defining the pressure problem
45 -  P - the linear operator to be preconditioned
46 
47    Level: advanced
48 
49    Notes: The operators can be either passed in monolithic global ordering or in interface pressure ordering.
50           In the latter case, the interface pressure ordering of dofs needs to satisfy
51              pid_1 < pid_2  iff  gid_1 < gid_2
52           where pid_1 and pid_2 are two different pressure dof numbers and gid_1 and gid_2 the corresponding
53           id in the global ordering.
54 
55 .seealso: MATIS, PCBDDC, KSPFETIDPGetInnerBDDC, KSPFETIDPGetInnerKSP, KSPSetOperators
56 @*/
57 PetscErrorCode KSPFETIDPSetPressureOperators(KSP ksp, Mat A, Mat P)
58 {
59   PetscErrorCode ierr;
60 
61   PetscFunctionBegin;
62   PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
63   if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,2);
64   if (P) PetscValidHeaderSpecific(P,MAT_CLASSID,3);
65   ierr = PetscTryMethod(ksp,"KSPFETIDPSetPressureOperators_C",(KSP,Mat,Mat),(ksp,A,P));CHKERRQ(ierr);
66   PetscFunctionReturn(0);
67 }
68 
69 #undef __FUNCT__
70 #define __FUNCT__ "KSPFETIDPGetInnerKSP_FETIDP"
71 static PetscErrorCode KSPFETIDPGetInnerKSP_FETIDP(KSP ksp, KSP* innerksp)
72 {
73   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
74 
75   PetscFunctionBegin;
76   *innerksp = fetidp->innerksp;
77   PetscFunctionReturn(0);
78 }
79 
80 #undef __FUNCT__
81 #define __FUNCT__ "KSPFETIDPGetInnerKSP"
82 /*@
83  KSPFETIDPGetInnerKSP - Gets the KSP object for the Lagrange multipliers
84 
85    Input Parameters:
86 +  ksp - the FETI-DP KSP
87 -  innerksp - the KSP for the multipliers
88 
89    Level: advanced
90 
91    Notes:
92 
93 .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC, KSPFETIDPGetInnerBDDC
94 @*/
95 PetscErrorCode KSPFETIDPGetInnerKSP(KSP ksp, KSP* innerksp)
96 {
97   PetscErrorCode ierr;
98 
99   PetscFunctionBegin;
100   PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
101   PetscValidPointer(innerksp,2);
102   ierr = PetscUseMethod(ksp,"KSPFETIDPGetInnerKSP_C",(KSP,KSP*),(ksp,innerksp));CHKERRQ(ierr);
103   PetscFunctionReturn(0);
104 }
105 
106 #undef __FUNCT__
107 #define __FUNCT__ "KSPFETIDPGetInnerBDDC_FETIDP"
108 static PetscErrorCode KSPFETIDPGetInnerBDDC_FETIDP(KSP ksp, PC* pc)
109 {
110   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
111 
112   PetscFunctionBegin;
113   *pc = fetidp->innerbddc;
114   PetscFunctionReturn(0);
115 }
116 
117 #undef __FUNCT__
118 #define __FUNCT__ "KSPFETIDPGetInnerBDDC"
119 /*@
120  KSPFETIDPGetInnerBDDC - Gets the BDDC preconditioner used to setup the FETI-DP matrix for the Lagrange multipliers
121 
122    Input Parameters:
123 +  ksp - the FETI-DP Krylov solver
124 -  pc - the BDDC preconditioner
125 
126    Level: advanced
127 
128    Notes:
129 
130 .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC, KSPFETIDPGetInnerKSP
131 @*/
132 PetscErrorCode KSPFETIDPGetInnerBDDC(KSP ksp, PC* pc)
133 {
134   PetscErrorCode ierr;
135 
136   PetscFunctionBegin;
137   PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
138   PetscValidPointer(pc,2);
139   ierr = PetscUseMethod(ksp,"KSPFETIDPGetInnerBDDC_C",(KSP,PC*),(ksp,pc));CHKERRQ(ierr);
140   PetscFunctionReturn(0);
141 }
142 
143 #undef __FUNCT__
144 #define __FUNCT__ "KSPFETIDPSetInnerBDDC_FETIDP"
145 static PetscErrorCode KSPFETIDPSetInnerBDDC_FETIDP(KSP ksp, PC pc)
146 {
147   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
148   PetscErrorCode ierr;
149 
150   PetscFunctionBegin;
151   ierr = PetscObjectReference((PetscObject)pc);CHKERRQ(ierr);
152   ierr = PCDestroy(&fetidp->innerbddc);CHKERRQ(ierr);
153   fetidp->innerbddc = pc;
154   fetidp->userbddc  = PETSC_TRUE;
155   PetscFunctionReturn(0);
156 }
157 
158 #undef __FUNCT__
159 #define __FUNCT__ "KSPFETIDPSetInnerBDDC"
160 /*@
161  KSPFETIDPSetInnerBDDC - Sets the BDDC preconditioner used to setup the FETI-DP matrix for the Lagrange multipliers
162 
163    Collective on KSP
164 
165    Input Parameters:
166 +  ksp - the FETI-DP Krylov solver
167 -  pc - the BDDC preconditioner
168 
169    Level: advanced
170 
171    Notes:
172 
173 .seealso: MATIS, PCBDDC, KSPFETIDPGetInnerBDDC, KSPFETIDPGetInnerKSP
174 @*/
175 PetscErrorCode KSPFETIDPSetInnerBDDC(KSP ksp, PC pc)
176 {
177   PetscBool      isbddc;
178   PetscErrorCode ierr;
179 
180   PetscFunctionBegin;
181   PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
182   PetscValidHeaderSpecific(pc,PC_CLASSID,2);
183   ierr = PetscObjectTypeCompare((PetscObject)pc,PCBDDC,&isbddc);CHKERRQ(ierr);
184   if (!isbddc) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONG,"KSPFETIDPSetInnerBDDC need a PCBDDC preconditioner");
185   ierr = PetscTryMethod(ksp,"KSPFETIDPSetInnerBDDC_C",(KSP,PC),(ksp,pc));CHKERRQ(ierr);
186   PetscFunctionReturn(0);
187 }
188 
189 #undef __FUNCT__
190 #define __FUNCT__ "KSPBuildSolution_FETIDP"
191 static PetscErrorCode KSPBuildSolution_FETIDP(KSP ksp,Vec v,Vec *V)
192 {
193   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
194   Mat            F;
195   Vec            Xl;
196   PetscErrorCode ierr;
197 
198   PetscFunctionBegin;
199   ierr = KSPGetOperators(fetidp->innerksp,&F,NULL);CHKERRQ(ierr);
200   ierr = KSPBuildSolution(fetidp->innerksp,NULL,&Xl);CHKERRQ(ierr);
201   if (v) {
202     ierr = PCBDDCMatFETIDPGetSolution(F,Xl,v);CHKERRQ(ierr);
203     *V   = v;
204   } else {
205     ierr = PCBDDCMatFETIDPGetSolution(F,Xl,*V);CHKERRQ(ierr);
206   }
207   PetscFunctionReturn(0);
208 }
209 
210 #undef __FUNCT__
211 #define __FUNCT__ "KSPMonitor_FETIDP"
212 static PetscErrorCode KSPMonitor_FETIDP(KSP ksp,PetscInt it,PetscReal rnorm,void* ctx)
213 {
214   KSP_FETIDPMon  *monctx = (KSP_FETIDPMon*)ctx;
215   PetscErrorCode ierr;
216 
217   PetscFunctionBegin;
218   ierr = KSPMonitor(monctx->parentksp,it,rnorm);CHKERRQ(ierr);
219   PetscFunctionReturn(0);
220 }
221 
222 #undef __FUNCT__
223 #define __FUNCT__ "KSPComputeEigenvalues_FETIDP"
224 static PetscErrorCode KSPComputeEigenvalues_FETIDP(KSP ksp,PetscInt nmax,PetscReal *r,PetscReal *c,PetscInt *neig)
225 {
226   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
227   PetscErrorCode ierr;
228 
229   PetscFunctionBegin;
230   ierr = KSPComputeEigenvalues(fetidp->innerksp,nmax,r,c,neig);CHKERRQ(ierr);
231   PetscFunctionReturn(0);
232 }
233 
234 #undef __FUNCT__
235 #define __FUNCT__ "KSPComputeExtremeSingularValues_FETIDP"
236 static PetscErrorCode KSPComputeExtremeSingularValues_FETIDP(KSP ksp,PetscReal *emax,PetscReal *emin)
237 {
238   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
239   PetscErrorCode ierr;
240 
241   PetscFunctionBegin;
242   ierr = KSPComputeExtremeSingularValues(fetidp->innerksp,emax,emin);CHKERRQ(ierr);
243   PetscFunctionReturn(0);
244 }
245 
246 #undef __FUNCT__
247 #define __FUNCT__ "KSPFETIDPSetUpOperators"
248 static PetscErrorCode KSPFETIDPSetUpOperators(KSP ksp)
249 {
250   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
251   PC_BDDC        *pcbddc = (PC_BDDC*)fetidp->innerbddc->data;
252   Mat            A,Ap;
253   PetscInt       fid = -1;
254   PetscBool      ismatis,pisz;
255   PetscBool      flip; /* Usually, Stokes is written
256                          | A B'| | v | = | f |
257                          | B 0 | | p | = | g |
258                           If saddlepoint_flip is true, the code assumes it is written as
259                          | A B'| | v | = | f |
260                          |-B 0 | | p | = |-g |
261                        */
262   PetscErrorCode ierr;
263 
264   PetscFunctionBegin;
265   pisz = PETSC_FALSE;
266   flip = PETSC_FALSE;
267   ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->prefix,"FETI-DP options","PC");CHKERRQ(ierr);
268   ierr = PetscOptionsInt("-ksp_fetidp_pressure_field","Field id for pressures for saddle-point problems",NULL,fid,&fid,NULL);CHKERRQ(ierr);
269   ierr = PetscOptionsBool("-ksp_fetidp_pressure_iszero","Zero pressure block",NULL,pisz,&pisz,NULL);CHKERRQ(ierr);
270   ierr = PetscOptionsBool("-ksp_fetidp_saddlepoint_flip","Flip the sign of the pressure-velocity (lower-left) block",NULL,flip,&flip,NULL);CHKERRQ(ierr);
271   ierr = PetscOptionsEnd();CHKERRQ(ierr);
272 
273   fetidp->saddlepoint = (fid >= 0 ? PETSC_TRUE : fetidp->saddlepoint);
274 
275   ierr = KSPGetOperators(ksp,&A,&Ap);CHKERRQ(ierr);
276   ierr = PetscObjectTypeCompare((PetscObject)Ap,MATIS,&ismatis);CHKERRQ(ierr);
277   if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Pmat should be of type MATIS");
278 
279   /* see if MATIS has same fields attached */
280   if (!pcbddc->n_ISForDofsLocal && !pcbddc->n_ISForDofs) {
281     PetscContainer c;
282 
283     ierr = PetscObjectQuery((PetscObject)Ap,"_convert_nest_lfields",(PetscObject*)&c);CHKERRQ(ierr);
284     if (c) {
285       MatISLocalFields lf;
286       ierr = PetscContainerGetPointer(c,(void**)&lf);CHKERRQ(ierr);
287       ierr = PCBDDCSetDofsSplittingLocal(fetidp->innerbddc,lf->nr,lf->rf);CHKERRQ(ierr);
288     }
289   }
290 
291   if (!fetidp->saddlepoint) {
292     ierr = PCSetOperators(fetidp->innerbddc,A,Ap);CHKERRQ(ierr);
293   } else {
294     KSP kP;
295     Mat nA,lA;
296     Mat PAmat,PPmat;
297     Vec rhs_flip;
298     IS  pP,lP;
299 
300     ierr = MatISGetLocalMat(Ap,&lA);CHKERRQ(ierr);
301     ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",(PetscObject)lA);CHKERRQ(ierr);
302 
303     /* saved index sets */
304     pP  = NULL;
305     lP  = NULL;
306     ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP" ,(PetscObject*)&pP);CHKERRQ(ierr);
307     ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lP" ,(PetscObject*)&lP);CHKERRQ(ierr);
308 
309     /* vector for flipping rhs */
310     rhs_flip = NULL;
311     ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_flip" ,(PetscObject*)&rhs_flip);CHKERRQ(ierr);
312     if (!pP && !lP) { /* first time, need to compute boundary pressure dofs */
313       PC_IS                  *pcis = (PC_IS*)fetidp->innerbddc->data;
314       Mat_IS                 *matis = (Mat_IS*)(Ap->data);
315       ISLocalToGlobalMapping l2g;
316       IS                     II,pII,lPall,Pall,is1,is2;
317       const PetscInt         *idxs;
318       PetscInt               nl,ni,*widxs;
319       PetscInt               i,j,n_neigh,*neigh,*n_shared,**shared,*count;
320       PetscInt               rst,ren,n;
321       PetscBool              ploc;
322 
323       ierr = MatGetLocalSize(Ap,&nl,NULL);CHKERRQ(ierr);
324       ierr = MatGetOwnershipRange(Ap,&rst,&ren);CHKERRQ(ierr);
325       ierr = MatGetLocalSize(lA,&n,NULL);CHKERRQ(ierr);
326 
327       if (!pcis->is_I_local) { /* need to compute interior dofs */
328         ierr = PetscCalloc1(n,&count);CHKERRQ(ierr);
329         ierr = MatGetLocalToGlobalMapping(Ap,&l2g,NULL);CHKERRQ(ierr);
330         ierr = ISLocalToGlobalMappingGetInfo(l2g,&n_neigh,&neigh,&n_shared,&shared);CHKERRQ(ierr);
331         for (i=1;i<n_neigh;i++)
332           for (j=0;j<n_shared[i];j++)
333             count[shared[i][j]] += 1;
334         for (i=0,j=0;i<n;i++) if (!count[i]) count[j++] = i;
335         ierr = ISLocalToGlobalMappingRestoreInfo(l2g,&n_neigh,&neigh,&n_shared,&shared);CHKERRQ(ierr);
336         ierr = ISCreateGeneral(PETSC_COMM_SELF,j,count,PETSC_OWN_POINTER,&II);CHKERRQ(ierr);
337       } else {
338         ierr = PetscObjectReference((PetscObject)pcis->is_I_local);CHKERRQ(ierr);
339         II   = pcis->is_I_local;
340       }
341 
342       /* interior dofs in layout */
343       ierr = MatISSetUpSF(Ap);CHKERRQ(ierr);
344       ierr = PetscMemzero(matis->sf_leafdata,n*sizeof(PetscInt));CHKERRQ(ierr);
345       ierr = PetscMemzero(matis->sf_rootdata,nl*sizeof(PetscInt));CHKERRQ(ierr);
346       ierr = ISGetLocalSize(II,&ni);CHKERRQ(ierr);
347       ierr = ISGetIndices(II,&idxs);CHKERRQ(ierr);
348       for (i=0;i<ni;i++) matis->sf_leafdata[idxs[i]] = 1;
349       ierr = ISRestoreIndices(II,&idxs);CHKERRQ(ierr);
350       ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr);
351       ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr);
352       ierr = PetscMalloc1(PetscMax(nl,n),&widxs);CHKERRQ(ierr);
353       for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst;
354       ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&pII);CHKERRQ(ierr);
355 
356       /* pressure space at the interface */
357       Pall  = NULL;
358       lPall = NULL;
359       ploc  = PETSC_FALSE;
360       if (fid >= 0) {
361         if (pcbddc->n_ISForDofsLocal) {
362           PetscInt np;
363 
364           if (fid >= pcbddc->n_ISForDofsLocal) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Invalid field id for pressure %D, max %D",fid,pcbddc->n_ISForDofsLocal);
365           /* need a sequential IS */
366           ierr = ISGetLocalSize(pcbddc->ISForDofsLocal[fid],&np);CHKERRQ(ierr);
367           ierr = ISGetIndices(pcbddc->ISForDofsLocal[fid],&idxs);CHKERRQ(ierr);
368           ierr = ISCreateGeneral(PETSC_COMM_SELF,np,idxs,PETSC_COPY_VALUES,&lPall);CHKERRQ(ierr);
369           ierr = ISRestoreIndices(pcbddc->ISForDofsLocal[fid],&idxs);CHKERRQ(ierr);
370           ploc = PETSC_TRUE;
371         } else if (pcbddc->n_ISForDofs) {
372           if (fid >= pcbddc->n_ISForDofs) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Invalid field id for pressure %D, max %D",fid,pcbddc->n_ISForDofs);
373           ierr = PetscObjectReference((PetscObject)pcbddc->ISForDofs[fid]);CHKERRQ(ierr);
374           Pall = pcbddc->ISForDofs[fid];
375         } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Missing fields! Use PCBDDCSetDofsSplitting/Local");
376       } else { /* fallback to zero pressure block */
377         ierr = MatFindZeroDiagonals(Ap,&Pall);CHKERRQ(ierr);
378       }
379       if (ploc) {
380         ierr = ISDifference(lPall,II,&lP);CHKERRQ(ierr);
381         ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lP",(PetscObject)lP);CHKERRQ(ierr);
382       } else {
383         ierr = ISDifference(Pall,pII,&pP);CHKERRQ(ierr);
384         ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",(PetscObject)pP);CHKERRQ(ierr);
385         /* need all local pressure dofs */
386         ierr = PetscMemzero(matis->sf_leafdata,n*sizeof(PetscInt));CHKERRQ(ierr);
387         ierr = PetscMemzero(matis->sf_rootdata,nl*sizeof(PetscInt));CHKERRQ(ierr);
388         ierr = ISGetLocalSize(Pall,&ni);CHKERRQ(ierr);
389         ierr = ISGetIndices(Pall,&idxs);CHKERRQ(ierr);
390         for (i=0;i<ni;i++) matis->sf_rootdata[idxs[i]-rst] = 1;
391         ierr = ISRestoreIndices(Pall,&idxs);CHKERRQ(ierr);
392         ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr);
393         ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr);
394         for (i=0,ni=0;i<n;i++) if (matis->sf_leafdata[i]) widxs[ni++] = i;
395         ierr = ISCreateGeneral(PETSC_COMM_SELF,ni,widxs,PETSC_COPY_VALUES,&lPall);CHKERRQ(ierr);
396       }
397 
398       if (flip) {
399         PetscInt npl;
400         if (!Pall) {
401           ierr = PetscMemzero(matis->sf_leafdata,n*sizeof(PetscInt));CHKERRQ(ierr);
402           ierr = PetscMemzero(matis->sf_rootdata,nl*sizeof(PetscInt));CHKERRQ(ierr);
403           ierr = ISGetLocalSize(lPall,&ni);CHKERRQ(ierr);
404           ierr = ISGetIndices(lPall,&idxs);CHKERRQ(ierr);
405           for (i=0;i<ni;i++) matis->sf_leafdata[idxs[i]] = 1;
406           ierr = ISRestoreIndices(lPall,&idxs);CHKERRQ(ierr);
407           ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr);
408           ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr);
409           for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst;
410           ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&Pall);CHKERRQ(ierr);
411         }
412         ierr = ISGetLocalSize(Pall,&npl);CHKERRQ(ierr);
413         ierr = ISGetIndices(Pall,&idxs);CHKERRQ(ierr);
414         ierr = MatCreateVecs(Ap,NULL,&rhs_flip);CHKERRQ(ierr);
415         ierr = VecSet(rhs_flip,1.);CHKERRQ(ierr);
416         ierr = VecSetOption(rhs_flip,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
417         for (i=0;i<npl;i++) {
418           ierr = VecSetValue(rhs_flip,idxs[i],-1.,INSERT_VALUES);CHKERRQ(ierr);
419         }
420         ierr = VecAssemblyBegin(rhs_flip);CHKERRQ(ierr);
421         ierr = VecAssemblyEnd(rhs_flip);CHKERRQ(ierr);
422         ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_flip",(PetscObject)rhs_flip);CHKERRQ(ierr);
423         ierr = ISRestoreIndices(Pall,&idxs);CHKERRQ(ierr);
424       }
425       ierr = ISDestroy(&Pall);CHKERRQ(ierr);
426       ierr = ISDestroy(&pII);CHKERRQ(ierr);
427 
428       /* local interface pressures in subdomain-wise and global ordering */
429       ierr = PetscMemzero(matis->sf_leafdata,n*sizeof(PetscInt));CHKERRQ(ierr);
430       ierr = PetscMemzero(matis->sf_rootdata,nl*sizeof(PetscInt));CHKERRQ(ierr);
431       if (pP) {
432         ierr = ISGetLocalSize(pP,&ni);CHKERRQ(ierr);
433         ierr = ISGetIndices(pP,&idxs);CHKERRQ(ierr);
434         for (i=0;i<ni;i++) matis->sf_rootdata[idxs[i]-rst] = 1;
435         ierr = ISRestoreIndices(pP,&idxs);CHKERRQ(ierr);
436         ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr);
437         ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr);
438         for (i=0,ni=0;i<n;i++) if (matis->sf_leafdata[i]) widxs[ni++] = i;
439         ierr = ISLocalToGlobalMappingApply(l2g,ni,widxs,widxs+ni);CHKERRQ(ierr);
440         ierr = ISCreateGeneral(PETSC_COMM_SELF,ni,widxs,PETSC_COPY_VALUES,&lP);CHKERRQ(ierr);
441         ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lP",(PetscObject)lP);CHKERRQ(ierr);
442         ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs+ni,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr);
443         ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_gP",(PetscObject)is1);CHKERRQ(ierr);
444         ierr = ISDestroy(&is1);CHKERRQ(ierr);
445       } else {
446         ierr = ISGetLocalSize(lP,&ni);CHKERRQ(ierr);
447         ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr);
448         for (i=0;i<ni;i++)
449           if (idxs[i] >=0 && idxs[i] < n)
450             matis->sf_leafdata[idxs[i]] = 1;
451         ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr);
452         ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr);
453         ierr = ISLocalToGlobalMappingApply(l2g,ni,idxs,widxs);CHKERRQ(ierr);
454         ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr);
455         ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_gP",(PetscObject)is1);CHKERRQ(ierr);
456         ierr = ISDestroy(&is1);CHKERRQ(ierr);
457         ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr);
458         for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst;
459         ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&pP);CHKERRQ(ierr);
460         ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",(PetscObject)pP);CHKERRQ(ierr);
461       }
462       ierr = PetscFree(widxs);CHKERRQ(ierr);
463 
464       /* We may want to use a discrete harmonic solver instead
465          of a Stokes harmonic for the Dirichler preconditioner
466          Need to extract the interior pressure dofs in interior dofs ordering */
467       ierr = ISDifference(lPall,lP,&is1);CHKERRQ(ierr);
468       ierr = ISDifference(II,is1,&is2);CHKERRQ(ierr);
469       ierr = ISDestroy(&is1);CHKERRQ(ierr);
470       ierr = ISLocalToGlobalMappingCreateIS(II,&l2g);CHKERRQ(ierr);
471       ierr = ISGlobalToLocalMappingApplyIS(l2g,IS_GTOLM_DROP,is2,&is1);CHKERRQ(ierr);
472       ierr = ISGetLocalSize(is1,&i);CHKERRQ(ierr);
473       ierr = ISGetLocalSize(is2,&j);CHKERRQ(ierr);
474       if (i != j) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent local sizes %D and %D for iP",i,j);
475       ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_iP",(PetscObject)is1);CHKERRQ(ierr);
476       ierr = ISLocalToGlobalMappingDestroy(&l2g);CHKERRQ(ierr);
477       ierr = ISDestroy(&is1);CHKERRQ(ierr);
478       ierr = ISDestroy(&is2);CHKERRQ(ierr);
479       ierr = ISDestroy(&lPall);CHKERRQ(ierr);
480       ierr = ISDestroy(&II);CHKERRQ(ierr);
481 
482       /* exclude interface pressures from the inner BDDC */
483       if (pcbddc->DirichletBoundariesLocal) {
484         IS       list[2],plP,isout;
485         PetscInt np;
486 
487         /* need a parallel IS */
488         ierr = ISGetLocalSize(lP,&np);CHKERRQ(ierr);
489         ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr);
490         ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),np,idxs,PETSC_USE_POINTER,&plP);CHKERRQ(ierr);
491         list[0] = plP;
492         list[1] = pcbddc->DirichletBoundariesLocal;
493         ierr = ISConcatenate(PetscObjectComm((PetscObject)ksp),2,list,&isout);CHKERRQ(ierr);
494         ierr = ISDestroy(&plP);CHKERRQ(ierr);
495         ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr);
496         ierr = PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc,isout);CHKERRQ(ierr);
497         ierr = ISDestroy(&isout);CHKERRQ(ierr);
498       } else if (pcbddc->DirichletBoundaries) {
499         IS list[2],isout;
500 
501         list[0] = pP;
502         list[1] = pcbddc->DirichletBoundaries;
503         ierr = ISConcatenate(PetscObjectComm((PetscObject)ksp),2,list,&isout);CHKERRQ(ierr);
504         ierr = PCBDDCSetDirichletBoundaries(fetidp->innerbddc,isout);CHKERRQ(ierr);
505         ierr = ISDestroy(&isout);CHKERRQ(ierr);
506       } else {
507         IS       plP;
508         PetscInt np;
509 
510         /* need a parallel IS */
511         ierr = ISGetLocalSize(lP,&np);CHKERRQ(ierr);
512         ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr);
513         ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),np,idxs,PETSC_COPY_VALUES,&plP);CHKERRQ(ierr);
514         ierr = PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc,plP);CHKERRQ(ierr);
515         ierr = ISDestroy(&plP);CHKERRQ(ierr);
516         ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr);
517       }
518     } else {
519       if (!pP) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Missing global interface pressure field");
520       if (!lP) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Missing local interface pressure field");
521       ierr = PetscObjectReference((PetscObject)pP);CHKERRQ(ierr);
522       ierr = PetscObjectReference((PetscObject)lP);CHKERRQ(ierr);
523       if (rhs_flip) {
524         ierr = PetscObjectReference((PetscObject)rhs_flip);CHKERRQ(ierr);
525       }
526     }
527 
528     /* Set operator for inner BDDC */
529     ierr = MatDuplicate(Ap,MAT_COPY_VALUES,&nA);CHKERRQ(ierr);
530     if (rhs_flip) {
531       Mat lA2;
532 
533       ierr = MatDiagonalScale(nA,rhs_flip,NULL);CHKERRQ(ierr);
534       ierr = MatISGetLocalMat(nA,&lA);CHKERRQ(ierr);
535       ierr = MatDuplicate(lA,MAT_COPY_VALUES,&lA2);CHKERRQ(ierr);
536       ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",(PetscObject)lA2);CHKERRQ(ierr);
537       ierr = MatDestroy(&lA2);CHKERRQ(ierr);
538     }
539     ierr = MatSetOption(nA,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
540     ierr = MatZeroRowsColumnsIS(nA,pP,1.,NULL,NULL);CHKERRQ(ierr);
541     ierr = PCSetOperators(fetidp->innerbddc,nA,nA);CHKERRQ(ierr);
542     ierr = MatDestroy(&nA);CHKERRQ(ierr);
543     ierr = VecDestroy(&rhs_flip);CHKERRQ(ierr);
544 
545     /* non-zero rhs on interior dofs when applying the preconditioner */
546     pcbddc->switch_static = PETSC_TRUE;
547 
548     /* Operators for pressure preconditioner */
549     ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PAmat",(PetscObject*)&PAmat);CHKERRQ(ierr);
550     ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject*)&PPmat);CHKERRQ(ierr);
551     if (PAmat) {
552       ierr = PetscObjectReference((PetscObject)PAmat);CHKERRQ(ierr);
553     }
554     if (PPmat) {
555       ierr = PetscObjectReference((PetscObject)PPmat);CHKERRQ(ierr);
556     }
557 
558     /* Extract pressure block */
559     if (!pisz) {
560       Mat C;
561 
562       ierr = MatGetSubMatrix(Ap,pP,pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr);
563       ierr = MatScale(C,-1.);CHKERRQ(ierr);
564       ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_C",(PetscObject)C);CHKERRQ(ierr);
565       /* default operators for the interface pressure solver */
566       if (!PAmat) {
567         ierr  = PetscObjectReference((PetscObject)C);CHKERRQ(ierr);
568         PAmat = C;
569         ierr  = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PAmat",(PetscObject)PAmat);CHKERRQ(ierr);
570       }
571       if (!PPmat) {
572         ierr  = PetscObjectReference((PetscObject)C);CHKERRQ(ierr);
573         PPmat = C;
574         ierr  = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject)PPmat);CHKERRQ(ierr);
575       }
576       ierr = MatDestroy(&C);CHKERRQ(ierr);
577     }
578 
579     if (!PAmat) { /* Just use the identity */
580       PetscInt nl;
581 
582       ierr = ISGetLocalSize(pP,&nl);CHKERRQ(ierr);
583       ierr = MatCreate(PetscObjectComm((PetscObject)ksp),&PAmat);CHKERRQ(ierr);
584       ierr = MatSetSizes(PAmat,nl,nl,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
585       ierr = MatSetType(PAmat,MATAIJ);CHKERRQ(ierr);
586       ierr = MatMPIAIJSetPreallocation(PAmat,1,NULL,0,NULL);CHKERRQ(ierr);
587       ierr = MatSeqAIJSetPreallocation(PAmat,1,NULL);CHKERRQ(ierr);
588       ierr = MatAssemblyBegin(PAmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
589       ierr = MatAssemblyEnd(PAmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
590       ierr = MatShift(PAmat,1.);CHKERRQ(ierr);
591       ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PAmat",(PetscObject)PAmat);CHKERRQ(ierr);
592     } else { /* check size of pressure operator and restrict if needed */
593       PetscInt AM,PAM,PAN,pam,pan,am,an,pl;
594 
595       ierr = MatGetSize(Ap,&AM,NULL);CHKERRQ(ierr);
596       ierr = MatGetSize(PAmat,&PAM,&PAN);CHKERRQ(ierr);
597       ierr = MatGetLocalSize(PAmat,&pam,&pan);CHKERRQ(ierr);
598       ierr = MatGetLocalSize(Ap,&am,&an);CHKERRQ(ierr);
599       ierr = ISGetLocalSize(pP,&pl);CHKERRQ(ierr);
600       if (PAM != PAN) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Pressure matrix must be square, unsupported %D x %D",PAM,PAN);
601       if (pam != pan) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Local sizes of pressure matrix must be equal, unsupported %D x %D",pam,pan);
602       if (pam != am && pam != pl) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_USER,"Invalid number of local rows %D for pressure matrix! Supported are %D or %D",pam,am,pl);
603       if (pan != an && pan != pl) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_USER,"Invalid number of local columns %D for pressure matrix! Supported are %D or %D",pan,an,pl);
604       if (PAM == AM) { /* monolithic ordering, restrict to interface pressure */
605         Mat C;
606 
607         ierr  = MatGetSubMatrix(PAmat,pP,pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr);
608         ierr  = MatDestroy(&PAmat);CHKERRQ(ierr);
609         PAmat = C;
610         ierr  = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PAmat",(PetscObject)PAmat);CHKERRQ(ierr);
611       }
612     }
613     if (PPmat) {
614       PetscInt AM,PAM,PAN,pam,pan,am,an,pl;
615 
616       ierr = MatGetSize(Ap,&AM,NULL);CHKERRQ(ierr);
617       ierr = MatGetSize(PPmat,&PAM,&PAN);CHKERRQ(ierr);
618       ierr = MatGetLocalSize(PPmat,&pam,&pan);CHKERRQ(ierr);
619       ierr = MatGetLocalSize(Ap,&am,&an);CHKERRQ(ierr);
620       ierr = ISGetLocalSize(pP,&pl);CHKERRQ(ierr);
621       if (PAM != PAN) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Pressure matrix must be square, unsupported %D x %D",PAM,PAN);
622       if (pam != pan) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Local sizes of pressure matrix must be equal, unsupported %D x %D",pam,pan);
623       if (pam != am && pam != pl) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_USER,"Invalid number of local rows %D for pressure matrix! Supported are %D or %D",pam,am,pl);
624       if (pan != an && pan != pl) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_USER,"Invalid number of local columns %D for pressure matrix! Supported are %D or %D",pan,an,pl);
625       if (PAM == AM) { /* monolithic ordering, restrict to interface pressure */
626         Mat C;
627 
628         ierr  = MatGetSubMatrix(PPmat,pP,pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr);
629         ierr  = MatDestroy(&PPmat);CHKERRQ(ierr);
630         PPmat = C;
631         ierr  = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject)PPmat);CHKERRQ(ierr);
632       }
633     } else {
634       ierr  = PetscObjectReference((PetscObject)PAmat);CHKERRQ(ierr);
635       PPmat = PAmat;
636     }
637 
638     /* create pressure solver */
639     ierr = KSPCreate(PetscObjectComm((PetscObject)ksp),&kP);CHKERRQ(ierr);
640     ierr = KSPSetOperators(kP,PAmat,PPmat);CHKERRQ(ierr);
641     ierr = KSPSetOptionsPrefix(kP,((PetscObject)ksp)->prefix);CHKERRQ(ierr);
642     ierr = KSPAppendOptionsPrefix(kP,"fetidp_p_");CHKERRQ(ierr);
643     ierr = KSPSetFromOptions(kP);CHKERRQ(ierr);
644     ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PKSP",(PetscObject)kP);CHKERRQ(ierr);
645     ierr = MatDestroy(&PAmat);CHKERRQ(ierr);
646     ierr = MatDestroy(&PPmat);CHKERRQ(ierr);
647     ierr = KSPDestroy(&kP);CHKERRQ(ierr);
648 
649     ierr = ISDestroy(&lP);CHKERRQ(ierr);
650     ierr = ISDestroy(&pP);CHKERRQ(ierr);
651   }
652   PetscFunctionReturn(0);
653 }
654 
655 #undef __FUNCT__
656 #define __FUNCT__ "KSPSetUp_FETIDP"
657 static PetscErrorCode KSPSetUp_FETIDP(KSP ksp)
658 {
659   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
660   PC_BDDC        *pcbddc = (PC_BDDC*)fetidp->innerbddc->data;
661   PetscBool      flg;
662   PetscErrorCode ierr;
663 
664   PetscFunctionBegin;
665   ierr = KSPFETIDPSetUpOperators(ksp);CHKERRQ(ierr);
666   /* set up BDDC */
667   ierr = PCSetUp(fetidp->innerbddc);CHKERRQ(ierr);
668 
669   /* FETI-DP as it is implemented needs an exact coarse solver */
670   if (pcbddc->coarse_ksp) {
671     ierr = KSPSetTolerances(pcbddc->coarse_ksp,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,1000);CHKERRQ(ierr);
672     ierr = KSPSetNormType(pcbddc->coarse_ksp,KSP_NORM_DEFAULT);CHKERRQ(ierr);
673   }
674   /* FETI-DP as it is implemented needs exact local Neumann solvers */
675   ierr = KSPSetTolerances(pcbddc->ksp_R,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,1000);CHKERRQ(ierr);
676   ierr = KSPSetNormType(pcbddc->ksp_R,KSP_NORM_DEFAULT);CHKERRQ(ierr);
677 
678   /* if the primal space is changed, setup F */
679   if (pcbddc->new_primal_space || !pcbddc->coarse_size) {
680     Mat F; /* the FETI-DP matrix */
681     PC  D; /* the FETI-DP preconditioner */
682     ierr = KSPReset(fetidp->innerksp);CHKERRQ(ierr);
683     ierr = PCBDDCCreateFETIDPOperators(fetidp->innerbddc,fetidp->fully_redundant,((PetscObject)ksp)->prefix,&F,&D);CHKERRQ(ierr);
684     ierr = KSPSetOperators(fetidp->innerksp,F,F);CHKERRQ(ierr);
685     ierr = KSPSetTolerances(fetidp->innerksp,ksp->rtol,ksp->abstol,ksp->divtol,ksp->max_it);CHKERRQ(ierr);
686     ierr = KSPSetPC(fetidp->innerksp,D);CHKERRQ(ierr);
687     ierr = MatCreateVecs(F,&(fetidp->innerksp)->vec_rhs,&(fetidp->innerksp)->vec_sol);CHKERRQ(ierr);
688     ierr = MatDestroy(&F);CHKERRQ(ierr);
689     ierr = PCDestroy(&D);CHKERRQ(ierr);
690   }
691 
692   /* propagate settings to the inner solve */
693   ierr = KSPGetComputeSingularValues(ksp,&flg);CHKERRQ(ierr);
694   ierr = KSPSetComputeSingularValues(fetidp->innerksp,flg);CHKERRQ(ierr);
695   if (ksp->res_hist) {
696     ierr = KSPSetResidualHistory(fetidp->innerksp,ksp->res_hist,ksp->res_hist_max,ksp->res_hist_reset);CHKERRQ(ierr);
697   }
698   ierr = KSPSetUp(fetidp->innerksp);CHKERRQ(ierr);
699   PetscFunctionReturn(0);
700 }
701 
702 #undef __FUNCT__
703 #define __FUNCT__ "KSPSolve_FETIDP"
704 static PetscErrorCode KSPSolve_FETIDP(KSP ksp)
705 {
706   PetscErrorCode ierr;
707   Mat            F;
708   Vec            X,B,Xl,Bl;
709   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
710   PC_BDDC        *pcbddc = (PC_BDDC*)fetidp->innerbddc->data;
711 
712   PetscFunctionBegin;
713   ierr = KSPGetRhs(ksp,&B);CHKERRQ(ierr);
714   ierr = KSPGetSolution(ksp,&X);CHKERRQ(ierr);
715   ierr = KSPGetOperators(fetidp->innerksp,&F,NULL);CHKERRQ(ierr);
716   ierr = KSPGetRhs(fetidp->innerksp,&Bl);CHKERRQ(ierr);
717   ierr = KSPGetSolution(fetidp->innerksp,&Xl);CHKERRQ(ierr);
718   ierr = PCBDDCMatFETIDPGetRHS(F,B,Bl);CHKERRQ(ierr);
719   if (ksp->transpose_solve) {
720     ierr = KSPSolveTranspose(fetidp->innerksp,Bl,Xl);CHKERRQ(ierr);
721   } else {
722     ierr = KSPSolve(fetidp->innerksp,Bl,Xl);CHKERRQ(ierr);
723   }
724   ierr = PCBDDCMatFETIDPGetSolution(F,Xl,X);CHKERRQ(ierr);
725   /* update ksp with stats from inner ksp */
726   ierr = KSPGetConvergedReason(fetidp->innerksp,&ksp->reason);CHKERRQ(ierr);
727   ierr = KSPGetIterationNumber(fetidp->innerksp,&ksp->its);CHKERRQ(ierr);
728   ksp->totalits += ksp->its;
729   ierr = KSPGetResidualHistory(fetidp->innerksp,NULL,&ksp->res_hist_len);CHKERRQ(ierr);
730   /* restore defaults for inner BDDC (Pre/PostSolve flags) */
731   pcbddc->temp_solution_used        = PETSC_FALSE;
732   pcbddc->rhs_change                = PETSC_FALSE;
733   pcbddc->exact_dirichlet_trick_app = PETSC_FALSE;
734   PetscFunctionReturn(0);
735 }
736 
737 #undef __FUNCT__
738 #define __FUNCT__ "KSPDestroy_FETIDP"
739 static PetscErrorCode KSPDestroy_FETIDP(KSP ksp)
740 {
741   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
742   PetscErrorCode ierr;
743 
744   PetscFunctionBegin;
745   ierr = PCDestroy(&fetidp->innerbddc);CHKERRQ(ierr);
746   ierr = KSPDestroy(&fetidp->innerksp);CHKERRQ(ierr);
747   ierr = PetscFree(fetidp->monctx);CHKERRQ(ierr);
748   ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetInnerBDDC_C",NULL);CHKERRQ(ierr);
749   ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerBDDC_C",NULL);CHKERRQ(ierr);
750   ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerKSP_C",NULL);CHKERRQ(ierr);
751   ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetPressureOperators_C",NULL);CHKERRQ(ierr);
752   ierr = PetscFree(ksp->data);CHKERRQ(ierr);
753   PetscFunctionReturn(0);
754 }
755 
756 #undef __FUNCT__
757 #define __FUNCT__ "KSPView_FETIDP"
758 static PetscErrorCode KSPView_FETIDP(KSP ksp,PetscViewer viewer)
759 {
760   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
761   PetscErrorCode ierr;
762   PetscBool      iascii;
763 
764   PetscFunctionBegin;
765   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
766   if (iascii) {
767     ierr = PetscViewerASCIIPrintf(viewer,"  FETI-DP: fully redundant: %D\n",fetidp->fully_redundant);CHKERRQ(ierr);
768     ierr = PetscViewerASCIIPrintf(viewer,"  FETI-DP: saddle point:    %D\n",fetidp->saddlepoint);CHKERRQ(ierr);
769     ierr = PetscViewerASCIIPrintf(viewer,"  FETI-DP: inner solver details\n");CHKERRQ(ierr);
770     ierr = PetscViewerASCIIAddTab(viewer,2);CHKERRQ(ierr);
771   }
772   ierr = KSPView(fetidp->innerksp,viewer);CHKERRQ(ierr);
773   if (iascii) {
774     ierr = PetscViewerASCIISubtractTab(viewer,2);CHKERRQ(ierr);
775     ierr = PetscViewerASCIIPrintf(viewer,"  FETI-DP: BDDC solver details\n");CHKERRQ(ierr);
776     ierr = PetscViewerASCIIAddTab(viewer,2);CHKERRQ(ierr);
777   }
778   ierr = PCView(fetidp->innerbddc,viewer);CHKERRQ(ierr);
779   if (iascii) {
780     ierr = PetscViewerASCIISubtractTab(viewer,2);CHKERRQ(ierr);
781   }
782   PetscFunctionReturn(0);
783 }
784 
785 #undef __FUNCT__
786 #define __FUNCT__ "KSPSetFromOptions_FETIDP"
787 static PetscErrorCode KSPSetFromOptions_FETIDP(PetscOptionItems *PetscOptionsObject,KSP ksp)
788 {
789   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
790   PetscErrorCode ierr;
791 
792   PetscFunctionBegin;
793   /* set options prefixes for the inner objects, since the parent prefix will be valid at this point */
794   ierr = PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerksp,((PetscObject)ksp)->prefix);CHKERRQ(ierr);
795   ierr = PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerksp,"fetidp_");CHKERRQ(ierr);
796   if (!fetidp->userbddc) {
797     ierr = PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerbddc,((PetscObject)ksp)->prefix);CHKERRQ(ierr);
798     ierr = PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerbddc,"fetidp_bddc_");CHKERRQ(ierr);
799   }
800   ierr = PetscOptionsHead(PetscOptionsObject,"KSP FETIDP options");CHKERRQ(ierr);
801   ierr = PetscOptionsBool("-ksp_fetidp_fullyredundant","Use fully redundant multipliers","none",fetidp->fully_redundant,&fetidp->fully_redundant,NULL);CHKERRQ(ierr);
802   ierr = PetscOptionsBool("-ksp_fetidp_saddlepoint","Activates support for saddle-point problems",NULL,fetidp->saddlepoint,&fetidp->saddlepoint,NULL);CHKERRQ(ierr);
803   ierr = PetscOptionsTail();CHKERRQ(ierr);
804   ierr = PCSetFromOptions(fetidp->innerbddc);CHKERRQ(ierr);
805   ierr = KSPSetFromOptions(fetidp->innerksp);CHKERRQ(ierr);
806   PetscFunctionReturn(0);
807 }
808 
809 /*MC
810      KSPFETIDP - The FETI-DP method
811 
812    This class implements the FETI-DP method [1].
813    The preconditioning matrix for the KSP must be of type MATIS.
814    The FETI-DP linear system (automatically generated constructing an internal PCBDDC object) is solved using an inner KSP object.
815 
816    Options Database Keys:
817 +   -ksp_fetidp_fullyredundant <false> : use a fully redundant set of Lagrange multipliers
818 -   -ksp_fetidp_saddlepoint <false>    : activates support for saddle point problems, see [2]
819 
820    Level: Advanced
821 
822    Notes: Options for the inner KSP and for the customization of the PCBDDC object can be specified at command line by using the prefixes -fetidp_ and -fetidp_bddc_. E.g.,
823 .vb
824       -fetidp_ksp_type gmres -fetidp_bddc_pc_bddc_symmetric false
825 .ve
826    will use GMRES for the solution of the linear system on the Lagrange multipliers, generated using a non-symmetric PCBDDC.
827    For saddle point problems with continuous pressures, the operators for the interface pressure solver can be specified with KSPFETIDPSetPressureOperators().
828 
829    References:
830 .vb
831 .  [1] - C. Farhat, M. Lesoinne, P. LeTallec, K. Pierson, and D. Rixen, FETI-DP: a dual-primal unified FETI method. I. A faster alternative to the two-level FETI method, Internat. J. Numer. Methods Engrg., 50 (2001), pp. 1523--1544
832 .  [2] - X. Tu, J. Li, A FETI-DP type domain decomposition algorithm for three-dimensional incompressible Stokes equations, SIAM J. Numer. Anal., 53 (2015), pp. 720-742
833 .ve
834 
835 .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC, KSPFETIDPGetInnerBDDC, KSPFETIDPGetInnerKSP
836 M*/
837 #undef __FUNCT__
838 #define __FUNCT__ "KSPCreate_FETIDP"
839 PETSC_EXTERN PetscErrorCode KSPCreate_FETIDP(KSP ksp)
840 {
841   PetscErrorCode ierr;
842   KSP_FETIDP     *fetidp;
843   KSP_FETIDPMon  *monctx;
844   PC_BDDC        *pcbddc;
845   PC             pc;
846 
847   PetscFunctionBegin;
848   ierr = PetscNewLog(ksp,&fetidp);CHKERRQ(ierr);
849   ksp->data = (void*)fetidp;
850   ksp->ops->setup                        = KSPSetUp_FETIDP;
851   ksp->ops->solve                        = KSPSolve_FETIDP;
852   ksp->ops->destroy                      = KSPDestroy_FETIDP;
853   ksp->ops->computeeigenvalues           = KSPComputeEigenvalues_FETIDP;
854   ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_FETIDP;
855   ksp->ops->view                         = KSPView_FETIDP;
856   ksp->ops->setfromoptions               = KSPSetFromOptions_FETIDP;
857   ksp->ops->buildsolution                = KSPBuildSolution_FETIDP;
858   ksp->ops->buildresidual                = KSPBuildResidualDefault;
859   /* create the inner KSP for the Lagrange multipliers */
860   ierr = KSPCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerksp);CHKERRQ(ierr);
861   ierr = KSPGetPC(fetidp->innerksp,&pc);CHKERRQ(ierr);
862   ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr);
863   ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerksp);CHKERRQ(ierr);
864   /* monitor */
865   ierr = PetscNew(&monctx);CHKERRQ(ierr);
866   monctx->parentksp = ksp;
867   fetidp->monctx = monctx;
868   ierr = KSPMonitorSet(fetidp->innerksp,KSPMonitor_FETIDP,fetidp->monctx,NULL);CHKERRQ(ierr);
869   /* create the inner BDDC */
870   ierr = PCCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerbddc);CHKERRQ(ierr);
871   ierr = PCSetType(fetidp->innerbddc,PCBDDC);CHKERRQ(ierr);
872   /* make sure we always obtain a consistent FETI-DP matrix
873      for symmetric problems, the user can always customize it through the command line */
874   pcbddc = (PC_BDDC*)fetidp->innerbddc->data;
875   pcbddc->symmetric_primal = PETSC_FALSE;
876   ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerbddc);CHKERRQ(ierr);
877   /* composed functions */
878   ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetInnerBDDC_C",KSPFETIDPSetInnerBDDC_FETIDP);CHKERRQ(ierr);
879   ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerBDDC_C",KSPFETIDPGetInnerBDDC_FETIDP);CHKERRQ(ierr);
880   ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerKSP_C",KSPFETIDPGetInnerKSP_FETIDP);CHKERRQ(ierr);
881   ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetPressureOperators_C",KSPFETIDPSetPressureOperators_FETIDP);CHKERRQ(ierr);
882   /* need to call KSPSetUp_FETIDP even with KSP_SETUP_NEWMATRIX */
883   ksp->setupnewmatrix = PETSC_TRUE;
884   PetscFunctionReturn(0);
885 }
886