1 #include <petsc/private/kspimpl.h> /*I <petscksp.h> I*/
2 #include <../src/ksp/pc/impls/bddc/bddc.h>
3 #include <../src/ksp/pc/impls/bddc/bddcprivate.h>
4 #include <petscdm.h>
5 
6 static PetscBool  cited  = PETSC_FALSE;
7 static PetscBool  cited2 = PETSC_FALSE;
8 static const char citation[] =
9 "@article{ZampiniPCBDDC,\n"
10 "author = {Stefano Zampini},\n"
11 "title = {{PCBDDC}: A Class of Robust Dual-Primal Methods in {PETS}c},\n"
12 "journal = {SIAM Journal on Scientific Computing},\n"
13 "volume = {38},\n"
14 "number = {5},\n"
15 "pages = {S282-S306},\n"
16 "year = {2016},\n"
17 "doi = {10.1137/15M1025785},\n"
18 "URL = {http://dx.doi.org/10.1137/15M1025785},\n"
19 "eprint = {http://dx.doi.org/10.1137/15M1025785}\n"
20 "}\n"
21 "@article{ZampiniDualPrimal,\n"
22 "author = {Stefano Zampini},\n"
23 "title = {{D}ual-{P}rimal methods for the cardiac {B}idomain model},\n"
24 "volume = {24},\n"
25 "number = {04},\n"
26 "pages = {667-696},\n"
27 "year = {2014},\n"
28 "doi = {10.1142/S0218202513500632},\n"
29 "URL = {https://www.worldscientific.com/doi/abs/10.1142/S0218202513500632},\n"
30 "eprint = {https://www.worldscientific.com/doi/pdf/10.1142/S0218202513500632}\n"
31 "}\n";
32 static const char citation2[] =
33 "@article{li2013nonoverlapping,\n"
34 "title={A nonoverlapping domain decomposition method for incompressible Stokes equations with continuous pressures},\n"
35 "author={Li, Jing and Tu, Xuemin},\n"
36 "journal={SIAM Journal on Numerical Analysis},\n"
37 "volume={51},\n"
38 "number={2},\n"
39 "pages={1235--1253},\n"
40 "year={2013},\n"
41 "publisher={Society for Industrial and Applied Mathematics}\n"
42 "}\n";
43 
44 /*
45     This file implements the FETI-DP method in PETSc as part of KSP.
46 */
47 typedef struct {
48   KSP parentksp;
49 } KSP_FETIDPMon;
50 
51 typedef struct {
52   KSP              innerksp;         /* the KSP for the Lagrange multipliers */
53   PC               innerbddc;        /* the inner BDDC object */
54   PetscBool        fully_redundant;  /* true for using a fully redundant set of multipliers */
55   PetscBool        userbddc;         /* true if the user provided the PCBDDC object */
56   PetscBool        saddlepoint;      /* support for saddle point problems */
57   IS               pP;               /* index set for pressure variables */
58   Vec              rhs_flip;         /* see KSPFETIDPSetUpOperators */
59   KSP_FETIDPMon    *monctx;          /* monitor context, used to pass user defined monitors
60                                         in the physical space */
61   PetscObjectState matstate;         /* these are needed just in the saddle point case */
62   PetscObjectState matnnzstate;      /* where we are going to use MatZeroRows on pmat */
63   PetscBool        statechanged;
64   PetscBool        check;
65 } KSP_FETIDP;
66 
KSPFETIDPSetPressureOperator_FETIDP(KSP ksp,Mat P)67 static PetscErrorCode KSPFETIDPSetPressureOperator_FETIDP(KSP ksp, Mat P)
68 {
69   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
70   PetscErrorCode ierr;
71 
72   PetscFunctionBegin;
73   if (P) fetidp->saddlepoint = PETSC_TRUE;
74   ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject)P);CHKERRQ(ierr);
75   PetscFunctionReturn(0);
76 }
77 
78 /*@
79  KSPFETIDPSetPressureOperator - Sets the operator used to setup the pressure preconditioner for saddle point FETI-DP.
80 
81    Collective on ksp
82 
83    Input Parameters:
84 +  ksp - the FETI-DP Krylov solver
85 -  P - the linear operator to be preconditioned, usually the mass matrix.
86 
87    Level: advanced
88 
89    Notes:
90     The operator can be either passed in a) monolithic global ordering, b) pressure-only global ordering
91           or c) interface pressure ordering (if -ksp_fetidp_pressure_all false).
92           In cases b) and c), the pressure ordering of dofs needs to satisfy
93              pid_1 < pid_2  iff  gid_1 < gid_2
94           where pid_1 and pid_2 are two different pressure dof numbers and gid_1 and gid_2 the corresponding
95           id in the monolithic global ordering.
96 
97 .seealso: MATIS, PCBDDC, KSPFETIDPGetInnerBDDC, KSPFETIDPGetInnerKSP, KSPSetOperators
98 @*/
KSPFETIDPSetPressureOperator(KSP ksp,Mat P)99 PetscErrorCode KSPFETIDPSetPressureOperator(KSP ksp, Mat P)
100 {
101   PetscErrorCode ierr;
102 
103   PetscFunctionBegin;
104   PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
105   if (P) PetscValidHeaderSpecific(P,MAT_CLASSID,2);
106   ierr = PetscTryMethod(ksp,"KSPFETIDPSetPressureOperator_C",(KSP,Mat),(ksp,P));CHKERRQ(ierr);
107   PetscFunctionReturn(0);
108 }
109 
KSPFETIDPGetInnerKSP_FETIDP(KSP ksp,KSP * innerksp)110 static PetscErrorCode KSPFETIDPGetInnerKSP_FETIDP(KSP ksp, KSP* innerksp)
111 {
112   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
113 
114   PetscFunctionBegin;
115   *innerksp = fetidp->innerksp;
116   PetscFunctionReturn(0);
117 }
118 
119 /*@
120  KSPFETIDPGetInnerKSP - Gets the KSP object for the Lagrange multipliers
121 
122    Input Parameters:
123 +  ksp - the FETI-DP KSP
124 -  innerksp - the KSP for the multipliers
125 
126    Level: advanced
127 
128    Notes:
129 
130 .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC, KSPFETIDPGetInnerBDDC
131 @*/
KSPFETIDPGetInnerKSP(KSP ksp,KSP * innerksp)132 PetscErrorCode KSPFETIDPGetInnerKSP(KSP ksp, KSP* innerksp)
133 {
134   PetscErrorCode ierr;
135 
136   PetscFunctionBegin;
137   PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
138   PetscValidPointer(innerksp,2);
139   ierr = PetscUseMethod(ksp,"KSPFETIDPGetInnerKSP_C",(KSP,KSP*),(ksp,innerksp));CHKERRQ(ierr);
140   PetscFunctionReturn(0);
141 }
142 
KSPFETIDPGetInnerBDDC_FETIDP(KSP ksp,PC * pc)143 static PetscErrorCode KSPFETIDPGetInnerBDDC_FETIDP(KSP ksp, PC* pc)
144 {
145   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
146 
147   PetscFunctionBegin;
148   *pc = fetidp->innerbddc;
149   PetscFunctionReturn(0);
150 }
151 
152 /*@
153  KSPFETIDPGetInnerBDDC - Gets the BDDC preconditioner used to setup the FETI-DP matrix for the Lagrange multipliers
154 
155    Input Parameters:
156 +  ksp - the FETI-DP Krylov solver
157 -  pc - the BDDC preconditioner
158 
159    Level: advanced
160 
161    Notes:
162 
163 .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC, KSPFETIDPGetInnerKSP
164 @*/
KSPFETIDPGetInnerBDDC(KSP ksp,PC * pc)165 PetscErrorCode KSPFETIDPGetInnerBDDC(KSP ksp, PC* pc)
166 {
167   PetscErrorCode ierr;
168 
169   PetscFunctionBegin;
170   PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
171   PetscValidPointer(pc,2);
172   ierr = PetscUseMethod(ksp,"KSPFETIDPGetInnerBDDC_C",(KSP,PC*),(ksp,pc));CHKERRQ(ierr);
173   PetscFunctionReturn(0);
174 }
175 
KSPFETIDPSetInnerBDDC_FETIDP(KSP ksp,PC pc)176 static PetscErrorCode KSPFETIDPSetInnerBDDC_FETIDP(KSP ksp, PC pc)
177 {
178   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
179   PetscErrorCode ierr;
180 
181   PetscFunctionBegin;
182   ierr = PetscObjectReference((PetscObject)pc);CHKERRQ(ierr);
183   ierr = PCDestroy(&fetidp->innerbddc);CHKERRQ(ierr);
184   fetidp->innerbddc = pc;
185   fetidp->userbddc  = PETSC_TRUE;
186   PetscFunctionReturn(0);
187 }
188 
189 /*@
190  KSPFETIDPSetInnerBDDC - Sets the BDDC preconditioner used to setup the FETI-DP matrix for the Lagrange multipliers
191 
192    Collective on ksp
193 
194    Input Parameters:
195 +  ksp - the FETI-DP Krylov solver
196 -  pc - the BDDC preconditioner
197 
198    Level: advanced
199 
200    Notes:
201 
202 .seealso: MATIS, PCBDDC, KSPFETIDPGetInnerBDDC, KSPFETIDPGetInnerKSP
203 @*/
KSPFETIDPSetInnerBDDC(KSP ksp,PC pc)204 PetscErrorCode KSPFETIDPSetInnerBDDC(KSP ksp, PC pc)
205 {
206   PetscBool      isbddc;
207   PetscErrorCode ierr;
208 
209   PetscFunctionBegin;
210   PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
211   PetscValidHeaderSpecific(pc,PC_CLASSID,2);
212   ierr = PetscObjectTypeCompare((PetscObject)pc,PCBDDC,&isbddc);CHKERRQ(ierr);
213   if (!isbddc) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONG,"KSPFETIDPSetInnerBDDC need a PCBDDC preconditioner");
214   ierr = PetscTryMethod(ksp,"KSPFETIDPSetInnerBDDC_C",(KSP,PC),(ksp,pc));CHKERRQ(ierr);
215   PetscFunctionReturn(0);
216 }
217 
KSPBuildSolution_FETIDP(KSP ksp,Vec v,Vec * V)218 static PetscErrorCode KSPBuildSolution_FETIDP(KSP ksp,Vec v,Vec *V)
219 {
220   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
221   Mat            F;
222   Vec            Xl;
223   PetscErrorCode ierr;
224 
225   PetscFunctionBegin;
226   ierr = KSPGetOperators(fetidp->innerksp,&F,NULL);CHKERRQ(ierr);
227   ierr = KSPBuildSolution(fetidp->innerksp,NULL,&Xl);CHKERRQ(ierr);
228   if (v) {
229     ierr = PCBDDCMatFETIDPGetSolution(F,Xl,v);CHKERRQ(ierr);
230     *V   = v;
231   } else {
232     ierr = PCBDDCMatFETIDPGetSolution(F,Xl,*V);CHKERRQ(ierr);
233   }
234   PetscFunctionReturn(0);
235 }
236 
KSPMonitor_FETIDP(KSP ksp,PetscInt it,PetscReal rnorm,void * ctx)237 static PetscErrorCode KSPMonitor_FETIDP(KSP ksp,PetscInt it,PetscReal rnorm,void* ctx)
238 {
239   KSP_FETIDPMon  *monctx = (KSP_FETIDPMon*)ctx;
240   PetscErrorCode ierr;
241 
242   PetscFunctionBegin;
243   ierr = KSPMonitor(monctx->parentksp,it,rnorm);CHKERRQ(ierr);
244   PetscFunctionReturn(0);
245 }
246 
KSPComputeEigenvalues_FETIDP(KSP ksp,PetscInt nmax,PetscReal * r,PetscReal * c,PetscInt * neig)247 static PetscErrorCode KSPComputeEigenvalues_FETIDP(KSP ksp,PetscInt nmax,PetscReal *r,PetscReal *c,PetscInt *neig)
248 {
249   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
250   PetscErrorCode ierr;
251 
252   PetscFunctionBegin;
253   ierr = KSPComputeEigenvalues(fetidp->innerksp,nmax,r,c,neig);CHKERRQ(ierr);
254   PetscFunctionReturn(0);
255 }
256 
KSPComputeExtremeSingularValues_FETIDP(KSP ksp,PetscReal * emax,PetscReal * emin)257 static PetscErrorCode KSPComputeExtremeSingularValues_FETIDP(KSP ksp,PetscReal *emax,PetscReal *emin)
258 {
259   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
260   PetscErrorCode ierr;
261 
262   PetscFunctionBegin;
263   ierr = KSPComputeExtremeSingularValues(fetidp->innerksp,emax,emin);CHKERRQ(ierr);
264   PetscFunctionReturn(0);
265 }
266 
KSPFETIDPCheckOperators(KSP ksp,PetscViewer viewer)267 static PetscErrorCode KSPFETIDPCheckOperators(KSP ksp, PetscViewer viewer)
268 {
269   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
270   PC_BDDC        *pcbddc = (PC_BDDC*)fetidp->innerbddc->data;
271   PC_IS          *pcis = (PC_IS*)fetidp->innerbddc->data;
272   Mat_IS         *matis = (Mat_IS*)fetidp->innerbddc->pmat->data;
273   Mat            F;
274   FETIDPMat_ctx  fetidpmat_ctx;
275   Vec            test_vec,test_vec_p = NULL,fetidp_global;
276   IS             dirdofs,isvert;
277   MPI_Comm       comm = PetscObjectComm((PetscObject)ksp);
278   PetscScalar    sval,*array;
279   PetscReal      val,rval;
280   const PetscInt *vertex_indices;
281   PetscInt       i,n_vertices;
282   PetscBool      isascii;
283   PetscErrorCode ierr;
284 
285   PetscFunctionBegin;
286   PetscCheckSameComm(ksp,1,viewer,2);
287   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
288   if (!isascii) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported viewer");
289   ierr = PetscViewerASCIIPrintf(viewer,"----------FETI-DP MAT  --------------\n");CHKERRQ(ierr);
290   ierr = PetscViewerASCIIAddTab(viewer,2);CHKERRQ(ierr);
291   ierr = KSPGetOperators(fetidp->innerksp,&F,NULL);CHKERRQ(ierr);
292   ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
293   ierr = MatView(F,viewer);CHKERRQ(ierr);
294   ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
295   ierr = PetscViewerASCIISubtractTab(viewer,2);CHKERRQ(ierr);
296   ierr = MatShellGetContext(F,(void**)&fetidpmat_ctx);CHKERRQ(ierr);
297   ierr = PetscViewerASCIIPrintf(viewer,"----------FETI-DP TESTS--------------\n");CHKERRQ(ierr);
298   ierr = PetscViewerASCIIPrintf(viewer,"All tests should return zero!\n");CHKERRQ(ierr);
299   ierr = PetscViewerASCIIPrintf(viewer,"FETIDP MAT context in the ");CHKERRQ(ierr);
300   if (fetidp->fully_redundant) {
301     ierr = PetscViewerASCIIPrintf(viewer,"fully redundant case for lagrange multipliers.\n");CHKERRQ(ierr);
302   } else {
303     ierr = PetscViewerASCIIPrintf(viewer,"Non-fully redundant case for lagrange multiplier.\n");CHKERRQ(ierr);
304   }
305   ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
306 
307   /* Get Vertices used to define the BDDC */
308   ierr = PCBDDCGraphGetCandidatesIS(pcbddc->mat_graph,NULL,NULL,NULL,NULL,&isvert);CHKERRQ(ierr);
309   ierr = ISGetLocalSize(isvert,&n_vertices);CHKERRQ(ierr);
310   ierr = ISGetIndices(isvert,&vertex_indices);CHKERRQ(ierr);
311 
312   /******************************************************************/
313   /* TEST A/B: Test numbering of global fetidp dofs                 */
314   /******************************************************************/
315   ierr = MatCreateVecs(F,&fetidp_global,NULL);CHKERRQ(ierr);
316   ierr = VecDuplicate(fetidpmat_ctx->lambda_local,&test_vec);CHKERRQ(ierr);
317   ierr = VecSet(fetidp_global,1.0);CHKERRQ(ierr);
318   ierr = VecSet(test_vec,1.);CHKERRQ(ierr);
319   ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
320   ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
321   if (fetidpmat_ctx->l2g_p) {
322     ierr = VecDuplicate(fetidpmat_ctx->vP,&test_vec_p);CHKERRQ(ierr);
323     ierr = VecSet(test_vec_p,1.);CHKERRQ(ierr);
324     ierr = VecScatterBegin(fetidpmat_ctx->l2g_p,fetidp_global,fetidpmat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
325     ierr = VecScatterEnd(fetidpmat_ctx->l2g_p,fetidp_global,fetidpmat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
326   }
327   ierr = VecAXPY(test_vec,-1.0,fetidpmat_ctx->lambda_local);CHKERRQ(ierr);
328   ierr = VecNorm(test_vec,NORM_INFINITY,&val);CHKERRQ(ierr);
329   ierr = VecDestroy(&test_vec);CHKERRQ(ierr);
330   ierr = MPI_Reduce(&val,&rval,1,MPIU_REAL,MPI_MAX,0,comm);CHKERRQ(ierr);
331   ierr = PetscViewerASCIIPrintf(viewer,"A: CHECK glob to loc: % 1.14e\n",rval);CHKERRQ(ierr);
332 
333   if (fetidpmat_ctx->l2g_p) {
334     ierr = VecAXPY(test_vec_p,-1.0,fetidpmat_ctx->vP);CHKERRQ(ierr);
335     ierr = VecNorm(test_vec_p,NORM_INFINITY,&val);CHKERRQ(ierr);
336     ierr = MPI_Reduce(&val,&rval,1,MPIU_REAL,MPI_MAX,0,comm);CHKERRQ(ierr);
337     ierr = PetscViewerASCIIPrintf(viewer,"A: CHECK glob to loc (p): % 1.14e\n",rval);CHKERRQ(ierr);
338   }
339 
340   if (fetidp->fully_redundant) {
341     ierr = VecSet(fetidp_global,0.0);CHKERRQ(ierr);
342     ierr = VecSet(fetidpmat_ctx->lambda_local,0.5);CHKERRQ(ierr);
343     ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
344     ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
345     ierr = VecSum(fetidp_global,&sval);CHKERRQ(ierr);
346     val  = PetscRealPart(sval)-fetidpmat_ctx->n_lambda;
347     ierr = MPI_Reduce(&val,&rval,1,MPIU_REAL,MPI_MAX,0,comm);CHKERRQ(ierr);
348     ierr = PetscViewerASCIIPrintf(viewer,"B: CHECK loc to glob: % 1.14e\n",rval);CHKERRQ(ierr);
349   }
350 
351   if (fetidpmat_ctx->l2g_p) {
352     ierr = VecSet(pcis->vec1_N,1.0);CHKERRQ(ierr);
353     ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr);
354     ierr = VecScatterBegin(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
355     ierr = VecScatterEnd(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
356 
357     ierr = VecSet(fetidp_global,0.0);CHKERRQ(ierr);
358     ierr = VecSet(fetidpmat_ctx->vP,-1.0);CHKERRQ(ierr);
359     ierr = VecScatterBegin(fetidpmat_ctx->l2g_p,fetidpmat_ctx->vP,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
360     ierr = VecScatterEnd(fetidpmat_ctx->l2g_p,fetidpmat_ctx->vP,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
361     ierr = VecScatterBegin(fetidpmat_ctx->g2g_p,fetidp_global,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
362     ierr = VecScatterEnd(fetidpmat_ctx->g2g_p,fetidp_global,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
363     ierr = VecScatterBegin(fetidpmat_ctx->g2g_p,pcis->vec1_global,fetidp_global,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
364     ierr = VecScatterEnd(fetidpmat_ctx->g2g_p,pcis->vec1_global,fetidp_global,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
365     ierr = VecSum(fetidp_global,&sval);CHKERRQ(ierr);
366     val  = PetscRealPart(sval);
367     ierr = MPI_Reduce(&val,&rval,1,MPIU_REAL,MPI_MAX,0,comm);CHKERRQ(ierr);
368     ierr = PetscViewerASCIIPrintf(viewer,"B: CHECK loc to glob (p): % 1.14e\n",rval);CHKERRQ(ierr);
369   }
370 
371   /******************************************************************/
372   /* TEST C: It should hold B_delta*w=0, w\in\widehat{W}            */
373   /* This is the meaning of the B matrix                            */
374   /******************************************************************/
375 
376   ierr = VecSetRandom(pcis->vec1_N,NULL);CHKERRQ(ierr);
377   ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr);
378   ierr = VecScatterBegin(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
379   ierr = VecScatterEnd(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
380   ierr = VecScatterBegin(matis->rctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
381   ierr = VecScatterEnd(matis->rctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
382   ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
383   ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
384   /* Action of B_delta */
385   ierr = MatMult(fetidpmat_ctx->B_delta,pcis->vec1_B,fetidpmat_ctx->lambda_local);CHKERRQ(ierr);
386   ierr = VecSet(fetidp_global,0.0);CHKERRQ(ierr);
387   ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
388   ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
389   ierr = VecNorm(fetidp_global,NORM_INFINITY,&val);CHKERRQ(ierr);
390   ierr = PetscViewerASCIIPrintf(viewer,"C: CHECK infty norm of B_delta*w (w continuous): % 1.14e\n",val);CHKERRQ(ierr);
391 
392   /******************************************************************/
393   /* TEST D: It should hold E_Dw = w - P_Dw w\in\widetilde{W}       */
394   /* E_D = R_D^TR                                                   */
395   /* P_D = B_{D,delta}^T B_{delta}                                  */
396   /* eq.44 Mandel Tezaur and Dohrmann 2005                          */
397   /******************************************************************/
398 
399   /* compute a random vector in \widetilde{W} */
400   ierr = VecSetRandom(pcis->vec1_N,NULL);CHKERRQ(ierr);
401   /* set zero at vertices and essential dofs */
402   ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr);
403   for (i=0;i<n_vertices;i++) array[vertex_indices[i]] = 0.0;
404   ierr = PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph,&dirdofs);CHKERRQ(ierr);
405   if (dirdofs) {
406     const PetscInt *idxs;
407     PetscInt       ndir;
408 
409     ierr = ISGetLocalSize(dirdofs,&ndir);CHKERRQ(ierr);
410     ierr = ISGetIndices(dirdofs,&idxs);CHKERRQ(ierr);
411     for (i=0;i<ndir;i++) array[idxs[i]] = 0.0;
412     ierr = ISRestoreIndices(dirdofs,&idxs);CHKERRQ(ierr);
413   }
414   ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr);
415   /* store w for final comparison */
416   ierr = VecDuplicate(pcis->vec1_B,&test_vec);CHKERRQ(ierr);
417   ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,test_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
418   ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_N,test_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
419 
420   /* Jump operator P_D : results stored in pcis->vec1_B */
421   /* Action of B_delta */
422   ierr = MatMult(fetidpmat_ctx->B_delta,test_vec,fetidpmat_ctx->lambda_local);CHKERRQ(ierr);
423   ierr = VecSet(fetidp_global,0.0);CHKERRQ(ierr);
424   ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
425   ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
426   /* Action of B_Ddelta^T */
427   ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
428   ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
429   ierr = MatMultTranspose(fetidpmat_ctx->B_Ddelta,fetidpmat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr);
430 
431   /* Average operator E_D : results stored in pcis->vec2_B */
432   ierr = PCBDDCScalingExtension(fetidpmat_ctx->pc,test_vec,pcis->vec1_global);CHKERRQ(ierr);
433   ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec2_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
434   ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec2_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
435 
436   /* test E_D=I-P_D */
437   ierr = VecAXPY(pcis->vec1_B,1.0,pcis->vec2_B);CHKERRQ(ierr);
438   ierr = VecAXPY(pcis->vec1_B,-1.0,test_vec);CHKERRQ(ierr);
439   ierr = VecNorm(pcis->vec1_B,NORM_INFINITY,&val);CHKERRQ(ierr);
440   ierr = VecDestroy(&test_vec);CHKERRQ(ierr);
441   ierr = MPI_Reduce(&val,&rval,1,MPIU_REAL,MPI_MAX,0,comm);CHKERRQ(ierr);
442   ierr = PetscViewerASCIIPrintf(viewer,"D: CHECK infty norm of E_D + P_D - I: % 1.14e\n",PetscGlobalRank,val);CHKERRQ(ierr);
443 
444   /******************************************************************/
445   /* TEST E: It should hold R_D^TP_Dw=0 w\in\widetilde{W}           */
446   /* eq.48 Mandel Tezaur and Dohrmann 2005                          */
447   /******************************************************************/
448 
449   ierr = VecSetRandom(pcis->vec1_N,NULL);CHKERRQ(ierr);
450   /* set zero at vertices and essential dofs */
451   ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr);
452   for (i=0;i<n_vertices;i++) array[vertex_indices[i]] = 0.0;
453   if (dirdofs) {
454     const PetscInt *idxs;
455     PetscInt       ndir;
456 
457     ierr = ISGetLocalSize(dirdofs,&ndir);CHKERRQ(ierr);
458     ierr = ISGetIndices(dirdofs,&idxs);CHKERRQ(ierr);
459     for (i=0;i<ndir;i++) array[idxs[i]] = 0.0;
460     ierr = ISRestoreIndices(dirdofs,&idxs);CHKERRQ(ierr);
461   }
462   ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr);
463 
464   /* Jump operator P_D : results stored in pcis->vec1_B */
465 
466   ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
467   ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
468   /* Action of B_delta */
469   ierr = MatMult(fetidpmat_ctx->B_delta,pcis->vec1_B,fetidpmat_ctx->lambda_local);CHKERRQ(ierr);
470   ierr = VecSet(fetidp_global,0.0);CHKERRQ(ierr);
471   ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
472   ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
473   /* Action of B_Ddelta^T */
474   ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
475   ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
476   ierr = MatMultTranspose(fetidpmat_ctx->B_Ddelta,fetidpmat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr);
477   /* scaling */
478   ierr = PCBDDCScalingExtension(fetidpmat_ctx->pc,pcis->vec1_B,pcis->vec1_global);CHKERRQ(ierr);
479   ierr = VecNorm(pcis->vec1_global,NORM_INFINITY,&val);CHKERRQ(ierr);
480   ierr = PetscViewerASCIIPrintf(viewer,"E: CHECK infty norm of R^T_D P_D: % 1.14e\n",val);CHKERRQ(ierr);
481 
482   if (!fetidp->fully_redundant) {
483     /******************************************************************/
484     /* TEST F: It should holds B_{delta}B^T_{D,delta}=I               */
485     /* Corollary thm 14 Mandel Tezaur and Dohrmann 2005               */
486     /******************************************************************/
487     ierr = VecDuplicate(fetidp_global,&test_vec);CHKERRQ(ierr);
488     ierr = VecSetRandom(fetidp_global,NULL);CHKERRQ(ierr);
489     if (fetidpmat_ctx->l2g_p) {
490       ierr = VecSet(fetidpmat_ctx->vP,0.);CHKERRQ(ierr);
491       ierr = VecScatterBegin(fetidpmat_ctx->l2g_p,fetidpmat_ctx->vP,fetidp_global,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
492       ierr = VecScatterEnd(fetidpmat_ctx->l2g_p,fetidpmat_ctx->vP,fetidp_global,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
493     }
494     /* Action of B_Ddelta^T */
495     ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
496     ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
497     ierr = MatMultTranspose(fetidpmat_ctx->B_Ddelta,fetidpmat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr);
498     /* Action of B_delta */
499     ierr = MatMult(fetidpmat_ctx->B_delta,pcis->vec1_B,fetidpmat_ctx->lambda_local);CHKERRQ(ierr);
500     ierr = VecSet(test_vec,0.0);CHKERRQ(ierr);
501     ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,test_vec,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
502     ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,test_vec,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
503     ierr = VecAXPY(fetidp_global,-1.,test_vec);CHKERRQ(ierr);
504     ierr = VecNorm(fetidp_global,NORM_INFINITY,&val);CHKERRQ(ierr);
505     ierr = PetscViewerASCIIPrintf(viewer,"E: CHECK infty norm of P^T_D - I: % 1.14e\n",val);CHKERRQ(ierr);
506     ierr = VecDestroy(&test_vec);CHKERRQ(ierr);
507   }
508   ierr = PetscViewerASCIIPrintf(viewer,"-------------------------------------\n");CHKERRQ(ierr);
509   ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
510   ierr = VecDestroy(&test_vec_p);CHKERRQ(ierr);
511   ierr = ISDestroy(&dirdofs);CHKERRQ(ierr);
512   ierr = VecDestroy(&fetidp_global);CHKERRQ(ierr);
513   ierr = ISRestoreIndices(isvert,&vertex_indices);CHKERRQ(ierr);
514   ierr = PCBDDCGraphRestoreCandidatesIS(pcbddc->mat_graph,NULL,NULL,NULL,NULL,&isvert);CHKERRQ(ierr);
515   PetscFunctionReturn(0);
516 }
517 
KSPFETIDPSetUpOperators(KSP ksp)518 static PetscErrorCode KSPFETIDPSetUpOperators(KSP ksp)
519 {
520   KSP_FETIDP       *fetidp = (KSP_FETIDP*)ksp->data;
521   PC_BDDC          *pcbddc = (PC_BDDC*)fetidp->innerbddc->data;
522   Mat              A,Ap;
523   PetscInt         fid = -1;
524   PetscMPIInt      size;
525   PetscBool        ismatis,pisz,allp,schp;
526   PetscBool        flip; /* Usually, Stokes is written (B = -\int_\Omega \nabla \cdot u q)
527                            | A B'| | v | = | f |
528                            | B 0 | | p | = | g |
529                             If -ksp_fetidp_saddlepoint_flip is true, the code assumes it is written as
530                            | A B'| | v | = | f |
531                            |-B 0 | | p | = |-g |
532                          */
533   PetscObjectState matstate, matnnzstate;
534   PetscErrorCode   ierr;
535 
536   PetscFunctionBegin;
537   pisz = PETSC_FALSE;
538   flip = PETSC_FALSE;
539   allp = PETSC_FALSE;
540   schp = PETSC_FALSE;
541   ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->prefix,"FETI-DP options","PC");CHKERRQ(ierr);
542   ierr = PetscOptionsInt("-ksp_fetidp_pressure_field","Field id for pressures for saddle-point problems",NULL,fid,&fid,NULL);CHKERRQ(ierr);
543   ierr = PetscOptionsBool("-ksp_fetidp_pressure_all","Use the whole pressure set instead of just that at the interface",NULL,allp,&allp,NULL);CHKERRQ(ierr);
544   ierr = PetscOptionsBool("-ksp_fetidp_saddlepoint_flip","Flip the sign of the pressure-velocity (lower-left) block",NULL,flip,&flip,NULL);CHKERRQ(ierr);
545   ierr = PetscOptionsBool("-ksp_fetidp_pressure_schur","Use a BDDC solver for pressure",NULL,schp,&schp,NULL);CHKERRQ(ierr);
546   ierr = PetscOptionsEnd();CHKERRQ(ierr);
547 
548   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)ksp),&size);CHKERRQ(ierr);
549   fetidp->saddlepoint = (fid >= 0 ? PETSC_TRUE : fetidp->saddlepoint);
550   if (size == 1) fetidp->saddlepoint = PETSC_FALSE;
551 
552   ierr = KSPGetOperators(ksp,&A,&Ap);CHKERRQ(ierr);
553   ierr = PetscObjectTypeCompare((PetscObject)A,MATIS,&ismatis);CHKERRQ(ierr);
554   if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Amat should be of type MATIS");
555 
556   /* Quiet return if the matrix states are unchanged.
557      Needed only for the saddle point case since it uses MatZeroRows
558      on a matrix that may not have changed */
559   ierr = PetscObjectStateGet((PetscObject)A,&matstate);CHKERRQ(ierr);
560   ierr = MatGetNonzeroState(A,&matnnzstate);CHKERRQ(ierr);
561   if (matstate == fetidp->matstate && matnnzstate == fetidp->matnnzstate) PetscFunctionReturn(0);
562   fetidp->matstate     = matstate;
563   fetidp->matnnzstate  = matnnzstate;
564   fetidp->statechanged = fetidp->saddlepoint;
565 
566   /* see if we have some fields attached */
567   if (!pcbddc->n_ISForDofsLocal && !pcbddc->n_ISForDofs) {
568     DM             dm;
569     PetscContainer c;
570 
571     ierr = KSPGetDM(ksp,&dm);CHKERRQ(ierr);
572     ierr = PetscObjectQuery((PetscObject)A,"_convert_nest_lfields",(PetscObject*)&c);CHKERRQ(ierr);
573     if (dm) {
574       IS      *fields;
575       PetscInt nf,i;
576 
577       ierr = DMCreateFieldDecomposition(dm,&nf,NULL,&fields,NULL);CHKERRQ(ierr);
578       ierr = PCBDDCSetDofsSplitting(fetidp->innerbddc,nf,fields);CHKERRQ(ierr);
579       for (i=0;i<nf;i++) {
580         ierr = ISDestroy(&fields[i]);CHKERRQ(ierr);
581       }
582       ierr = PetscFree(fields);CHKERRQ(ierr);
583     } else if (c) {
584       MatISLocalFields lf;
585 
586       ierr = PetscContainerGetPointer(c,(void**)&lf);CHKERRQ(ierr);
587       ierr = PCBDDCSetDofsSplittingLocal(fetidp->innerbddc,lf->nr,lf->rf);CHKERRQ(ierr);
588     }
589   }
590 
591   if (!fetidp->saddlepoint) {
592     ierr = PCSetOperators(fetidp->innerbddc,A,A);CHKERRQ(ierr);
593   } else {
594     Mat          nA,lA,PPmat;
595     MatNullSpace nnsp;
596     IS           pP;
597     PetscInt     totP;
598 
599     ierr = MatISGetLocalMat(A,&lA);CHKERRQ(ierr);
600     ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",(PetscObject)lA);CHKERRQ(ierr);
601 
602     pP = fetidp->pP;
603     if (!pP) { /* first time, need to compute pressure dofs */
604       PC_IS                  *pcis = (PC_IS*)fetidp->innerbddc->data;
605       Mat_IS                 *matis = (Mat_IS*)(A->data);
606       ISLocalToGlobalMapping l2g;
607       IS                     lP = NULL,II,pII,lPall,Pall,is1,is2;
608       const PetscInt         *idxs;
609       PetscInt               nl,ni,*widxs;
610       PetscInt               i,j,n_neigh,*neigh,*n_shared,**shared,*count;
611       PetscInt               rst,ren,n;
612       PetscBool              ploc;
613 
614       ierr = MatGetLocalSize(A,&nl,NULL);CHKERRQ(ierr);
615       ierr = MatGetOwnershipRange(A,&rst,&ren);CHKERRQ(ierr);
616       ierr = MatGetLocalSize(lA,&n,NULL);CHKERRQ(ierr);
617       ierr = MatGetLocalToGlobalMapping(A,&l2g,NULL);CHKERRQ(ierr);
618 
619       if (!pcis->is_I_local) { /* need to compute interior dofs */
620         ierr = PetscCalloc1(n,&count);CHKERRQ(ierr);
621         ierr = ISLocalToGlobalMappingGetInfo(l2g,&n_neigh,&neigh,&n_shared,&shared);CHKERRQ(ierr);
622         for (i=1;i<n_neigh;i++)
623           for (j=0;j<n_shared[i];j++)
624             count[shared[i][j]] += 1;
625         for (i=0,j=0;i<n;i++) if (!count[i]) count[j++] = i;
626         ierr = ISLocalToGlobalMappingRestoreInfo(l2g,&n_neigh,&neigh,&n_shared,&shared);CHKERRQ(ierr);
627         ierr = ISCreateGeneral(PETSC_COMM_SELF,j,count,PETSC_OWN_POINTER,&II);CHKERRQ(ierr);
628       } else {
629         ierr = PetscObjectReference((PetscObject)pcis->is_I_local);CHKERRQ(ierr);
630         II   = pcis->is_I_local;
631       }
632 
633       /* interior dofs in layout */
634       ierr = PetscArrayzero(matis->sf_leafdata,n);CHKERRQ(ierr);
635       ierr = PetscArrayzero(matis->sf_rootdata,nl);CHKERRQ(ierr);
636       ierr = ISGetLocalSize(II,&ni);CHKERRQ(ierr);
637       ierr = ISGetIndices(II,&idxs);CHKERRQ(ierr);
638       for (i=0;i<ni;i++) matis->sf_leafdata[idxs[i]] = 1;
639       ierr = ISRestoreIndices(II,&idxs);CHKERRQ(ierr);
640       ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr);
641       ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr);
642       ierr = PetscMalloc1(PetscMax(nl,n),&widxs);CHKERRQ(ierr);
643       for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst;
644       ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&pII);CHKERRQ(ierr);
645 
646       /* pressure dofs */
647       Pall  = NULL;
648       lPall = NULL;
649       ploc  = PETSC_FALSE;
650       if (fid < 0) { /* zero pressure block */
651         PetscInt np;
652 
653         ierr = MatFindZeroDiagonals(A,&Pall);CHKERRQ(ierr);
654         ierr = ISGetSize(Pall,&np);CHKERRQ(ierr);
655         if (!np) { /* zero-block not found, defaults to last field (if set) */
656           fid  = pcbddc->n_ISForDofsLocal ? pcbddc->n_ISForDofsLocal - 1 : pcbddc->n_ISForDofs - 1;
657           ierr = ISDestroy(&Pall);CHKERRQ(ierr);
658         } else if (!pcbddc->n_ISForDofsLocal && !pcbddc->n_ISForDofs) {
659           ierr = PCBDDCSetDofsSplitting(fetidp->innerbddc,1,&Pall);CHKERRQ(ierr);
660         }
661       }
662       if (!Pall) { /* look for registered fields */
663         if (pcbddc->n_ISForDofsLocal) {
664           PetscInt np;
665 
666           if (fid < 0 || fid >= pcbddc->n_ISForDofsLocal) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Invalid field id for pressure %D, max %D",fid,pcbddc->n_ISForDofsLocal);
667           /* need a sequential IS */
668           ierr = ISGetLocalSize(pcbddc->ISForDofsLocal[fid],&np);CHKERRQ(ierr);
669           ierr = ISGetIndices(pcbddc->ISForDofsLocal[fid],&idxs);CHKERRQ(ierr);
670           ierr = ISCreateGeneral(PETSC_COMM_SELF,np,idxs,PETSC_COPY_VALUES,&lPall);CHKERRQ(ierr);
671           ierr = ISRestoreIndices(pcbddc->ISForDofsLocal[fid],&idxs);CHKERRQ(ierr);
672           ploc = PETSC_TRUE;
673         } else if (pcbddc->n_ISForDofs) {
674           if (fid < 0 || fid >= pcbddc->n_ISForDofs) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Invalid field id for pressure %D, max %D",fid,pcbddc->n_ISForDofs);
675           ierr = PetscObjectReference((PetscObject)pcbddc->ISForDofs[fid]);CHKERRQ(ierr);
676           Pall = pcbddc->ISForDofs[fid];
677         } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Cannot detect pressure field! Use KSPFETIDPGetInnerBDDC() + PCBDDCSetDofsSplitting or PCBDDCSetDofsSplittingLocal");
678       }
679 
680       /* if the user requested the entire pressure,
681          remove the interior pressure dofs from II (or pII) */
682       if (allp) {
683         if (ploc) {
684           IS nII;
685           ierr = ISDifference(II,lPall,&nII);CHKERRQ(ierr);
686           ierr = ISDestroy(&II);CHKERRQ(ierr);
687           II   = nII;
688         } else {
689           IS nII;
690           ierr = ISDifference(pII,Pall,&nII);CHKERRQ(ierr);
691           ierr = ISDestroy(&pII);CHKERRQ(ierr);
692           pII  = nII;
693         }
694       }
695       if (ploc) {
696         ierr = ISDifference(lPall,II,&lP);CHKERRQ(ierr);
697         ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lP",(PetscObject)lP);CHKERRQ(ierr);
698       } else {
699         ierr = ISDifference(Pall,pII,&pP);CHKERRQ(ierr);
700         ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",(PetscObject)pP);CHKERRQ(ierr);
701         /* need all local pressure dofs */
702         ierr = PetscArrayzero(matis->sf_leafdata,n);CHKERRQ(ierr);
703         ierr = PetscArrayzero(matis->sf_rootdata,nl);CHKERRQ(ierr);
704         ierr = ISGetLocalSize(Pall,&ni);CHKERRQ(ierr);
705         ierr = ISGetIndices(Pall,&idxs);CHKERRQ(ierr);
706         for (i=0;i<ni;i++) matis->sf_rootdata[idxs[i]-rst] = 1;
707         ierr = ISRestoreIndices(Pall,&idxs);CHKERRQ(ierr);
708         ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr);
709         ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr);
710         for (i=0,ni=0;i<n;i++) if (matis->sf_leafdata[i]) widxs[ni++] = i;
711         ierr = ISCreateGeneral(PETSC_COMM_SELF,ni,widxs,PETSC_COPY_VALUES,&lPall);CHKERRQ(ierr);
712       }
713 
714       if (!Pall) {
715         ierr = PetscArrayzero(matis->sf_leafdata,n);CHKERRQ(ierr);
716         ierr = PetscArrayzero(matis->sf_rootdata,nl);CHKERRQ(ierr);
717         ierr = ISGetLocalSize(lPall,&ni);CHKERRQ(ierr);
718         ierr = ISGetIndices(lPall,&idxs);CHKERRQ(ierr);
719         for (i=0;i<ni;i++) matis->sf_leafdata[idxs[i]] = 1;
720         ierr = ISRestoreIndices(lPall,&idxs);CHKERRQ(ierr);
721         ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr);
722         ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr);
723         for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst;
724         ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&Pall);CHKERRQ(ierr);
725       }
726       ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject)Pall);CHKERRQ(ierr);
727 
728       if (flip) {
729         PetscInt npl;
730         ierr = ISGetLocalSize(Pall,&npl);CHKERRQ(ierr);
731         ierr = ISGetIndices(Pall,&idxs);CHKERRQ(ierr);
732         ierr = MatCreateVecs(A,NULL,&fetidp->rhs_flip);CHKERRQ(ierr);
733         ierr = VecSet(fetidp->rhs_flip,1.);CHKERRQ(ierr);
734         ierr = VecSetOption(fetidp->rhs_flip,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
735         for (i=0;i<npl;i++) {
736           ierr = VecSetValue(fetidp->rhs_flip,idxs[i],-1.,INSERT_VALUES);CHKERRQ(ierr);
737         }
738         ierr = VecAssemblyBegin(fetidp->rhs_flip);CHKERRQ(ierr);
739         ierr = VecAssemblyEnd(fetidp->rhs_flip);CHKERRQ(ierr);
740         ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_flip",(PetscObject)fetidp->rhs_flip);CHKERRQ(ierr);
741         ierr = ISRestoreIndices(Pall,&idxs);CHKERRQ(ierr);
742       }
743       ierr = ISDestroy(&Pall);CHKERRQ(ierr);
744       ierr = ISDestroy(&pII);CHKERRQ(ierr);
745 
746       /* local selected pressures in subdomain-wise and global ordering */
747       ierr = PetscArrayzero(matis->sf_leafdata,n);CHKERRQ(ierr);
748       ierr = PetscArrayzero(matis->sf_rootdata,nl);CHKERRQ(ierr);
749       if (!ploc) {
750         PetscInt *widxs2;
751 
752         if (!pP) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Missing parallel pressure IS");
753         ierr = ISGetLocalSize(pP,&ni);CHKERRQ(ierr);
754         ierr = ISGetIndices(pP,&idxs);CHKERRQ(ierr);
755         for (i=0;i<ni;i++) matis->sf_rootdata[idxs[i]-rst] = 1;
756         ierr = ISRestoreIndices(pP,&idxs);CHKERRQ(ierr);
757         ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr);
758         ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr);
759         for (i=0,ni=0;i<n;i++) if (matis->sf_leafdata[i]) widxs[ni++] = i;
760         ierr = PetscMalloc1(ni,&widxs2);CHKERRQ(ierr);
761         ierr = ISLocalToGlobalMappingApply(l2g,ni,widxs,widxs2);CHKERRQ(ierr);
762         ierr = ISCreateGeneral(PETSC_COMM_SELF,ni,widxs,PETSC_COPY_VALUES,&lP);CHKERRQ(ierr);
763         ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lP",(PetscObject)lP);CHKERRQ(ierr);
764         ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs2,PETSC_OWN_POINTER,&is1);CHKERRQ(ierr);
765         ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_gP",(PetscObject)is1);CHKERRQ(ierr);
766         ierr = ISDestroy(&is1);CHKERRQ(ierr);
767       } else {
768         if (!lP) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing sequential pressure IS");
769         ierr = ISGetLocalSize(lP,&ni);CHKERRQ(ierr);
770         ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr);
771         for (i=0;i<ni;i++)
772           if (idxs[i] >=0 && idxs[i] < n)
773             matis->sf_leafdata[idxs[i]] = 1;
774         ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr);
775         ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr);
776         ierr = ISLocalToGlobalMappingApply(l2g,ni,idxs,widxs);CHKERRQ(ierr);
777         ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr);
778         ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_gP",(PetscObject)is1);CHKERRQ(ierr);
779         ierr = ISDestroy(&is1);CHKERRQ(ierr);
780         ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr);
781         for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst;
782         ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&pP);CHKERRQ(ierr);
783         ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",(PetscObject)pP);CHKERRQ(ierr);
784       }
785       ierr = PetscFree(widxs);CHKERRQ(ierr);
786 
787       /* If there's any "interior pressure",
788          we may want to use a discrete harmonic solver instead
789          of a Stokes harmonic for the Dirichlet preconditioner
790          Need to extract the interior velocity dofs in interior dofs ordering (iV)
791          and interior pressure dofs in local ordering (iP) */
792       if (!allp) {
793         ISLocalToGlobalMapping l2g_t;
794 
795         ierr = ISDifference(lPall,lP,&is1);CHKERRQ(ierr);
796         ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_iP",(PetscObject)is1);CHKERRQ(ierr);
797         ierr = ISDifference(II,is1,&is2);CHKERRQ(ierr);
798         ierr = ISDestroy(&is1);CHKERRQ(ierr);
799         ierr = ISLocalToGlobalMappingCreateIS(II,&l2g_t);CHKERRQ(ierr);
800         ierr = ISGlobalToLocalMappingApplyIS(l2g_t,IS_GTOLM_DROP,is2,&is1);CHKERRQ(ierr);
801         ierr = ISGetLocalSize(is1,&i);CHKERRQ(ierr);
802         ierr = ISGetLocalSize(is2,&j);CHKERRQ(ierr);
803         if (i != j) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent local sizes %D and %D for iV",i,j);
804         ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_iV",(PetscObject)is1);CHKERRQ(ierr);
805         ierr = ISLocalToGlobalMappingDestroy(&l2g_t);CHKERRQ(ierr);
806         ierr = ISDestroy(&is1);CHKERRQ(ierr);
807         ierr = ISDestroy(&is2);CHKERRQ(ierr);
808       }
809       ierr = ISDestroy(&II);CHKERRQ(ierr);
810 
811       /* exclude selected pressures from the inner BDDC */
812       if (pcbddc->DirichletBoundariesLocal) {
813         IS       list[2],plP,isout;
814         PetscInt np;
815 
816         /* need a parallel IS */
817         ierr = ISGetLocalSize(lP,&np);CHKERRQ(ierr);
818         ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr);
819         ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),np,idxs,PETSC_USE_POINTER,&plP);CHKERRQ(ierr);
820         list[0] = plP;
821         list[1] = pcbddc->DirichletBoundariesLocal;
822         ierr = ISConcatenate(PetscObjectComm((PetscObject)ksp),2,list,&isout);CHKERRQ(ierr);
823         ierr = ISSortRemoveDups(isout);CHKERRQ(ierr);
824         ierr = ISDestroy(&plP);CHKERRQ(ierr);
825         ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr);
826         ierr = PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc,isout);CHKERRQ(ierr);
827         ierr = ISDestroy(&isout);CHKERRQ(ierr);
828       } else if (pcbddc->DirichletBoundaries) {
829         IS list[2],isout;
830 
831         list[0] = pP;
832         list[1] = pcbddc->DirichletBoundaries;
833         ierr = ISConcatenate(PetscObjectComm((PetscObject)ksp),2,list,&isout);CHKERRQ(ierr);
834         ierr = ISSortRemoveDups(isout);CHKERRQ(ierr);
835         ierr = PCBDDCSetDirichletBoundaries(fetidp->innerbddc,isout);CHKERRQ(ierr);
836         ierr = ISDestroy(&isout);CHKERRQ(ierr);
837       } else {
838         IS       plP;
839         PetscInt np;
840 
841         /* need a parallel IS */
842         ierr = ISGetLocalSize(lP,&np);CHKERRQ(ierr);
843         ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr);
844         ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),np,idxs,PETSC_COPY_VALUES,&plP);CHKERRQ(ierr);
845         ierr = PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc,plP);CHKERRQ(ierr);
846         ierr = ISDestroy(&plP);CHKERRQ(ierr);
847         ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr);
848       }
849 
850       /* save CSR information for the pressure BDDC solver (if any) */
851       if (schp) {
852         PetscInt np,nt;
853 
854         ierr = MatGetSize(matis->A,&nt,NULL);CHKERRQ(ierr);
855         ierr = ISGetLocalSize(lP,&np);CHKERRQ(ierr);
856         if (np) {
857           PetscInt *xadj = pcbddc->mat_graph->xadj;
858           PetscInt *adjn = pcbddc->mat_graph->adjncy;
859           PetscInt nv = pcbddc->mat_graph->nvtxs_csr;
860 
861           if (nv && nv == nt) {
862             ISLocalToGlobalMapping pmap;
863             PetscInt               *schp_csr,*schp_xadj,*schp_adjn,p;
864             PetscContainer         c;
865 
866             ierr = ISLocalToGlobalMappingCreateIS(lPall,&pmap);CHKERRQ(ierr);
867             ierr = ISGetIndices(lPall,&idxs);CHKERRQ(ierr);
868             for (p = 0, nv = 0; p < np; p++) {
869               PetscInt x,n = idxs[p];
870 
871               ierr = ISGlobalToLocalMappingApply(pmap,IS_GTOLM_DROP,xadj[n+1]-xadj[n],adjn+xadj[n],&x,NULL);CHKERRQ(ierr);
872               nv  += x;
873             }
874             ierr = PetscMalloc1(np + 1 + nv,&schp_csr);CHKERRQ(ierr);
875             schp_xadj = schp_csr;
876             schp_adjn = schp_csr + np + 1;
877             for (p = 0, schp_xadj[0] = 0; p < np; p++) {
878               PetscInt x,n = idxs[p];
879 
880               ierr = ISGlobalToLocalMappingApply(pmap,IS_GTOLM_DROP,xadj[n+1]-xadj[n],adjn+xadj[n],&x,schp_adjn + schp_xadj[p]);CHKERRQ(ierr);
881               schp_xadj[p+1] = schp_xadj[p] + x;
882             }
883             ierr = ISRestoreIndices(lPall,&idxs);CHKERRQ(ierr);
884             ierr = ISLocalToGlobalMappingDestroy(&pmap);CHKERRQ(ierr);
885             ierr = PetscContainerCreate(PETSC_COMM_SELF,&c);CHKERRQ(ierr);
886             ierr = PetscContainerSetPointer(c,schp_csr);CHKERRQ(ierr);
887             ierr = PetscContainerSetUserDestroy(c,PetscContainerUserDestroyDefault);CHKERRQ(ierr);
888             ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pCSR",(PetscObject)c);CHKERRQ(ierr);
889             ierr = PetscContainerDestroy(&c);CHKERRQ(ierr);
890 
891           }
892         }
893       }
894       ierr = ISDestroy(&lPall);CHKERRQ(ierr);
895       ierr = ISDestroy(&lP);CHKERRQ(ierr);
896       fetidp->pP = pP;
897     }
898 
899     /* total number of selected pressure dofs */
900     ierr = ISGetSize(fetidp->pP,&totP);CHKERRQ(ierr);
901 
902     /* Set operator for inner BDDC */
903     if (totP || fetidp->rhs_flip) {
904       ierr = MatDuplicate(A,MAT_COPY_VALUES,&nA);CHKERRQ(ierr);
905     } else {
906       ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
907       nA   = A;
908     }
909     if (fetidp->rhs_flip) {
910       ierr = MatDiagonalScale(nA,fetidp->rhs_flip,NULL);CHKERRQ(ierr);
911       if (totP) {
912         Mat lA2;
913 
914         ierr = MatISGetLocalMat(nA,&lA);CHKERRQ(ierr);
915         ierr = MatDuplicate(lA,MAT_COPY_VALUES,&lA2);CHKERRQ(ierr);
916         ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",(PetscObject)lA2);CHKERRQ(ierr);
917         ierr = MatDestroy(&lA2);CHKERRQ(ierr);
918       }
919     }
920 
921     if (totP) {
922       ierr = MatSetOption(nA,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
923       ierr = MatZeroRowsColumnsIS(nA,fetidp->pP,1.,NULL,NULL);CHKERRQ(ierr);
924     } else {
925       ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",NULL);CHKERRQ(ierr);
926     }
927     ierr = MatGetNearNullSpace(Ap,&nnsp);CHKERRQ(ierr);
928     if (!nnsp) {
929       ierr = MatGetNullSpace(Ap,&nnsp);CHKERRQ(ierr);
930     }
931     if (!nnsp) {
932       ierr = MatGetNearNullSpace(A,&nnsp);CHKERRQ(ierr);
933     }
934     if (!nnsp) {
935       ierr = MatGetNullSpace(A,&nnsp);CHKERRQ(ierr);
936     }
937     ierr = MatSetNearNullSpace(nA,nnsp);CHKERRQ(ierr);
938     ierr = PCSetOperators(fetidp->innerbddc,nA,nA);CHKERRQ(ierr);
939     ierr = MatDestroy(&nA);CHKERRQ(ierr);
940 
941     /* non-zero rhs on interior dofs when applying the preconditioner */
942     if (totP) pcbddc->switch_static = PETSC_TRUE;
943 
944     /* if there are no interface pressures, set inner bddc flag for benign saddle point */
945     if (!totP) {
946       pcbddc->benign_saddle_point = PETSC_TRUE;
947       pcbddc->compute_nonetflux   = PETSC_TRUE;
948     }
949 
950     /* Operators for pressure preconditioner */
951     if (totP) {
952       /* Extract pressure block if needed */
953       if (!pisz) {
954         Mat C;
955         IS  nzrows = NULL;
956 
957         ierr = MatCreateSubMatrix(A,fetidp->pP,fetidp->pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr);
958         ierr = MatFindNonzeroRows(C,&nzrows);CHKERRQ(ierr);
959         if (nzrows) {
960           PetscInt i;
961 
962           ierr = ISGetSize(nzrows,&i);CHKERRQ(ierr);
963           ierr = ISDestroy(&nzrows);CHKERRQ(ierr);
964           if (!i) pisz = PETSC_TRUE;
965         }
966         if (!pisz) {
967           ierr = MatScale(C,-1.);CHKERRQ(ierr); /* i.e. Almost Incompressible Elasticity, Stokes discretized with Q1xQ1_stabilized */
968           ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_C",(PetscObject)C);CHKERRQ(ierr);
969         }
970         ierr = MatDestroy(&C);CHKERRQ(ierr);
971       }
972       /* Divergence mat */
973       if (!pcbddc->divudotp) {
974         Mat       B;
975         IS        P;
976         IS        l2l = NULL;
977         PetscBool save;
978 
979         ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject*)&P);CHKERRQ(ierr);
980         if (!pisz) {
981           IS       F,V;
982           PetscInt m,M;
983 
984           ierr = MatGetOwnershipRange(A,&m,&M);CHKERRQ(ierr);
985           ierr = ISCreateStride(PetscObjectComm((PetscObject)A),M-m,m,1,&F);CHKERRQ(ierr);
986           ierr = ISComplement(P,m,M,&V);CHKERRQ(ierr);
987           ierr = MatCreateSubMatrix(A,P,V,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr);
988           {
989             Mat_IS *Bmatis = (Mat_IS*)B->data;
990             ierr = PetscObjectReference((PetscObject)Bmatis->getsub_cis);CHKERRQ(ierr);
991             l2l  = Bmatis->getsub_cis;
992           }
993           ierr = ISDestroy(&V);CHKERRQ(ierr);
994           ierr = ISDestroy(&F);CHKERRQ(ierr);
995         } else {
996           ierr = MatCreateSubMatrix(A,P,NULL,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr);
997         }
998         save = pcbddc->compute_nonetflux; /* SetDivergenceMat activates nonetflux computation */
999         ierr = PCBDDCSetDivergenceMat(fetidp->innerbddc,B,PETSC_FALSE,l2l);CHKERRQ(ierr);
1000         pcbddc->compute_nonetflux = save;
1001         ierr = MatDestroy(&B);CHKERRQ(ierr);
1002         ierr = ISDestroy(&l2l);CHKERRQ(ierr);
1003       }
1004       if (A != Ap) { /* user has provided a different Pmat, this always superseeds the setter (TODO: is it OK?) */
1005         /* use monolithic operator, we restrict later */
1006         ierr = KSPFETIDPSetPressureOperator(ksp,Ap);CHKERRQ(ierr);
1007       }
1008       ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject*)&PPmat);CHKERRQ(ierr);
1009 
1010       /* PPmat not present, use some default choice */
1011       if (!PPmat) {
1012         Mat C;
1013 
1014         ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_C",(PetscObject*)&C);CHKERRQ(ierr);
1015         if (!schp && C) { /* non-zero pressure block, most likely Almost Incompressible Elasticity */
1016           ierr = KSPFETIDPSetPressureOperator(ksp,C);CHKERRQ(ierr);
1017         } else if (!pisz && schp) { /* we need the whole pressure mass matrix to define the interface BDDC */
1018           IS  P;
1019 
1020           ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject*)&P);CHKERRQ(ierr);
1021           ierr = MatCreateSubMatrix(A,P,P,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr);
1022           ierr = MatScale(C,-1.);CHKERRQ(ierr);
1023           ierr = KSPFETIDPSetPressureOperator(ksp,C);CHKERRQ(ierr);
1024           ierr = MatDestroy(&C);CHKERRQ(ierr);
1025         } else { /* identity (need to be scaled properly by the user using e.g. a Richardson method */
1026           PetscInt nl;
1027 
1028           ierr = ISGetLocalSize(fetidp->pP,&nl);CHKERRQ(ierr);
1029           ierr = MatCreate(PetscObjectComm((PetscObject)ksp),&C);CHKERRQ(ierr);
1030           ierr = MatSetSizes(C,nl,nl,totP,totP);CHKERRQ(ierr);
1031           ierr = MatSetType(C,MATAIJ);CHKERRQ(ierr);
1032           ierr = MatMPIAIJSetPreallocation(C,1,NULL,0,NULL);CHKERRQ(ierr);
1033           ierr = MatSeqAIJSetPreallocation(C,1,NULL);CHKERRQ(ierr);
1034           ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1035           ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1036           ierr = MatShift(C,1.);CHKERRQ(ierr);
1037           ierr = KSPFETIDPSetPressureOperator(ksp,C);CHKERRQ(ierr);
1038           ierr = MatDestroy(&C);CHKERRQ(ierr);
1039         }
1040       }
1041 
1042       /* Preconditioned operator for the pressure block */
1043       ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject*)&PPmat);CHKERRQ(ierr);
1044       if (PPmat) {
1045         Mat      C;
1046         IS       Pall;
1047         PetscInt AM,PAM,PAN,pam,pan,am,an,pl,pIl,pAg,pIg;
1048 
1049         ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject*)&Pall);CHKERRQ(ierr);
1050         ierr = MatGetSize(A,&AM,NULL);CHKERRQ(ierr);
1051         ierr = MatGetSize(PPmat,&PAM,&PAN);CHKERRQ(ierr);
1052         ierr = ISGetSize(Pall,&pAg);CHKERRQ(ierr);
1053         ierr = ISGetSize(fetidp->pP,&pIg);CHKERRQ(ierr);
1054         ierr = MatGetLocalSize(PPmat,&pam,&pan);CHKERRQ(ierr);
1055         ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr);
1056         ierr = ISGetLocalSize(Pall,&pIl);CHKERRQ(ierr);
1057         ierr = ISGetLocalSize(fetidp->pP,&pl);CHKERRQ(ierr);
1058         if (PAM != PAN) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Pressure matrix must be square, unsupported %D x %D",PAM,PAN);
1059         if (pam != pan) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Local sizes of pressure matrix must be equal, unsupported %D x %D",pam,pan);
1060         if (pam != am && pam != pl && pam != pIl) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_USER,"Invalid number of local rows %D for pressure matrix! Supported are %D, %D or %D",pam,am,pl,pIl);
1061         if (pan != an && pan != pl && pan != pIl) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_USER,"Invalid number of local columns %D for pressure matrix! Supported are %D, %D or %D",pan,an,pl,pIl);
1062         if (PAM == AM) { /* monolithic ordering, restrict to pressure */
1063           if (schp) {
1064             ierr = MatCreateSubMatrix(PPmat,Pall,Pall,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr);
1065           } else {
1066             ierr = MatCreateSubMatrix(PPmat,fetidp->pP,fetidp->pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr);
1067           }
1068         } else if (pAg == PAM) { /* global ordering for pressure only */
1069           if (!allp && !schp) { /* solving for interface pressure only */
1070             IS restr;
1071 
1072             ierr = ISRenumber(fetidp->pP,NULL,NULL,&restr);CHKERRQ(ierr);
1073             ierr = MatCreateSubMatrix(PPmat,restr,restr,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr);
1074             ierr = ISDestroy(&restr);CHKERRQ(ierr);
1075           } else {
1076             ierr = PetscObjectReference((PetscObject)PPmat);CHKERRQ(ierr);
1077             C    = PPmat;
1078           }
1079         } else if (pIg == PAM) { /* global ordering for selected pressure only */
1080           if (schp) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Need the entire matrix");
1081           ierr = PetscObjectReference((PetscObject)PPmat);CHKERRQ(ierr);
1082           C    = PPmat;
1083         } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Unable to use the pressure matrix");
1084 
1085         ierr = KSPFETIDPSetPressureOperator(ksp,C);CHKERRQ(ierr);
1086         ierr = MatDestroy(&C);CHKERRQ(ierr);
1087       } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Missing Pmat for pressure block");
1088     } else { /* totP == 0 */
1089       ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",NULL);CHKERRQ(ierr);
1090     }
1091   }
1092   PetscFunctionReturn(0);
1093 }
1094 
KSPSetUp_FETIDP(KSP ksp)1095 static PetscErrorCode KSPSetUp_FETIDP(KSP ksp)
1096 {
1097   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
1098   PC_BDDC        *pcbddc = (PC_BDDC*)fetidp->innerbddc->data;
1099   PetscBool      flg;
1100   PetscErrorCode ierr;
1101 
1102   PetscFunctionBegin;
1103   ierr = KSPFETIDPSetUpOperators(ksp);CHKERRQ(ierr);
1104   /* set up BDDC */
1105   ierr = PCSetErrorIfFailure(fetidp->innerbddc,ksp->errorifnotconverged);CHKERRQ(ierr);
1106   ierr = PCSetUp(fetidp->innerbddc);CHKERRQ(ierr);
1107   /* FETI-DP as it is implemented needs an exact coarse solver */
1108   if (pcbddc->coarse_ksp) {
1109     ierr = KSPSetTolerances(pcbddc->coarse_ksp,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,1000);CHKERRQ(ierr);
1110     ierr = KSPSetNormType(pcbddc->coarse_ksp,KSP_NORM_DEFAULT);CHKERRQ(ierr);
1111   }
1112   /* FETI-DP as it is implemented needs exact local Neumann solvers */
1113   ierr = KSPSetTolerances(pcbddc->ksp_R,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,1000);CHKERRQ(ierr);
1114   ierr = KSPSetNormType(pcbddc->ksp_R,KSP_NORM_DEFAULT);CHKERRQ(ierr);
1115 
1116   /* setup FETI-DP operators
1117      If fetidp->statechanged is true, we need to update the operators
1118      needed in the saddle-point case. This should be replaced
1119      by a better logic when the FETI-DP matrix and preconditioner will
1120      have their own classes */
1121   if (pcbddc->new_primal_space || fetidp->statechanged) {
1122     Mat F; /* the FETI-DP matrix */
1123     PC  D; /* the FETI-DP preconditioner */
1124     ierr = KSPReset(fetidp->innerksp);CHKERRQ(ierr);
1125     ierr = PCBDDCCreateFETIDPOperators(fetidp->innerbddc,fetidp->fully_redundant,((PetscObject)ksp)->prefix,&F,&D);CHKERRQ(ierr);
1126     ierr = KSPSetOperators(fetidp->innerksp,F,F);CHKERRQ(ierr);
1127     ierr = KSPSetTolerances(fetidp->innerksp,ksp->rtol,ksp->abstol,ksp->divtol,ksp->max_it);CHKERRQ(ierr);
1128     ierr = KSPSetPC(fetidp->innerksp,D);CHKERRQ(ierr);
1129     ierr = PetscObjectIncrementTabLevel((PetscObject)D,(PetscObject)fetidp->innerksp,0);CHKERRQ(ierr);
1130     ierr = KSPSetFromOptions(fetidp->innerksp);CHKERRQ(ierr);
1131     ierr = MatCreateVecs(F,&(fetidp->innerksp)->vec_rhs,&(fetidp->innerksp)->vec_sol);CHKERRQ(ierr);
1132     ierr = MatDestroy(&F);CHKERRQ(ierr);
1133     ierr = PCDestroy(&D);CHKERRQ(ierr);
1134     if (fetidp->check) {
1135       PetscViewer viewer;
1136 
1137       if (!pcbddc->dbg_viewer) {
1138         viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp));
1139       } else {
1140         viewer = pcbddc->dbg_viewer;
1141       }
1142       ierr = KSPFETIDPCheckOperators(ksp,viewer);CHKERRQ(ierr);
1143     }
1144   }
1145   fetidp->statechanged     = PETSC_FALSE;
1146   pcbddc->new_primal_space = PETSC_FALSE;
1147 
1148   /* propagate settings to the inner solve */
1149   ierr = KSPGetComputeSingularValues(ksp,&flg);CHKERRQ(ierr);
1150   ierr = KSPSetComputeSingularValues(fetidp->innerksp,flg);CHKERRQ(ierr);
1151   if (ksp->res_hist) {
1152     ierr = KSPSetResidualHistory(fetidp->innerksp,ksp->res_hist,ksp->res_hist_max,ksp->res_hist_reset);CHKERRQ(ierr);
1153   }
1154   ierr = KSPSetErrorIfNotConverged(fetidp->innerksp,ksp->errorifnotconverged);CHKERRQ(ierr);
1155   ierr = KSPSetUp(fetidp->innerksp);CHKERRQ(ierr);
1156   PetscFunctionReturn(0);
1157 }
1158 
KSPSolve_FETIDP(KSP ksp)1159 static PetscErrorCode KSPSolve_FETIDP(KSP ksp)
1160 {
1161   PetscErrorCode     ierr;
1162   Mat                F,A;
1163   MatNullSpace       nsp;
1164   Vec                X,B,Xl,Bl;
1165   KSP_FETIDP         *fetidp = (KSP_FETIDP*)ksp->data;
1166   PC_BDDC            *pcbddc = (PC_BDDC*)fetidp->innerbddc->data;
1167   KSPConvergedReason reason;
1168   PC                 pc;
1169   PCFailedReason     pcreason;
1170 
1171   PetscFunctionBegin;
1172   ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr);
1173   if (fetidp->saddlepoint) {
1174     ierr = PetscCitationsRegister(citation2,&cited2);CHKERRQ(ierr);
1175   }
1176   ierr = KSPGetOperators(ksp,&A,NULL);CHKERRQ(ierr);
1177   ierr = KSPGetRhs(ksp,&B);CHKERRQ(ierr);
1178   ierr = KSPGetSolution(ksp,&X);CHKERRQ(ierr);
1179   ierr = KSPGetOperators(fetidp->innerksp,&F,NULL);CHKERRQ(ierr);
1180   ierr = KSPGetRhs(fetidp->innerksp,&Bl);CHKERRQ(ierr);
1181   ierr = KSPGetSolution(fetidp->innerksp,&Xl);CHKERRQ(ierr);
1182   ierr = PCBDDCMatFETIDPGetRHS(F,B,Bl);CHKERRQ(ierr);
1183   if (ksp->transpose_solve) {
1184     ierr = KSPSolveTranspose(fetidp->innerksp,Bl,Xl);CHKERRQ(ierr);
1185   } else {
1186     ierr = KSPSolve(fetidp->innerksp,Bl,Xl);CHKERRQ(ierr);
1187   }
1188   ierr = KSPGetConvergedReason(fetidp->innerksp,&reason);CHKERRQ(ierr);
1189   ierr = KSPGetPC(fetidp->innerksp,&pc);CHKERRQ(ierr);
1190   ierr = PCGetFailedReason(pc,&pcreason);CHKERRQ(ierr);
1191   if ((reason < 0 && reason != KSP_DIVERGED_ITS) || pcreason) {
1192     PetscInt its;
1193     ierr = KSPGetIterationNumber(fetidp->innerksp,&its);CHKERRQ(ierr);
1194     ksp->reason = KSP_DIVERGED_PC_FAILED;
1195     ierr = VecSetInf(Xl);CHKERRQ(ierr);
1196     ierr = PetscInfo3(ksp,"Inner KSP solve failed: %s %s at iteration %D",KSPConvergedReasons[reason],PCFailedReasons[pcreason],its);CHKERRQ(ierr);
1197   }
1198   ierr = PCBDDCMatFETIDPGetSolution(F,Xl,X);CHKERRQ(ierr);
1199   ierr = MatGetNullSpace(A,&nsp);CHKERRQ(ierr);
1200   if (nsp) {
1201     ierr = MatNullSpaceRemove(nsp,X);CHKERRQ(ierr);
1202   }
1203   /* update ksp with stats from inner ksp */
1204   ierr = KSPGetConvergedReason(fetidp->innerksp,&ksp->reason);CHKERRQ(ierr);
1205   ierr = KSPGetIterationNumber(fetidp->innerksp,&ksp->its);CHKERRQ(ierr);
1206   ksp->totalits += ksp->its;
1207   ierr = KSPGetResidualHistory(fetidp->innerksp,NULL,&ksp->res_hist_len);CHKERRQ(ierr);
1208   /* restore defaults for inner BDDC (Pre/PostSolve flags) */
1209   pcbddc->temp_solution_used        = PETSC_FALSE;
1210   pcbddc->rhs_change                = PETSC_FALSE;
1211   pcbddc->exact_dirichlet_trick_app = PETSC_FALSE;
1212   PetscFunctionReturn(0);
1213 }
1214 
KSPReset_FETIDP(KSP ksp)1215 static PetscErrorCode KSPReset_FETIDP(KSP ksp)
1216 {
1217   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
1218   PC_BDDC        *pcbddc;
1219   PetscErrorCode ierr;
1220 
1221   PetscFunctionBegin;
1222   ierr = ISDestroy(&fetidp->pP);CHKERRQ(ierr);
1223   ierr = VecDestroy(&fetidp->rhs_flip);CHKERRQ(ierr);
1224   /* avoid PCReset that does not take into account ref counting */
1225   ierr = PCDestroy(&fetidp->innerbddc);CHKERRQ(ierr);
1226   ierr = PCCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerbddc);CHKERRQ(ierr);
1227   ierr = PCSetType(fetidp->innerbddc,PCBDDC);CHKERRQ(ierr);
1228   pcbddc = (PC_BDDC*)fetidp->innerbddc->data;
1229   pcbddc->symmetric_primal = PETSC_FALSE;
1230   ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerbddc);CHKERRQ(ierr);
1231   ierr = KSPDestroy(&fetidp->innerksp);CHKERRQ(ierr);
1232   fetidp->saddlepoint  = PETSC_FALSE;
1233   fetidp->matstate     = -1;
1234   fetidp->matnnzstate  = -1;
1235   fetidp->statechanged = PETSC_TRUE;
1236   PetscFunctionReturn(0);
1237 }
1238 
KSPDestroy_FETIDP(KSP ksp)1239 static PetscErrorCode KSPDestroy_FETIDP(KSP ksp)
1240 {
1241   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
1242   PetscErrorCode ierr;
1243 
1244   PetscFunctionBegin;
1245   ierr = KSPReset_FETIDP(ksp);CHKERRQ(ierr);
1246   ierr = PCDestroy(&fetidp->innerbddc);CHKERRQ(ierr);
1247   ierr = KSPDestroy(&fetidp->innerksp);CHKERRQ(ierr);
1248   ierr = PetscFree(fetidp->monctx);CHKERRQ(ierr);
1249   ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetInnerBDDC_C",NULL);CHKERRQ(ierr);
1250   ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerBDDC_C",NULL);CHKERRQ(ierr);
1251   ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerKSP_C",NULL);CHKERRQ(ierr);
1252   ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetPressureOperator_C",NULL);CHKERRQ(ierr);
1253   ierr = PetscFree(ksp->data);CHKERRQ(ierr);
1254   PetscFunctionReturn(0);
1255 }
1256 
KSPView_FETIDP(KSP ksp,PetscViewer viewer)1257 static PetscErrorCode KSPView_FETIDP(KSP ksp,PetscViewer viewer)
1258 {
1259   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
1260   PetscErrorCode ierr;
1261   PetscBool      iascii;
1262 
1263   PetscFunctionBegin;
1264   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
1265   if (iascii) {
1266     ierr = PetscViewerASCIIPrintf(viewer,"  fully redundant: %d\n",fetidp->fully_redundant);CHKERRQ(ierr);
1267     ierr = PetscViewerASCIIPrintf(viewer,"  saddle point:    %d\n",fetidp->saddlepoint);CHKERRQ(ierr);
1268     ierr = PetscViewerASCIIPrintf(viewer,"Inner KSP solver details\n");CHKERRQ(ierr);
1269   }
1270   ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
1271   ierr = KSPView(fetidp->innerksp,viewer);CHKERRQ(ierr);
1272   ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
1273   if (iascii) {
1274     ierr = PetscViewerASCIIPrintf(viewer,"Inner BDDC solver details\n");CHKERRQ(ierr);
1275   }
1276   ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
1277   ierr = PCView(fetidp->innerbddc,viewer);CHKERRQ(ierr);
1278   ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
1279   PetscFunctionReturn(0);
1280 }
1281 
KSPSetFromOptions_FETIDP(PetscOptionItems * PetscOptionsObject,KSP ksp)1282 static PetscErrorCode KSPSetFromOptions_FETIDP(PetscOptionItems *PetscOptionsObject,KSP ksp)
1283 {
1284   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
1285   PetscErrorCode ierr;
1286 
1287   PetscFunctionBegin;
1288   /* set options prefixes for the inner objects, since the parent prefix will be valid at this point */
1289   ierr = PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerksp,((PetscObject)ksp)->prefix);CHKERRQ(ierr);
1290   ierr = PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerksp,"fetidp_");CHKERRQ(ierr);
1291   if (!fetidp->userbddc) {
1292     ierr = PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerbddc,((PetscObject)ksp)->prefix);CHKERRQ(ierr);
1293     ierr = PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerbddc,"fetidp_bddc_");CHKERRQ(ierr);
1294   }
1295   ierr = PetscOptionsHead(PetscOptionsObject,"KSP FETIDP options");CHKERRQ(ierr);
1296   ierr = PetscOptionsBool("-ksp_fetidp_fullyredundant","Use fully redundant multipliers","none",fetidp->fully_redundant,&fetidp->fully_redundant,NULL);CHKERRQ(ierr);
1297   ierr = PetscOptionsBool("-ksp_fetidp_saddlepoint","Activates support for saddle-point problems",NULL,fetidp->saddlepoint,&fetidp->saddlepoint,NULL);CHKERRQ(ierr);
1298   ierr = PetscOptionsBool("-ksp_fetidp_check","Activates verbose debugging output FETI-DP operators",NULL,fetidp->check,&fetidp->check,NULL);CHKERRQ(ierr);
1299   ierr = PetscOptionsTail();CHKERRQ(ierr);
1300   ierr = PCSetFromOptions(fetidp->innerbddc);CHKERRQ(ierr);
1301   PetscFunctionReturn(0);
1302 }
1303 
1304 /*MC
1305      KSPFETIDP - The FETI-DP method
1306 
1307    This class implements the FETI-DP method [1].
1308    The matrix for the KSP must be of type MATIS.
1309    The FETI-DP linear system (automatically generated constructing an internal PCBDDC object) is solved using an internal KSP object.
1310 
1311    Options Database Keys:
1312 +   -ksp_fetidp_fullyredundant <false>   : use a fully redundant set of Lagrange multipliers
1313 .   -ksp_fetidp_saddlepoint <false>      : activates support for saddle point problems, see [2]
1314 .   -ksp_fetidp_saddlepoint_flip <false> : usually, an incompressible Stokes problem is written as
1315                                            | A B^T | | v | = | f |
1316                                            | B 0   | | p | = | g |
1317                                            with B representing -\int_\Omega \nabla \cdot u q.
1318                                            If -ksp_fetidp_saddlepoint_flip is true, the code assumes that the user provides it as
1319                                            | A B^T | | v | = | f |
1320                                            |-B 0   | | p | = |-g |
1321 .   -ksp_fetidp_pressure_field <-1>      : activates support for saddle point problems, and identifies the pressure field id.
1322                                            If this information is not provided, the pressure field is detected by using MatFindZeroDiagonals().
1323 -   -ksp_fetidp_pressure_all <false>     : if false, uses the interface pressures, as described in [2]. If true, uses the entire pressure field.
1324 
1325    Level: Advanced
1326 
1327    Notes:
1328     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.,
1329 .vb
1330       -fetidp_ksp_type gmres -fetidp_bddc_pc_bddc_symmetric false
1331 .ve
1332    will use GMRES for the solution of the linear system on the Lagrange multipliers, generated using a non-symmetric PCBDDC.
1333 
1334    For saddle point problems with continuous pressures, the preconditioned operator for the pressure solver can be specified with KSPFETIDPSetPressureOperator().
1335    Alternatively, the pressure operator is extracted from the precondioned matrix (if it is different from the linear solver matrix).
1336    If none of the above, an identity matrix will be created; the user then needs to scale it through a Richardson solver.
1337    Options for the pressure solver can be prefixed with -fetidp_fielsplit_p_, E.g.
1338 .vb
1339       -fetidp_fielsplit_p_ksp_type preonly -fetidp_fielsplit_p_pc_type lu -fetidp_fielsplit_p_pc_factor_mat_solver_type mumps
1340 .ve
1341    In order to use the deluxe version of FETI-DP, you must customize the inner BDDC operator with -fetidp_bddc_pc_bddc_use_deluxe_scaling -fetidp_bddc_pc_bddc_deluxe_singlemat and use
1342    non-redundant multipliers, i.e. -ksp_fetidp_fullyredundant false. Options for the scaling solver are prefixed by -fetidp_bddelta_, E.g.
1343 .vb
1344       -fetidp_bddelta_pc_factor_mat_solver_type mumps -fetidp_bddelta_pc_type lu
1345 .ve
1346 
1347    Some of the basic options such as the maximum number of iterations and tolerances are automatically passed from this KSP to the inner KSP that actually performs the iterations.
1348 
1349    The converged reason and number of iterations computed are passed from the inner KSP to this KSP at the end of the solution.
1350 
1351    Developer Notes:
1352     Even though this method does not directly use any norms, the user is allowed to set the KSPNormType to any value.
1353     This is so users do not have to change KSPNormType options when they switch from other KSP methods to this one.
1354 
1355    References:
1356 .vb
1357 .  [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
1358 .  [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
1359 .ve
1360 
1361 .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC(), KSPFETIDPGetInnerBDDC(), KSPFETIDPGetInnerKSP()
1362 M*/
KSPCreate_FETIDP(KSP ksp)1363 PETSC_EXTERN PetscErrorCode KSPCreate_FETIDP(KSP ksp)
1364 {
1365   PetscErrorCode ierr;
1366   KSP_FETIDP     *fetidp;
1367   KSP_FETIDPMon  *monctx;
1368   PC_BDDC        *pcbddc;
1369   PC             pc;
1370 
1371   PetscFunctionBegin;
1372   ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,3);CHKERRQ(ierr);
1373   ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,2);CHKERRQ(ierr);
1374   ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr);
1375   ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_RIGHT,2);CHKERRQ(ierr);
1376   ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr);
1377   ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);CHKERRQ(ierr);
1378   ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);CHKERRQ(ierr);
1379 
1380   ierr = PetscNewLog(ksp,&fetidp);CHKERRQ(ierr);
1381   fetidp->matstate     = -1;
1382   fetidp->matnnzstate  = -1;
1383   fetidp->statechanged = PETSC_TRUE;
1384 
1385   ksp->data = (void*)fetidp;
1386   ksp->ops->setup                        = KSPSetUp_FETIDP;
1387   ksp->ops->solve                        = KSPSolve_FETIDP;
1388   ksp->ops->destroy                      = KSPDestroy_FETIDP;
1389   ksp->ops->computeeigenvalues           = KSPComputeEigenvalues_FETIDP;
1390   ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_FETIDP;
1391   ksp->ops->view                         = KSPView_FETIDP;
1392   ksp->ops->setfromoptions               = KSPSetFromOptions_FETIDP;
1393   ksp->ops->buildsolution                = KSPBuildSolution_FETIDP;
1394   ksp->ops->buildresidual                = KSPBuildResidualDefault;
1395   ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
1396   ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr);
1397   /* create the inner KSP for the Lagrange multipliers */
1398   ierr = KSPCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerksp);CHKERRQ(ierr);
1399   ierr = KSPGetPC(fetidp->innerksp,&pc);CHKERRQ(ierr);
1400   ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr);
1401   ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerksp);CHKERRQ(ierr);
1402   /* monitor */
1403   ierr = PetscNew(&monctx);CHKERRQ(ierr);
1404   monctx->parentksp = ksp;
1405   fetidp->monctx = monctx;
1406   ierr = KSPMonitorSet(fetidp->innerksp,KSPMonitor_FETIDP,fetidp->monctx,NULL);CHKERRQ(ierr);
1407   /* create the inner BDDC */
1408   ierr = PCCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerbddc);CHKERRQ(ierr);
1409   ierr = PCSetType(fetidp->innerbddc,PCBDDC);CHKERRQ(ierr);
1410   /* make sure we always obtain a consistent FETI-DP matrix
1411      for symmetric problems, the user can always customize it through the command line */
1412   pcbddc = (PC_BDDC*)fetidp->innerbddc->data;
1413   pcbddc->symmetric_primal = PETSC_FALSE;
1414   ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerbddc);CHKERRQ(ierr);
1415   /* composed functions */
1416   ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetInnerBDDC_C",KSPFETIDPSetInnerBDDC_FETIDP);CHKERRQ(ierr);
1417   ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerBDDC_C",KSPFETIDPGetInnerBDDC_FETIDP);CHKERRQ(ierr);
1418   ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerKSP_C",KSPFETIDPGetInnerKSP_FETIDP);CHKERRQ(ierr);
1419   ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetPressureOperator_C",KSPFETIDPSetPressureOperator_FETIDP);CHKERRQ(ierr);
1420   /* need to call KSPSetUp_FETIDP even with KSP_SETUP_NEWMATRIX */
1421   ksp->setupnewmatrix = PETSC_TRUE;
1422   PetscFunctionReturn(0);
1423 }
1424