xref: /petsc/src/ksp/pc/impls/bddc/bddc.c (revision eab890ac70e3f8cf2d8ed246aa467cd84d643ced) !
1 /* TODOLIST
2 
3    Solvers
4    - Add support for cholesky for coarse solver (similar to local solvers)
5    - Propagate ksp prefixes for solvers to mat objects?
6 
7    User interface
8    - ** DM attached to pc?
9 
10    Debugging output
11    - * Better management of verbosity levels of debugging output
12 
13    Extra
14    - *** Is it possible to work with PCBDDCGraph on boundary indices only (less memory consumed)?
15    - BDDC with MG framework?
16 
17    FETIDP
18    - Move FETIDP code to its own classes
19 
20    MATIS related operations contained in BDDC code
21    - Provide general case for subassembling
22 
23 */
24 
25 #include <../src/ksp/pc/impls/bddc/bddc.h> /*I "petscpc.h" I*/  /* includes for fortran wrappers */
26 #include <../src/ksp/pc/impls/bddc/bddcprivate.h>
27 #include <petscblaslapack.h>
28 
29 /* -------------------------------------------------------------------------- */
30 #undef __FUNCT__
31 #define __FUNCT__ "PCSetFromOptions_BDDC"
32 PetscErrorCode PCSetFromOptions_BDDC(PetscOptions *PetscOptionsObject,PC pc)
33 {
34   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
35   PetscErrorCode ierr;
36 
37   PetscFunctionBegin;
38   ierr = PetscOptionsHead(PetscOptionsObject,"BDDC options");CHKERRQ(ierr);
39   /* Verbose debugging */
40   ierr = PetscOptionsInt("-pc_bddc_check_level","Verbose output for PCBDDC (intended for debug)","none",pcbddc->dbg_flag,&pcbddc->dbg_flag,NULL);CHKERRQ(ierr);
41   /* Primal space cumstomization */
42   ierr = PetscOptionsBool("-pc_bddc_use_local_mat_graph","Use or not adjacency graph of local mat for interface analysis","none",pcbddc->use_local_adj,&pcbddc->use_local_adj,NULL);CHKERRQ(ierr);
43   ierr = PetscOptionsBool("-pc_bddc_use_vertices","Use or not corner dofs in coarse space","none",pcbddc->use_vertices,&pcbddc->use_vertices,NULL);CHKERRQ(ierr);
44   ierr = PetscOptionsBool("-pc_bddc_use_edges","Use or not edge constraints in coarse space","none",pcbddc->use_edges,&pcbddc->use_edges,NULL);CHKERRQ(ierr);
45   ierr = PetscOptionsBool("-pc_bddc_use_faces","Use or not face constraints in coarse space","none",pcbddc->use_faces,&pcbddc->use_faces,NULL);CHKERRQ(ierr);
46   ierr = PetscOptionsBool("-pc_bddc_use_true_nnsp","Use near null space attached to the matrix without modifications","none",pcbddc->use_nnsp_true,&pcbddc->use_nnsp_true,NULL);CHKERRQ(ierr);
47   ierr = PetscOptionsBool("-pc_bddc_use_qr_single","Use QR factorization for single constraints on cc (QR is used when multiple constraints are present)","none",pcbddc->use_qr_single,&pcbddc->use_qr_single,NULL);CHKERRQ(ierr);
48   /* Change of basis */
49   ierr = PetscOptionsBool("-pc_bddc_use_change_of_basis","Use or not internal change of basis on local edge nodes","none",pcbddc->use_change_of_basis,&pcbddc->use_change_of_basis,NULL);CHKERRQ(ierr);
50   ierr = PetscOptionsBool("-pc_bddc_use_change_on_faces","Use or not internal change of basis on local face nodes","none",pcbddc->use_change_on_faces,&pcbddc->use_change_on_faces,NULL);CHKERRQ(ierr);
51   if (!pcbddc->use_change_of_basis) {
52     pcbddc->use_change_on_faces = PETSC_FALSE;
53   }
54   /* Switch between M_2 (default) and M_3 preconditioners (as defined by C. Dohrmann in the ref. article) */
55   ierr = PetscOptionsBool("-pc_bddc_switch_static","Switch on static condensation ops around the interface preconditioner","none",pcbddc->switch_static,&pcbddc->switch_static,NULL);CHKERRQ(ierr);
56   ierr = PetscOptionsInt("-pc_bddc_coarse_redistribute","Number of procs where to redistribute coarse problem","none",pcbddc->redistribute_coarse,&pcbddc->redistribute_coarse,NULL);CHKERRQ(ierr);
57   ierr = PetscOptionsInt("-pc_bddc_coarsening_ratio","Set coarsening ratio used in multilevel coarsening","none",pcbddc->coarsening_ratio,&pcbddc->coarsening_ratio,NULL);CHKERRQ(ierr);
58   ierr = PetscOptionsInt("-pc_bddc_levels","Set maximum number of levels for multilevel","none",pcbddc->max_levels,&pcbddc->max_levels,NULL);CHKERRQ(ierr);
59   ierr = PetscOptionsBool("-pc_bddc_use_coarse_estimates","Use estimated eigenvalues for coarse problem","none",pcbddc->use_coarse_estimates,&pcbddc->use_coarse_estimates,NULL);CHKERRQ(ierr);
60   ierr = PetscOptionsBool("-pc_bddc_use_deluxe_scaling","Use deluxe scaling for BDDC","none",pcbddc->use_deluxe_scaling,&pcbddc->use_deluxe_scaling,NULL);CHKERRQ(ierr);
61   ierr = PetscOptionsBool("-pc_bddc_schur_rebuild","Whether or not the interface graph for Schur principal minors has to be rebuilt (i.e. define the interface without any adjacency)","none",pcbddc->sub_schurs_rebuild,&pcbddc->sub_schurs_rebuild,NULL);CHKERRQ(ierr);
62   ierr = PetscOptionsInt("-pc_bddc_schur_layers","Number of dofs' layers for the computation of principal minors (i.e. -1 uses all dofs)","none",pcbddc->sub_schurs_layers,&pcbddc->sub_schurs_layers,NULL);CHKERRQ(ierr);
63   ierr = PetscOptionsBool("-pc_bddc_schur_use_useradj","Whether or not the CSR graph specified by the user should be used for computing successive layers (default is to use adj of local mat)","none",pcbddc->sub_schurs_use_useradj,&pcbddc->sub_schurs_use_useradj,NULL);CHKERRQ(ierr);
64   ierr = PetscOptionsBool("-pc_bddc_deluxe_faster","Faster application of deluxe scaling (requires extra work during setup)","none",pcbddc->faster_deluxe,&pcbddc->faster_deluxe,NULL);CHKERRQ(ierr);
65   ierr = PetscOptionsReal("-pc_bddc_adaptive_threshold","Threshold to be used for adaptive selection of constraints","none",pcbddc->adaptive_threshold,&pcbddc->adaptive_threshold,NULL);CHKERRQ(ierr);
66   ierr = PetscOptionsInt("-pc_bddc_adaptive_nmin","Minimum number of constraints per connected components","none",pcbddc->adaptive_nmin,&pcbddc->adaptive_nmin,NULL);CHKERRQ(ierr);
67   ierr = PetscOptionsInt("-pc_bddc_adaptive_nmax","Maximum number of constraints per connected components","none",pcbddc->adaptive_nmax,&pcbddc->adaptive_nmax,NULL);CHKERRQ(ierr);
68   ierr = PetscOptionsBool("-pc_bddc_symmetric","Symmetric computation of primal basis functions","none",pcbddc->symmetric_primal,&pcbddc->symmetric_primal,NULL);CHKERRQ(ierr);
69   ierr = PetscOptionsInt("-pc_bddc_coarse_adj","Number of processors where to map the coarse adjacency list","none",pcbddc->coarse_adj_red,&pcbddc->coarse_adj_red,NULL);CHKERRQ(ierr);
70   ierr = PetscOptionsTail();CHKERRQ(ierr);
71   PetscFunctionReturn(0);
72 }
73 /* -------------------------------------------------------------------------- */
74 #undef __FUNCT__
75 #define __FUNCT__ "PCBDDCSetChangeOfBasisMat_BDDC"
76 static PetscErrorCode PCBDDCSetChangeOfBasisMat_BDDC(PC pc, Mat change)
77 {
78   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
79   PetscErrorCode ierr;
80 
81   PetscFunctionBegin;
82   ierr = MatDestroy(&pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr);
83   ierr = PetscObjectReference((PetscObject)change);CHKERRQ(ierr);
84   pcbddc->user_ChangeOfBasisMatrix = change;
85   PetscFunctionReturn(0);
86 }
87 #undef __FUNCT__
88 #define __FUNCT__ "PCBDDCSetChangeOfBasisMat"
89 /*@
90  PCBDDCSetChangeOfBasisMat - Set user defined change of basis for dofs
91 
92    Collective on PC
93 
94    Input Parameters:
95 +  pc - the preconditioning context
96 -  change - the change of basis matrix
97 
98    Level: intermediate
99 
100    Notes:
101 
102 .seealso: PCBDDC
103 @*/
104 PetscErrorCode PCBDDCSetChangeOfBasisMat(PC pc, Mat change)
105 {
106   PetscErrorCode ierr;
107 
108   PetscFunctionBegin;
109   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
110   PetscValidHeaderSpecific(change,MAT_CLASSID,2);
111   PetscCheckSameComm(pc,1,change,2);
112   if (pc->mat) {
113     PetscInt rows_c,cols_c,rows,cols;
114     ierr = MatGetSize(pc->mat,&rows,&cols);CHKERRQ(ierr);
115     ierr = MatGetSize(change,&rows_c,&cols_c);CHKERRQ(ierr);
116     if (rows_c != rows) {
117       SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid number of rows for change of basis matrix! %d != %d",rows_c,rows);
118     }
119     if (cols_c != cols) {
120       SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid number of columns for change of basis matrix! %d != %d",cols_c,cols);
121     }
122     ierr = MatGetLocalSize(pc->mat,&rows,&cols);CHKERRQ(ierr);
123     ierr = MatGetLocalSize(change,&rows_c,&cols_c);CHKERRQ(ierr);
124     if (rows_c != rows) {
125       SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid number of local rows for change of basis matrix! %d != %d",rows_c,rows);
126     }
127     if (cols_c != cols) {
128       SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid number of local columns for change of basis matrix! %d != %d",cols_c,cols);
129     }
130   }
131   ierr = PetscTryMethod(pc,"PCBDDCSetChangeOfBasisMat_C",(PC,Mat),(pc,change));CHKERRQ(ierr);
132   PetscFunctionReturn(0);
133 }
134 /* -------------------------------------------------------------------------- */
135 #undef __FUNCT__
136 #define __FUNCT__ "PCBDDCSetPrimalVerticesLocalIS_BDDC"
137 static PetscErrorCode PCBDDCSetPrimalVerticesLocalIS_BDDC(PC pc, IS PrimalVertices)
138 {
139   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
140   PetscErrorCode ierr;
141 
142   PetscFunctionBegin;
143   ierr = ISDestroy(&pcbddc->user_primal_vertices);CHKERRQ(ierr);
144   ierr = PetscObjectReference((PetscObject)PrimalVertices);CHKERRQ(ierr);
145   pcbddc->user_primal_vertices = PrimalVertices;
146   PetscFunctionReturn(0);
147 }
148 #undef __FUNCT__
149 #define __FUNCT__ "PCBDDCSetPrimalVerticesLocalIS"
150 /*@
151  PCBDDCSetPrimalVerticesLocalIS - Set additional user defined primal vertices in PCBDDC
152 
153    Collective
154 
155    Input Parameters:
156 +  pc - the preconditioning context
157 -  PrimalVertices - index set of primal vertices in local numbering (can be empty)
158 
159    Level: intermediate
160 
161    Notes:
162 
163 .seealso: PCBDDC
164 @*/
165 PetscErrorCode PCBDDCSetPrimalVerticesLocalIS(PC pc, IS PrimalVertices)
166 {
167   PetscErrorCode ierr;
168 
169   PetscFunctionBegin;
170   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
171   PetscValidHeaderSpecific(PrimalVertices,IS_CLASSID,2);
172   PetscCheckSameComm(pc,1,PrimalVertices,2);
173   ierr = PetscTryMethod(pc,"PCBDDCSetPrimalVerticesLocalIS_C",(PC,IS),(pc,PrimalVertices));CHKERRQ(ierr);
174   PetscFunctionReturn(0);
175 }
176 /* -------------------------------------------------------------------------- */
177 #undef __FUNCT__
178 #define __FUNCT__ "PCBDDCSetCoarseningRatio_BDDC"
179 static PetscErrorCode PCBDDCSetCoarseningRatio_BDDC(PC pc,PetscInt k)
180 {
181   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
182 
183   PetscFunctionBegin;
184   pcbddc->coarsening_ratio = k;
185   PetscFunctionReturn(0);
186 }
187 
188 #undef __FUNCT__
189 #define __FUNCT__ "PCBDDCSetCoarseningRatio"
190 /*@
191  PCBDDCSetCoarseningRatio - Set coarsening ratio used in multilevel
192 
193    Logically collective on PC
194 
195    Input Parameters:
196 +  pc - the preconditioning context
197 -  k - coarsening ratio (H/h at the coarser level)
198 
199    Options Database Keys:
200 .    -pc_bddc_coarsening_ratio
201 
202    Level: intermediate
203 
204    Notes:
205      Approximatively k subdomains at the finer level will be aggregated into a single subdomain at the coarser level
206 
207 .seealso: PCBDDC, PCBDDCSetLevels()
208 @*/
209 PetscErrorCode PCBDDCSetCoarseningRatio(PC pc,PetscInt k)
210 {
211   PetscErrorCode ierr;
212 
213   PetscFunctionBegin;
214   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
215   PetscValidLogicalCollectiveInt(pc,k,2);
216   ierr = PetscTryMethod(pc,"PCBDDCSetCoarseningRatio_C",(PC,PetscInt),(pc,k));CHKERRQ(ierr);
217   PetscFunctionReturn(0);
218 }
219 
220 /* The following functions (PCBDDCSetUseExactDirichlet PCBDDCSetLevel) are not public */
221 #undef __FUNCT__
222 #define __FUNCT__ "PCBDDCSetUseExactDirichlet_BDDC"
223 static PetscErrorCode PCBDDCSetUseExactDirichlet_BDDC(PC pc,PetscBool flg)
224 {
225   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
226 
227   PetscFunctionBegin;
228   pcbddc->use_exact_dirichlet_trick = flg;
229   PetscFunctionReturn(0);
230 }
231 
232 #undef __FUNCT__
233 #define __FUNCT__ "PCBDDCSetUseExactDirichlet"
234 PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc,PetscBool flg)
235 {
236   PetscErrorCode ierr;
237 
238   PetscFunctionBegin;
239   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
240   PetscValidLogicalCollectiveBool(pc,flg,2);
241   ierr = PetscTryMethod(pc,"PCBDDCSetUseExactDirichlet_C",(PC,PetscBool),(pc,flg));CHKERRQ(ierr);
242   PetscFunctionReturn(0);
243 }
244 
245 #undef __FUNCT__
246 #define __FUNCT__ "PCBDDCSetLevel_BDDC"
247 static PetscErrorCode PCBDDCSetLevel_BDDC(PC pc,PetscInt level)
248 {
249   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
250 
251   PetscFunctionBegin;
252   pcbddc->current_level = level;
253   PetscFunctionReturn(0);
254 }
255 
256 #undef __FUNCT__
257 #define __FUNCT__ "PCBDDCSetLevel"
258 PetscErrorCode PCBDDCSetLevel(PC pc,PetscInt level)
259 {
260   PetscErrorCode ierr;
261 
262   PetscFunctionBegin;
263   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
264   PetscValidLogicalCollectiveInt(pc,level,2);
265   ierr = PetscTryMethod(pc,"PCBDDCSetLevel_C",(PC,PetscInt),(pc,level));CHKERRQ(ierr);
266   PetscFunctionReturn(0);
267 }
268 
269 #undef __FUNCT__
270 #define __FUNCT__ "PCBDDCSetLevels_BDDC"
271 static PetscErrorCode PCBDDCSetLevels_BDDC(PC pc,PetscInt levels)
272 {
273   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
274 
275   PetscFunctionBegin;
276   pcbddc->max_levels = levels;
277   PetscFunctionReturn(0);
278 }
279 
280 #undef __FUNCT__
281 #define __FUNCT__ "PCBDDCSetLevels"
282 /*@
283  PCBDDCSetLevels - Sets the maximum number of levels for multilevel
284 
285    Logically collective on PC
286 
287    Input Parameters:
288 +  pc - the preconditioning context
289 -  levels - the maximum number of levels (max 9)
290 
291    Options Database Keys:
292 .    -pc_bddc_levels
293 
294    Level: intermediate
295 
296    Notes:
297      Default value is 0, i.e. traditional one-level BDDC
298 
299 .seealso: PCBDDC, PCBDDCSetCoarseningRatio()
300 @*/
301 PetscErrorCode PCBDDCSetLevels(PC pc,PetscInt levels)
302 {
303   PetscErrorCode ierr;
304 
305   PetscFunctionBegin;
306   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
307   PetscValidLogicalCollectiveInt(pc,levels,2);
308   if (levels > 99) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Maximum number of levels for bddc is 99\n");
309   ierr = PetscTryMethod(pc,"PCBDDCSetLevels_C",(PC,PetscInt),(pc,levels));CHKERRQ(ierr);
310   PetscFunctionReturn(0);
311 }
312 /* -------------------------------------------------------------------------- */
313 
314 #undef __FUNCT__
315 #define __FUNCT__ "PCBDDCSetNullSpace_BDDC"
316 static PetscErrorCode PCBDDCSetNullSpace_BDDC(PC pc,MatNullSpace NullSpace)
317 {
318   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
319   PetscErrorCode ierr;
320 
321   PetscFunctionBegin;
322   ierr = PetscObjectReference((PetscObject)NullSpace);CHKERRQ(ierr);
323   ierr = MatNullSpaceDestroy(&pcbddc->NullSpace);CHKERRQ(ierr);
324   pcbddc->NullSpace = NullSpace;
325   PetscFunctionReturn(0);
326 }
327 
328 #undef __FUNCT__
329 #define __FUNCT__ "PCBDDCSetNullSpace"
330 /*@
331  PCBDDCSetNullSpace - Set nullspace for BDDC operator
332 
333    Logically collective on PC and MatNullSpace
334 
335    Input Parameters:
336 +  pc - the preconditioning context
337 -  NullSpace - Null space of the linear operator to be preconditioned (Pmat)
338 
339    Level: intermediate
340 
341    Notes:
342 
343 .seealso: PCBDDC
344 @*/
345 PetscErrorCode PCBDDCSetNullSpace(PC pc,MatNullSpace NullSpace)
346 {
347   PetscErrorCode ierr;
348 
349   PetscFunctionBegin;
350   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
351   PetscValidHeaderSpecific(NullSpace,MAT_NULLSPACE_CLASSID,2);
352   PetscCheckSameComm(pc,1,NullSpace,2);
353   ierr = PetscTryMethod(pc,"PCBDDCSetNullSpace_C",(PC,MatNullSpace),(pc,NullSpace));CHKERRQ(ierr);
354   PetscFunctionReturn(0);
355 }
356 /* -------------------------------------------------------------------------- */
357 
358 #undef __FUNCT__
359 #define __FUNCT__ "PCBDDCSetDirichletBoundaries_BDDC"
360 static PetscErrorCode PCBDDCSetDirichletBoundaries_BDDC(PC pc,IS DirichletBoundaries)
361 {
362   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
363   PetscErrorCode ierr;
364 
365   PetscFunctionBegin;
366   /* last user setting takes precendence -> destroy any other customization */
367   ierr = ISDestroy(&pcbddc->DirichletBoundariesLocal);CHKERRQ(ierr);
368   ierr = ISDestroy(&pcbddc->DirichletBoundaries);CHKERRQ(ierr);
369   ierr = PetscObjectReference((PetscObject)DirichletBoundaries);CHKERRQ(ierr);
370   pcbddc->DirichletBoundaries = DirichletBoundaries;
371   pcbddc->recompute_topography = PETSC_TRUE;
372   PetscFunctionReturn(0);
373 }
374 
375 #undef __FUNCT__
376 #define __FUNCT__ "PCBDDCSetDirichletBoundaries"
377 /*@
378  PCBDDCSetDirichletBoundaries - Set IS defining Dirichlet boundaries for the global problem.
379 
380    Collective
381 
382    Input Parameters:
383 +  pc - the preconditioning context
384 -  DirichletBoundaries - parallel IS defining the Dirichlet boundaries
385 
386    Level: intermediate
387 
388    Notes:
389      Provide the information if you used MatZeroRows/Columns routines. Any process can list any global node
390 
391 .seealso: PCBDDC, PCBDDCSetDirichletBoundariesLocal(), MatZeroRows(), MatZeroRowsColumns()
392 @*/
393 PetscErrorCode PCBDDCSetDirichletBoundaries(PC pc,IS DirichletBoundaries)
394 {
395   PetscErrorCode ierr;
396 
397   PetscFunctionBegin;
398   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
399   PetscValidHeaderSpecific(DirichletBoundaries,IS_CLASSID,2);
400   PetscCheckSameComm(pc,1,DirichletBoundaries,2);
401   ierr = PetscTryMethod(pc,"PCBDDCSetDirichletBoundaries_C",(PC,IS),(pc,DirichletBoundaries));CHKERRQ(ierr);
402   PetscFunctionReturn(0);
403 }
404 /* -------------------------------------------------------------------------- */
405 
406 #undef __FUNCT__
407 #define __FUNCT__ "PCBDDCSetDirichletBoundariesLocal_BDDC"
408 static PetscErrorCode PCBDDCSetDirichletBoundariesLocal_BDDC(PC pc,IS DirichletBoundaries)
409 {
410   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
411   PetscErrorCode ierr;
412 
413   PetscFunctionBegin;
414   /* last user setting takes precendence -> destroy any other customization */
415   ierr = ISDestroy(&pcbddc->DirichletBoundariesLocal);CHKERRQ(ierr);
416   ierr = ISDestroy(&pcbddc->DirichletBoundaries);CHKERRQ(ierr);
417   ierr = PetscObjectReference((PetscObject)DirichletBoundaries);CHKERRQ(ierr);
418   pcbddc->DirichletBoundariesLocal = DirichletBoundaries;
419   pcbddc->recompute_topography = PETSC_TRUE;
420   PetscFunctionReturn(0);
421 }
422 
423 #undef __FUNCT__
424 #define __FUNCT__ "PCBDDCSetDirichletBoundariesLocal"
425 /*@
426  PCBDDCSetDirichletBoundariesLocal - Set IS defining Dirichlet boundaries for the global problem in local ordering.
427 
428    Collective
429 
430    Input Parameters:
431 +  pc - the preconditioning context
432 -  DirichletBoundaries - parallel IS defining the Dirichlet boundaries (in local ordering)
433 
434    Level: intermediate
435 
436    Notes:
437 
438 .seealso: PCBDDC, PCBDDCSetDirichletBoundaries(), MatZeroRows(), MatZeroRowsColumns()
439 @*/
440 PetscErrorCode PCBDDCSetDirichletBoundariesLocal(PC pc,IS DirichletBoundaries)
441 {
442   PetscErrorCode ierr;
443 
444   PetscFunctionBegin;
445   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
446   PetscValidHeaderSpecific(DirichletBoundaries,IS_CLASSID,2);
447   PetscCheckSameComm(pc,1,DirichletBoundaries,2);
448   ierr = PetscTryMethod(pc,"PCBDDCSetDirichletBoundariesLocal_C",(PC,IS),(pc,DirichletBoundaries));CHKERRQ(ierr);
449   PetscFunctionReturn(0);
450 }
451 /* -------------------------------------------------------------------------- */
452 
453 #undef __FUNCT__
454 #define __FUNCT__ "PCBDDCSetNeumannBoundaries_BDDC"
455 static PetscErrorCode PCBDDCSetNeumannBoundaries_BDDC(PC pc,IS NeumannBoundaries)
456 {
457   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
458   PetscErrorCode ierr;
459 
460   PetscFunctionBegin;
461   /* last user setting takes precendence -> destroy any other customization */
462   ierr = ISDestroy(&pcbddc->NeumannBoundariesLocal);CHKERRQ(ierr);
463   ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr);
464   ierr = PetscObjectReference((PetscObject)NeumannBoundaries);CHKERRQ(ierr);
465   pcbddc->NeumannBoundaries = NeumannBoundaries;
466   pcbddc->recompute_topography = PETSC_TRUE;
467   PetscFunctionReturn(0);
468 }
469 
470 #undef __FUNCT__
471 #define __FUNCT__ "PCBDDCSetNeumannBoundaries"
472 /*@
473  PCBDDCSetNeumannBoundaries - Set IS defining Neumann boundaries for the global problem.
474 
475    Collective
476 
477    Input Parameters:
478 +  pc - the preconditioning context
479 -  NeumannBoundaries - parallel IS defining the Neumann boundaries
480 
481    Level: intermediate
482 
483    Notes:
484      Any process can list any global node
485 
486 .seealso: PCBDDC, PCBDDCSetNeumannBoundariesLocal()
487 @*/
488 PetscErrorCode PCBDDCSetNeumannBoundaries(PC pc,IS NeumannBoundaries)
489 {
490   PetscErrorCode ierr;
491 
492   PetscFunctionBegin;
493   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
494   PetscValidHeaderSpecific(NeumannBoundaries,IS_CLASSID,2);
495   PetscCheckSameComm(pc,1,NeumannBoundaries,2);
496   ierr = PetscTryMethod(pc,"PCBDDCSetNeumannBoundaries_C",(PC,IS),(pc,NeumannBoundaries));CHKERRQ(ierr);
497   PetscFunctionReturn(0);
498 }
499 /* -------------------------------------------------------------------------- */
500 
501 #undef __FUNCT__
502 #define __FUNCT__ "PCBDDCSetNeumannBoundariesLocal_BDDC"
503 static PetscErrorCode PCBDDCSetNeumannBoundariesLocal_BDDC(PC pc,IS NeumannBoundaries)
504 {
505   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
506   PetscErrorCode ierr;
507 
508   PetscFunctionBegin;
509   /* last user setting takes precendence -> destroy any other customization */
510   ierr = ISDestroy(&pcbddc->NeumannBoundariesLocal);CHKERRQ(ierr);
511   ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr);
512   ierr = PetscObjectReference((PetscObject)NeumannBoundaries);CHKERRQ(ierr);
513   pcbddc->NeumannBoundariesLocal = NeumannBoundaries;
514   pcbddc->recompute_topography = PETSC_TRUE;
515   PetscFunctionReturn(0);
516 }
517 
518 #undef __FUNCT__
519 #define __FUNCT__ "PCBDDCSetNeumannBoundariesLocal"
520 /*@
521  PCBDDCSetNeumannBoundariesLocal - Set IS defining Neumann boundaries for the global problem in local ordering.
522 
523    Collective
524 
525    Input Parameters:
526 +  pc - the preconditioning context
527 -  NeumannBoundaries - parallel IS defining the subdomain part of Neumann boundaries (in local ordering)
528 
529    Level: intermediate
530 
531    Notes:
532 
533 .seealso: PCBDDC, PCBDDCSetNeumannBoundaries()
534 @*/
535 PetscErrorCode PCBDDCSetNeumannBoundariesLocal(PC pc,IS NeumannBoundaries)
536 {
537   PetscErrorCode ierr;
538 
539   PetscFunctionBegin;
540   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
541   PetscValidHeaderSpecific(NeumannBoundaries,IS_CLASSID,2);
542   PetscCheckSameComm(pc,1,NeumannBoundaries,2);
543   ierr = PetscTryMethod(pc,"PCBDDCSetNeumannBoundariesLocal_C",(PC,IS),(pc,NeumannBoundaries));CHKERRQ(ierr);
544   PetscFunctionReturn(0);
545 }
546 /* -------------------------------------------------------------------------- */
547 
548 #undef __FUNCT__
549 #define __FUNCT__ "PCBDDCGetDirichletBoundaries_BDDC"
550 static PetscErrorCode PCBDDCGetDirichletBoundaries_BDDC(PC pc,IS *DirichletBoundaries)
551 {
552   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
553 
554   PetscFunctionBegin;
555   *DirichletBoundaries = pcbddc->DirichletBoundaries;
556   PetscFunctionReturn(0);
557 }
558 
559 #undef __FUNCT__
560 #define __FUNCT__ "PCBDDCGetDirichletBoundaries"
561 /*@
562  PCBDDCGetDirichletBoundaries - Get parallel IS for Dirichlet boundaries
563 
564    Collective
565 
566    Input Parameters:
567 .  pc - the preconditioning context
568 
569    Output Parameters:
570 .  DirichletBoundaries - index set defining the Dirichlet boundaries
571 
572    Level: intermediate
573 
574    Notes:
575      The IS returned (if any) is the same passed in earlier by the user with PCBDDCSetDirichletBoundaries
576 
577 .seealso: PCBDDC
578 @*/
579 PetscErrorCode PCBDDCGetDirichletBoundaries(PC pc,IS *DirichletBoundaries)
580 {
581   PetscErrorCode ierr;
582 
583   PetscFunctionBegin;
584   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
585   ierr = PetscUseMethod(pc,"PCBDDCGetDirichletBoundaries_C",(PC,IS*),(pc,DirichletBoundaries));CHKERRQ(ierr);
586   PetscFunctionReturn(0);
587 }
588 /* -------------------------------------------------------------------------- */
589 
590 #undef __FUNCT__
591 #define __FUNCT__ "PCBDDCGetDirichletBoundariesLocal_BDDC"
592 static PetscErrorCode PCBDDCGetDirichletBoundariesLocal_BDDC(PC pc,IS *DirichletBoundaries)
593 {
594   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
595 
596   PetscFunctionBegin;
597   *DirichletBoundaries = pcbddc->DirichletBoundariesLocal;
598   PetscFunctionReturn(0);
599 }
600 
601 #undef __FUNCT__
602 #define __FUNCT__ "PCBDDCGetDirichletBoundariesLocal"
603 /*@
604  PCBDDCGetDirichletBoundariesLocal - Get parallel IS for Dirichlet boundaries (in local ordering)
605 
606    Collective
607 
608    Input Parameters:
609 .  pc - the preconditioning context
610 
611    Output Parameters:
612 .  DirichletBoundaries - index set defining the subdomain part of Dirichlet boundaries
613 
614    Level: intermediate
615 
616    Notes:
617      The IS returned could be the same passed in earlier by the user (if provided with PCBDDCSetDirichletBoundariesLocal) or a global-to-local map of the global IS (if provided with PCBDDCSetDirichletBoundaries).
618           In the latter case, the IS will be available after PCSetUp.
619 
620 .seealso: PCBDDC
621 @*/
622 PetscErrorCode PCBDDCGetDirichletBoundariesLocal(PC pc,IS *DirichletBoundaries)
623 {
624   PetscErrorCode ierr;
625 
626   PetscFunctionBegin;
627   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
628   ierr = PetscUseMethod(pc,"PCBDDCGetDirichletBoundariesLocal_C",(PC,IS*),(pc,DirichletBoundaries));CHKERRQ(ierr);
629   PetscFunctionReturn(0);
630 }
631 /* -------------------------------------------------------------------------- */
632 
633 #undef __FUNCT__
634 #define __FUNCT__ "PCBDDCGetNeumannBoundaries_BDDC"
635 static PetscErrorCode PCBDDCGetNeumannBoundaries_BDDC(PC pc,IS *NeumannBoundaries)
636 {
637   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
638 
639   PetscFunctionBegin;
640   *NeumannBoundaries = pcbddc->NeumannBoundaries;
641   PetscFunctionReturn(0);
642 }
643 
644 #undef __FUNCT__
645 #define __FUNCT__ "PCBDDCGetNeumannBoundaries"
646 /*@
647  PCBDDCGetNeumannBoundaries - Get parallel IS for Neumann boundaries
648 
649    Collective
650 
651    Input Parameters:
652 .  pc - the preconditioning context
653 
654    Output Parameters:
655 .  NeumannBoundaries - index set defining the Neumann boundaries
656 
657    Level: intermediate
658 
659    Notes:
660      The IS returned (if any) is the same passed in earlier by the user with PCBDDCSetNeumannBoundaries
661 
662 .seealso: PCBDDC
663 @*/
664 PetscErrorCode PCBDDCGetNeumannBoundaries(PC pc,IS *NeumannBoundaries)
665 {
666   PetscErrorCode ierr;
667 
668   PetscFunctionBegin;
669   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
670   ierr = PetscUseMethod(pc,"PCBDDCGetNeumannBoundaries_C",(PC,IS*),(pc,NeumannBoundaries));CHKERRQ(ierr);
671   PetscFunctionReturn(0);
672 }
673 /* -------------------------------------------------------------------------- */
674 
675 #undef __FUNCT__
676 #define __FUNCT__ "PCBDDCGetNeumannBoundariesLocal_BDDC"
677 static PetscErrorCode PCBDDCGetNeumannBoundariesLocal_BDDC(PC pc,IS *NeumannBoundaries)
678 {
679   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
680 
681   PetscFunctionBegin;
682   *NeumannBoundaries = pcbddc->NeumannBoundariesLocal;
683   PetscFunctionReturn(0);
684 }
685 
686 #undef __FUNCT__
687 #define __FUNCT__ "PCBDDCGetNeumannBoundariesLocal"
688 /*@
689  PCBDDCGetNeumannBoundariesLocal - Get parallel IS for Neumann boundaries (in local ordering)
690 
691    Collective
692 
693    Input Parameters:
694 .  pc - the preconditioning context
695 
696    Output Parameters:
697 .  NeumannBoundaries - index set defining the subdomain part of Neumann boundaries
698 
699    Level: intermediate
700 
701    Notes:
702      The IS returned could be the same passed in earlier by the user (if provided with PCBDDCSetNeumannBoundariesLocal) or a global-to-local map of the global IS (if provided with PCBDDCSetNeumannBoundaries).
703           In the latter case, the IS will be available after PCSetUp.
704 
705 .seealso: PCBDDC
706 @*/
707 PetscErrorCode PCBDDCGetNeumannBoundariesLocal(PC pc,IS *NeumannBoundaries)
708 {
709   PetscErrorCode ierr;
710 
711   PetscFunctionBegin;
712   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
713   ierr = PetscUseMethod(pc,"PCBDDCGetNeumannBoundariesLocal_C",(PC,IS*),(pc,NeumannBoundaries));CHKERRQ(ierr);
714   PetscFunctionReturn(0);
715 }
716 /* -------------------------------------------------------------------------- */
717 
718 #undef __FUNCT__
719 #define __FUNCT__ "PCBDDCSetLocalAdjacencyGraph_BDDC"
720 static PetscErrorCode PCBDDCSetLocalAdjacencyGraph_BDDC(PC pc, PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode)
721 {
722   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
723   PCBDDCGraph    mat_graph = pcbddc->mat_graph;
724   PetscErrorCode ierr;
725 
726   PetscFunctionBegin;
727   /* free old CSR */
728   ierr = PCBDDCGraphResetCSR(mat_graph);CHKERRQ(ierr);
729   /* TODO: PCBDDCGraphSetAdjacency */
730   /* get CSR into graph structure */
731   if (copymode == PETSC_COPY_VALUES) {
732     ierr = PetscMalloc1(nvtxs+1,&mat_graph->xadj);CHKERRQ(ierr);
733     ierr = PetscMalloc1(xadj[nvtxs],&mat_graph->adjncy);CHKERRQ(ierr);
734     ierr = PetscMemcpy(mat_graph->xadj,xadj,(nvtxs+1)*sizeof(PetscInt));CHKERRQ(ierr);
735     ierr = PetscMemcpy(mat_graph->adjncy,adjncy,xadj[nvtxs]*sizeof(PetscInt));CHKERRQ(ierr);
736   } else if (copymode == PETSC_OWN_POINTER) {
737     mat_graph->xadj = (PetscInt*)xadj;
738     mat_graph->adjncy = (PetscInt*)adjncy;
739   }
740   mat_graph->nvtxs_csr = nvtxs;
741   PetscFunctionReturn(0);
742 }
743 
744 #undef __FUNCT__
745 #define __FUNCT__ "PCBDDCSetLocalAdjacencyGraph"
746 /*@
747  PCBDDCSetLocalAdjacencyGraph - Set adjacency structure (CSR graph) of the local matrix
748 
749    Not collective
750 
751    Input Parameters:
752 +  pc - the preconditioning context
753 .  nvtxs - number of local vertices of the graph (i.e., the size of the local problem)
754 .  xadj, adjncy - the CSR graph
755 -  copymode - either PETSC_COPY_VALUES or PETSC_OWN_POINTER.
756 
757    Level: intermediate
758 
759    Notes:
760 
761 .seealso: PCBDDC,PetscCopyMode
762 @*/
763 PetscErrorCode PCBDDCSetLocalAdjacencyGraph(PC pc,PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode)
764 {
765   void (*f)(void) = 0;
766   PetscErrorCode ierr;
767 
768   PetscFunctionBegin;
769   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
770   PetscValidIntPointer(xadj,3);
771   PetscValidIntPointer(adjncy,4);
772   if (copymode != PETSC_COPY_VALUES && copymode != PETSC_OWN_POINTER) {
773     SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported copy mode %d in %s\n",copymode,__FUNCT__);
774   }
775   ierr = PetscTryMethod(pc,"PCBDDCSetLocalAdjacencyGraph_C",(PC,PetscInt,const PetscInt[],const PetscInt[],PetscCopyMode),(pc,nvtxs,xadj,adjncy,copymode));CHKERRQ(ierr);
776   /* free arrays if PCBDDC is not the PC type */
777   ierr = PetscObjectQueryFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",&f);CHKERRQ(ierr);
778   if (!f && copymode == PETSC_OWN_POINTER) {
779     ierr = PetscFree(xadj);CHKERRQ(ierr);
780     ierr = PetscFree(adjncy);CHKERRQ(ierr);
781   }
782   PetscFunctionReturn(0);
783 }
784 /* -------------------------------------------------------------------------- */
785 
786 #undef __FUNCT__
787 #define __FUNCT__ "PCBDDCSetDofsSplittingLocal_BDDC"
788 static PetscErrorCode PCBDDCSetDofsSplittingLocal_BDDC(PC pc,PetscInt n_is, IS ISForDofs[])
789 {
790   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
791   PetscInt i;
792   PetscErrorCode ierr;
793 
794   PetscFunctionBegin;
795   /* Destroy ISes if they were already set */
796   for (i=0;i<pcbddc->n_ISForDofsLocal;i++) {
797     ierr = ISDestroy(&pcbddc->ISForDofsLocal[i]);CHKERRQ(ierr);
798   }
799   ierr = PetscFree(pcbddc->ISForDofsLocal);CHKERRQ(ierr);
800   /* last user setting takes precendence -> destroy any other customization */
801   for (i=0;i<pcbddc->n_ISForDofs;i++) {
802     ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr);
803   }
804   ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr);
805   pcbddc->n_ISForDofs = 0;
806   /* allocate space then set */
807   if (n_is) {
808     ierr = PetscMalloc1(n_is,&pcbddc->ISForDofsLocal);CHKERRQ(ierr);
809   }
810   for (i=0;i<n_is;i++) {
811     ierr = PetscObjectReference((PetscObject)ISForDofs[i]);CHKERRQ(ierr);
812     pcbddc->ISForDofsLocal[i]=ISForDofs[i];
813   }
814   pcbddc->n_ISForDofsLocal=n_is;
815   if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE;
816   pcbddc->recompute_topography = PETSC_TRUE;
817   PetscFunctionReturn(0);
818 }
819 
820 #undef __FUNCT__
821 #define __FUNCT__ "PCBDDCSetDofsSplittingLocal"
822 /*@
823  PCBDDCSetDofsSplittingLocal - Set index sets defining fields of the local subdomain matrix
824 
825    Collective
826 
827    Input Parameters:
828 +  pc - the preconditioning context
829 .  n_is - number of index sets defining the fields
830 -  ISForDofs - array of IS describing the fields in local ordering
831 
832    Level: intermediate
833 
834    Notes:
835      n_is should be the same among processes. Not all nodes need to be listed: unlisted nodes will belong to the complement field.
836 
837 .seealso: PCBDDC
838 @*/
839 PetscErrorCode PCBDDCSetDofsSplittingLocal(PC pc,PetscInt n_is, IS ISForDofs[])
840 {
841   PetscInt       i;
842   PetscErrorCode ierr;
843 
844   PetscFunctionBegin;
845   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
846   PetscValidLogicalCollectiveInt(pc,n_is,2);
847   for (i=0;i<n_is;i++) {
848     PetscCheckSameComm(pc,1,ISForDofs[i],3);
849     PetscValidHeaderSpecific(ISForDofs[i],IS_CLASSID,3);
850   }
851   ierr = PetscTryMethod(pc,"PCBDDCSetDofsSplittingLocal_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));CHKERRQ(ierr);
852   PetscFunctionReturn(0);
853 }
854 /* -------------------------------------------------------------------------- */
855 
856 #undef __FUNCT__
857 #define __FUNCT__ "PCBDDCSetDofsSplitting_BDDC"
858 static PetscErrorCode PCBDDCSetDofsSplitting_BDDC(PC pc,PetscInt n_is, IS ISForDofs[])
859 {
860   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;
861   PetscInt i;
862   PetscErrorCode ierr;
863 
864   PetscFunctionBegin;
865   /* Destroy ISes if they were already set */
866   for (i=0;i<pcbddc->n_ISForDofs;i++) {
867     ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr);
868   }
869   ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr);
870   /* last user setting takes precendence -> destroy any other customization */
871   for (i=0;i<pcbddc->n_ISForDofsLocal;i++) {
872     ierr = ISDestroy(&pcbddc->ISForDofsLocal[i]);CHKERRQ(ierr);
873   }
874   ierr = PetscFree(pcbddc->ISForDofsLocal);CHKERRQ(ierr);
875   pcbddc->n_ISForDofsLocal = 0;
876   /* allocate space then set */
877   if (n_is) {
878     ierr = PetscMalloc1(n_is,&pcbddc->ISForDofs);CHKERRQ(ierr);
879   }
880   for (i=0;i<n_is;i++) {
881     ierr = PetscObjectReference((PetscObject)ISForDofs[i]);CHKERRQ(ierr);
882     pcbddc->ISForDofs[i]=ISForDofs[i];
883   }
884   pcbddc->n_ISForDofs=n_is;
885   if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE;
886   pcbddc->recompute_topography = PETSC_TRUE;
887   PetscFunctionReturn(0);
888 }
889 
890 #undef __FUNCT__
891 #define __FUNCT__ "PCBDDCSetDofsSplitting"
892 /*@
893  PCBDDCSetDofsSplitting - Set index sets defining fields of the global matrix
894 
895    Collective
896 
897    Input Parameters:
898 +  pc - the preconditioning context
899 .  n_is - number of index sets defining the fields
900 -  ISForDofs - array of IS describing the fields in global ordering
901 
902    Level: intermediate
903 
904    Notes:
905      Any process can list any global node. Not all nodes need to be listed: unlisted nodes will belong to the complement field.
906 
907 .seealso: PCBDDC
908 @*/
909 PetscErrorCode PCBDDCSetDofsSplitting(PC pc,PetscInt n_is, IS ISForDofs[])
910 {
911   PetscInt       i;
912   PetscErrorCode ierr;
913 
914   PetscFunctionBegin;
915   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
916   PetscValidLogicalCollectiveInt(pc,n_is,2);
917   for (i=0;i<n_is;i++) {
918     PetscCheckSameComm(pc,1,ISForDofs[i],3);
919     PetscValidHeaderSpecific(ISForDofs[i],IS_CLASSID,3);
920   }
921   ierr = PetscTryMethod(pc,"PCBDDCSetDofsSplitting_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));CHKERRQ(ierr);
922   PetscFunctionReturn(0);
923 }
924 
925 /* -------------------------------------------------------------------------- */
926 #undef __FUNCT__
927 #define __FUNCT__ "PCPreSolve_BDDC"
928 /* -------------------------------------------------------------------------- */
929 /*
930    PCPreSolve_BDDC - Changes the right hand side and (if necessary) the initial
931                      guess if a transformation of basis approach has been selected.
932 
933    Input Parameter:
934 +  pc - the preconditioner contex
935 
936    Application Interface Routine: PCPreSolve()
937 
938    Notes:
939      The interface routine PCPreSolve() is not usually called directly by
940    the user, but instead is called by KSPSolve().
941 */
942 static PetscErrorCode PCPreSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x)
943 {
944   PetscErrorCode ierr;
945   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
946   PC_IS          *pcis = (PC_IS*)(pc->data);
947   Vec            used_vec;
948   PetscBool      copy_rhs = PETSC_TRUE;
949 
950   PetscFunctionBegin;
951   /* if we are working with cg, one dirichlet solve can be avoided during Krylov iterations */
952   if (ksp) {
953     PetscBool iscg;
954     ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPCG,&iscg);CHKERRQ(ierr);
955     if (!iscg) {
956       ierr = PCBDDCSetUseExactDirichlet(pc,PETSC_FALSE);CHKERRQ(ierr);
957     }
958   }
959   /* Creates parallel work vectors used in presolve */
960   if (!pcbddc->original_rhs) {
961     ierr = VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);CHKERRQ(ierr);
962   }
963   if (!pcbddc->temp_solution) {
964     ierr = VecDuplicate(pcis->vec1_global,&pcbddc->temp_solution);CHKERRQ(ierr);
965   }
966 
967   if (x) {
968     ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);
969     used_vec = x;
970   } else { /* it can only happen when calling PCBDDCMatFETIDPGetRHS */
971     ierr = PetscObjectReference((PetscObject)pcbddc->temp_solution);CHKERRQ(ierr);
972     used_vec = pcbddc->temp_solution;
973     ierr = VecSet(used_vec,0.0);CHKERRQ(ierr);
974   }
975 
976   /* hack into ksp data structure since PCPreSolve comes earlier than setting to zero the guess in src/ksp/ksp/interface/itfunc.c */
977   if (ksp) {
978     /* store the flag for the initial guess since it will be restored back during PCPostSolve_BDDC */
979     ierr = KSPGetInitialGuessNonzero(ksp,&pcbddc->ksp_guess_nonzero);CHKERRQ(ierr);
980     if (!pcbddc->ksp_guess_nonzero) {
981       ierr = VecSet(used_vec,0.0);CHKERRQ(ierr);
982     }
983   }
984 
985   pcbddc->rhs_change = PETSC_FALSE;
986 
987   /* Take into account zeroed rows -> change rhs and store solution removed */
988   if (rhs) {
989     IS dirIS = NULL;
990 
991     /* DirichletBoundariesLocal may not be consistent among neighbours; gets a dirichlet dofs IS from graph (may be cached) */
992     ierr = PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph,&dirIS);CHKERRQ(ierr);
993     if (dirIS) {
994       Mat_IS            *matis = (Mat_IS*)pc->pmat->data;
995       PetscInt          dirsize,i,*is_indices;
996       PetscScalar       *array_x;
997       const PetscScalar *array_diagonal;
998 
999       ierr = MatGetDiagonal(pc->pmat,pcis->vec1_global);CHKERRQ(ierr);
1000       ierr = VecPointwiseDivide(pcis->vec1_global,rhs,pcis->vec1_global);CHKERRQ(ierr);
1001       ierr = VecScatterBegin(matis->rctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1002       ierr = VecScatterEnd(matis->rctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1003       ierr = VecScatterBegin(matis->rctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1004       ierr = VecScatterEnd(matis->rctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1005       ierr = ISGetLocalSize(dirIS,&dirsize);CHKERRQ(ierr);
1006       ierr = VecGetArray(pcis->vec1_N,&array_x);CHKERRQ(ierr);
1007       ierr = VecGetArrayRead(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr);
1008       ierr = ISGetIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr);
1009       for (i=0; i<dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]];
1010       ierr = ISRestoreIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr);
1011       ierr = VecRestoreArrayRead(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr);
1012       ierr = VecRestoreArray(pcis->vec1_N,&array_x);CHKERRQ(ierr);
1013       ierr = VecScatterBegin(matis->rctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1014       ierr = VecScatterEnd(matis->rctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1015       pcbddc->rhs_change = PETSC_TRUE;
1016       ierr = ISDestroy(&dirIS);CHKERRQ(ierr);
1017     }
1018   }
1019 
1020   /* remove the computed solution or the initial guess from the rhs */
1021   if (pcbddc->rhs_change || (ksp && pcbddc->ksp_guess_nonzero) ) {
1022     /* store the original rhs */
1023     if (copy_rhs) {
1024       ierr = VecCopy(rhs,pcbddc->original_rhs);CHKERRQ(ierr);
1025       copy_rhs = PETSC_FALSE;
1026     }
1027     pcbddc->rhs_change = PETSC_TRUE;
1028     ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr);
1029     ierr = MatMultAdd(pc->pmat,used_vec,rhs,rhs);CHKERRQ(ierr);
1030     ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr);
1031     ierr = VecCopy(used_vec,pcbddc->temp_solution);CHKERRQ(ierr);
1032     if (ksp) {
1033       ierr = KSPSetInitialGuessNonzero(ksp,PETSC_FALSE);CHKERRQ(ierr);
1034     }
1035   }
1036   ierr = VecDestroy(&used_vec);CHKERRQ(ierr);
1037 
1038   /* store partially computed solution and set initial guess */
1039   if (x && pcbddc->use_exact_dirichlet_trick) {
1040     ierr = VecSet(x,0.0);CHKERRQ(ierr);
1041     ierr = VecScatterBegin(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1042     ierr = VecScatterEnd(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1043     ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr);
1044     ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1045     ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1046     if (ksp) {
1047       ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr);
1048     }
1049   }
1050 
1051   if (pcbddc->ChangeOfBasisMatrix) {
1052     PCBDDCChange_ctx change_ctx;
1053 
1054     /* get change ctx */
1055     ierr = MatShellGetContext(pcbddc->new_global_mat,&change_ctx);CHKERRQ(ierr);
1056 
1057     /* set current iteration matrix inside change context (change of basis has been already set into the ctx during PCSetUp) */
1058     ierr = MatDestroy(&change_ctx->original_mat);CHKERRQ(ierr);
1059     ierr = PetscObjectReference((PetscObject)pc->mat);CHKERRQ(ierr);
1060     change_ctx->original_mat = pc->mat;
1061 
1062     /* change iteration matrix */
1063     ierr = MatDestroy(&pc->mat);CHKERRQ(ierr);
1064     ierr = PetscObjectReference((PetscObject)pcbddc->new_global_mat);CHKERRQ(ierr);
1065     pc->mat = pcbddc->new_global_mat;
1066 
1067     /* store the original rhs */
1068     if (copy_rhs) {
1069       ierr = VecCopy(rhs,pcbddc->original_rhs);CHKERRQ(ierr);
1070       copy_rhs = PETSC_FALSE;
1071     }
1072 
1073     /* change rhs */
1074     ierr = MatMultTranspose(change_ctx->global_change,rhs,pcis->vec1_global);CHKERRQ(ierr);
1075     ierr = VecCopy(pcis->vec1_global,rhs);CHKERRQ(ierr);
1076     pcbddc->rhs_change = PETSC_TRUE;
1077   }
1078 
1079   /* remove nullspace if present */
1080   if (ksp && x && pcbddc->NullSpace) {
1081     ierr = MatNullSpaceRemove(pcbddc->NullSpace,x);CHKERRQ(ierr);
1082     /* store the original rhs */
1083     if (copy_rhs) {
1084       ierr = VecCopy(rhs,pcbddc->original_rhs);CHKERRQ(ierr);
1085       copy_rhs = PETSC_FALSE;
1086     }
1087     pcbddc->rhs_change = PETSC_TRUE;
1088     ierr = MatNullSpaceRemove(pcbddc->NullSpace,rhs);CHKERRQ(ierr);
1089   }
1090   PetscFunctionReturn(0);
1091 }
1092 
1093 /* -------------------------------------------------------------------------- */
1094 #undef __FUNCT__
1095 #define __FUNCT__ "PCPostSolve_BDDC"
1096 /* -------------------------------------------------------------------------- */
1097 /*
1098    PCPostSolve_BDDC - Changes the computed solution if a transformation of basis
1099                      approach has been selected. Also, restores rhs to its original state.
1100 
1101    Input Parameter:
1102 +  pc - the preconditioner contex
1103 
1104    Application Interface Routine: PCPostSolve()
1105 
1106    Notes:
1107      The interface routine PCPostSolve() is not usually called directly by
1108      the user, but instead is called by KSPSolve().
1109 */
1110 static PetscErrorCode PCPostSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x)
1111 {
1112   PetscErrorCode ierr;
1113   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
1114 
1115   PetscFunctionBegin;
1116   if (pcbddc->ChangeOfBasisMatrix) {
1117     PCBDDCChange_ctx change_ctx;
1118 
1119     /* get change ctx */
1120     ierr = MatShellGetContext(pcbddc->new_global_mat,&change_ctx);CHKERRQ(ierr);
1121 
1122     /* restore iteration matrix */
1123     ierr = MatDestroy(&pc->mat);CHKERRQ(ierr);
1124     ierr = PetscObjectReference((PetscObject)change_ctx->original_mat);CHKERRQ(ierr);
1125     pc->mat = change_ctx->original_mat;
1126 
1127     /* get solution in original basis */
1128     if (x) {
1129       PC_IS *pcis = (PC_IS*)(pc->data);
1130       ierr = MatMult(change_ctx->global_change,x,pcis->vec1_global);CHKERRQ(ierr);
1131       ierr = VecCopy(pcis->vec1_global,x);CHKERRQ(ierr);
1132     }
1133   }
1134 
1135   /* add solution removed in presolve */
1136   if (x && pcbddc->rhs_change) {
1137     ierr = VecAXPY(x,1.0,pcbddc->temp_solution);CHKERRQ(ierr);
1138   }
1139 
1140   /* restore rhs to its original state */
1141   if (rhs && pcbddc->rhs_change) {
1142     ierr = VecCopy(pcbddc->original_rhs,rhs);CHKERRQ(ierr);
1143   }
1144   pcbddc->rhs_change = PETSC_FALSE;
1145 
1146   /* restore ksp guess state */
1147   if (ksp) {
1148     ierr = KSPSetInitialGuessNonzero(ksp,pcbddc->ksp_guess_nonzero);CHKERRQ(ierr);
1149   }
1150   PetscFunctionReturn(0);
1151 }
1152 /* -------------------------------------------------------------------------- */
1153 #undef __FUNCT__
1154 #define __FUNCT__ "PCSetUp_BDDC"
1155 /* -------------------------------------------------------------------------- */
1156 /*
1157    PCSetUp_BDDC - Prepares for the use of the BDDC preconditioner
1158                   by setting data structures and options.
1159 
1160    Input Parameter:
1161 +  pc - the preconditioner context
1162 
1163    Application Interface Routine: PCSetUp()
1164 
1165    Notes:
1166      The interface routine PCSetUp() is not usually called directly by
1167      the user, but instead is called by PCApply() if necessary.
1168 */
1169 PetscErrorCode PCSetUp_BDDC(PC pc)
1170 {
1171   PetscErrorCode ierr;
1172   PC_BDDC*       pcbddc = (PC_BDDC*)pc->data;
1173   Mat_IS*        matis;
1174   MatNullSpace   nearnullspace;
1175   PetscInt       nrows,ncols;
1176   PetscBool      computetopography,computesolvers,computesubschurs;
1177   PetscBool      computeconstraintsmatrix;
1178   PetscBool      new_nearnullspace_provided,ismatis;
1179 
1180   PetscFunctionBegin;
1181   ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATIS,&ismatis);CHKERRQ(ierr);
1182   if (!ismatis) {
1183     SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PCBDDC preconditioner requires matrix of type MATIS");
1184   }
1185   ierr = MatGetSize(pc->pmat,&nrows,&ncols);CHKERRQ(ierr);
1186   if (nrows != ncols) {
1187     SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PCBDDC preconditioner requires a square preconditioning matrix");
1188   }
1189   matis = (Mat_IS*)pc->pmat->data;
1190   /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other future nonoverlapping preconditioners */
1191   /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup
1192      Also, BDDC directly build the Dirichlet problem */
1193   /* split work */
1194   if (pc->setupcalled) {
1195     if (pc->flag == SAME_NONZERO_PATTERN) {
1196       computetopography = PETSC_FALSE;
1197       computesolvers = PETSC_TRUE;
1198     } else { /* DIFFERENT_NONZERO_PATTERN */
1199       computetopography = PETSC_TRUE;
1200       computesolvers = PETSC_TRUE;
1201     }
1202   } else {
1203     computetopography = PETSC_TRUE;
1204     computesolvers = PETSC_TRUE;
1205   }
1206   if (pcbddc->recompute_topography) {
1207     computetopography = PETSC_TRUE;
1208   }
1209   computeconstraintsmatrix = PETSC_FALSE;
1210   if (pcbddc->adaptive_threshold > 0.0 && !pcbddc->use_deluxe_scaling) {
1211     SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot compute adaptive constraints without deluxe scaling. Rerun with -pc_bddc_use_deluxe_scaling");
1212   }
1213   pcbddc->adaptive_selection = (PetscBool)(pcbddc->adaptive_threshold > 0.0 && pcbddc->use_deluxe_scaling);
1214   if (pcbddc->adaptive_selection) pcbddc->use_faces = PETSC_TRUE;
1215 
1216   computesubschurs = (PetscBool)(pcbddc->adaptive_selection || pcbddc->use_deluxe_scaling);
1217   if (pcbddc->faster_deluxe && pcbddc->adaptive_selection && pcbddc->use_change_of_basis) {
1218     SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot compute faster deluxe if adaptivity and change of basis are both requested. Rerun with -pc_bddc_deluxe_faster false");
1219   }
1220   /* Get stdout for dbg */
1221   if (pcbddc->dbg_flag) {
1222     if (!pcbddc->dbg_viewer) {
1223       pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc));
1224       ierr = PetscViewerASCIISynchronizedAllow(pcbddc->dbg_viewer,PETSC_TRUE);CHKERRQ(ierr);
1225     }
1226     ierr = PetscViewerASCIIAddTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr);
1227   }
1228 
1229   if (pcbddc->user_ChangeOfBasisMatrix) {
1230     /* use_change_of_basis flag is used to automatically compute a change of basis from constraints */
1231     pcbddc->use_change_of_basis = PETSC_FALSE;
1232     ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr);
1233   } else {
1234     ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr);
1235     ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr);
1236     pcbddc->local_mat = matis->A;
1237   }
1238 
1239   /* workaround for reals */
1240 #if !defined(PETSC_USE_COMPLEX)
1241   if (matis->A->symmetric_set) {
1242     ierr = MatSetOption(pcbddc->local_mat,MAT_HERMITIAN,matis->A->symmetric);CHKERRQ(ierr);
1243   }
1244 #endif
1245 
1246   /* Set up all the "iterative substructuring" common block without computing solvers */
1247   {
1248     Mat temp_mat;
1249 
1250     temp_mat = matis->A;
1251     matis->A = pcbddc->local_mat;
1252     ierr = PCISSetUp(pc,PETSC_FALSE);CHKERRQ(ierr);
1253     pcbddc->local_mat = matis->A;
1254     matis->A = temp_mat;
1255   }
1256 
1257   /* Analyze interface and setup sub_schurs data */
1258   if (computetopography) {
1259     ierr = PCBDDCAnalyzeInterface(pc);CHKERRQ(ierr);
1260     computeconstraintsmatrix = PETSC_TRUE;
1261   }
1262 
1263   /* Setup local dirichlet solver ksp_D and sub_schurs solvers */
1264   if (computesolvers) {
1265     PCBDDCSubSchurs sub_schurs=pcbddc->sub_schurs;
1266 
1267     if (computesubschurs && computetopography) {
1268       ierr = PCBDDCInitSubSchurs(pc);CHKERRQ(ierr);
1269     }
1270     if (sub_schurs->use_mumps) {
1271       if (computesubschurs) {
1272         ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr);
1273       }
1274       ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr);
1275     } else {
1276       ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr);
1277       if (computesubschurs) {
1278         ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr);
1279       }
1280     }
1281     if (pcbddc->adaptive_selection) {
1282       ierr = PCBDDCAdaptiveSelection(pc);CHKERRQ(ierr);
1283       computeconstraintsmatrix = PETSC_TRUE;
1284     }
1285   }
1286 
1287   /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */
1288   new_nearnullspace_provided = PETSC_FALSE;
1289   ierr = MatGetNearNullSpace(pc->pmat,&nearnullspace);CHKERRQ(ierr);
1290   if (pcbddc->onearnullspace) { /* already used nearnullspace */
1291     if (!nearnullspace) { /* near null space attached to mat has been destroyed */
1292       new_nearnullspace_provided = PETSC_TRUE;
1293     } else {
1294       /* determine if the two nullspaces are different (should be lightweight) */
1295       if (nearnullspace != pcbddc->onearnullspace) {
1296         new_nearnullspace_provided = PETSC_TRUE;
1297       } else { /* maybe the user has changed the content of the nearnullspace so check vectors ObjectStateId */
1298         PetscInt         i;
1299         const Vec        *nearnullvecs;
1300         PetscObjectState state;
1301         PetscInt         nnsp_size;
1302         ierr = MatNullSpaceGetVecs(nearnullspace,NULL,&nnsp_size,&nearnullvecs);CHKERRQ(ierr);
1303         for (i=0;i<nnsp_size;i++) {
1304           ierr = PetscObjectStateGet((PetscObject)nearnullvecs[i],&state);CHKERRQ(ierr);
1305           if (pcbddc->onearnullvecs_state[i] != state) {
1306             new_nearnullspace_provided = PETSC_TRUE;
1307             break;
1308           }
1309         }
1310       }
1311     }
1312   } else {
1313     if (!nearnullspace) { /* both nearnullspaces are null */
1314       new_nearnullspace_provided = PETSC_FALSE;
1315     } else { /* nearnullspace attached later */
1316       new_nearnullspace_provided = PETSC_TRUE;
1317     }
1318   }
1319 
1320   /* Setup constraints and related work vectors */
1321   /* reset primal space flags */
1322   pcbddc->new_primal_space = PETSC_FALSE;
1323   pcbddc->new_primal_space_local = PETSC_FALSE;
1324   if (computeconstraintsmatrix || new_nearnullspace_provided) {
1325     /* It also sets the primal space flags */
1326     ierr = PCBDDCConstraintsSetUp(pc);CHKERRQ(ierr);
1327     /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */
1328     ierr = PCBDDCSetUpLocalWorkVectors(pc);CHKERRQ(ierr);
1329   }
1330 
1331   if (computesolvers || pcbddc->new_primal_space) {
1332     if (pcbddc->use_change_of_basis) {
1333       PC_IS *pcis = (PC_IS*)(pc->data);
1334 
1335       ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr);
1336       /* get submatrices */
1337       ierr = MatDestroy(&pcis->A_IB);CHKERRQ(ierr);
1338       ierr = MatDestroy(&pcis->A_BI);CHKERRQ(ierr);
1339       ierr = MatDestroy(&pcis->A_BB);CHKERRQ(ierr);
1340       ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_BB);CHKERRQ(ierr);
1341       ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_I_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_IB);CHKERRQ(ierr);
1342       ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_BI);CHKERRQ(ierr);
1343       /* set flag in pcis to not reuse submatrices during PCISCreate */
1344       pcis->reusesubmatrices = PETSC_FALSE;
1345     } else if (!pcbddc->user_ChangeOfBasisMatrix) {
1346       ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr);
1347       ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr);
1348       pcbddc->local_mat = matis->A;
1349     }
1350     /* SetUp coarse and local Neumann solvers */
1351     ierr = PCBDDCSetUpSolvers(pc);CHKERRQ(ierr);
1352     /* SetUp Scaling operator */
1353     ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr);
1354   }
1355 
1356   if (pcbddc->dbg_flag) {
1357     ierr = PetscViewerASCIISubtractTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr);
1358   }
1359   PetscFunctionReturn(0);
1360 }
1361 
1362 /* -------------------------------------------------------------------------- */
1363 /*
1364    PCApply_BDDC - Applies the BDDC operator to a vector.
1365 
1366    Input Parameters:
1367 +  pc - the preconditioner context
1368 -  r - input vector (global)
1369 
1370    Output Parameter:
1371 .  z - output vector (global)
1372 
1373    Application Interface Routine: PCApply()
1374  */
1375 #undef __FUNCT__
1376 #define __FUNCT__ "PCApply_BDDC"
1377 PetscErrorCode PCApply_BDDC(PC pc,Vec r,Vec z)
1378 {
1379   PC_IS             *pcis = (PC_IS*)(pc->data);
1380   PC_BDDC           *pcbddc = (PC_BDDC*)(pc->data);
1381   PetscInt          n_B = pcis->n_B, n_D = pcis->n - n_B;
1382   PetscErrorCode    ierr;
1383   const PetscScalar one = 1.0;
1384   const PetscScalar m_one = -1.0;
1385   const PetscScalar zero = 0.0;
1386 
1387 /* This code is similar to that provided in nn.c for PCNN
1388    NN interface preconditioner changed to BDDC
1389    Added support for M_3 preconditioner in the reference article (code is active if pcbddc->switch_static == PETSC_TRUE) */
1390 
1391   PetscFunctionBegin;
1392   if (!pcbddc->use_exact_dirichlet_trick) {
1393     ierr = VecCopy(r,z);CHKERRQ(ierr);
1394     /* First Dirichlet solve */
1395     ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1396     ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1397     /*
1398       Assembling right hand side for BDDC operator
1399       - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE)
1400       - pcis->vec1_B the interface part of the global vector z
1401     */
1402     if (n_D) {
1403       ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr);
1404       ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr);
1405       if (pcbddc->switch_static) { ierr = MatMultAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); }
1406       ierr = MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr);
1407     } else {
1408       ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr);
1409     }
1410     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1411     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1412     ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr);
1413   } else {
1414     if (pcbddc->switch_static) {
1415       ierr = VecSet(pcis->vec1_D,zero);CHKERRQ(ierr);
1416     }
1417     ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr);
1418   }
1419 
1420   /* Apply interface preconditioner
1421      input/output vecs: pcis->vec1_B and pcis->vec1_D */
1422   ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_FALSE);CHKERRQ(ierr);
1423 
1424   /* Apply transpose of partition of unity operator */
1425   ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr);
1426 
1427   /* Second Dirichlet solve and assembling of output */
1428   ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1429   ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1430   if (n_B) {
1431     ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr);
1432     if (pcbddc->switch_static) { ierr = MatMultAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D);CHKERRQ(ierr); }
1433   } else if (pcbddc->switch_static) {
1434     ierr = MatMult(pcis->A_II,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr);
1435   }
1436   ierr = KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr);
1437   if (!pcbddc->use_exact_dirichlet_trick) {
1438     if (pcbddc->switch_static) {
1439       ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr);
1440     } else {
1441       ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr);
1442     }
1443     ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1444     ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1445   } else {
1446     if (pcbddc->switch_static) {
1447       ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr);
1448     } else {
1449       ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr);
1450     }
1451     ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1452     ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1453   }
1454   PetscFunctionReturn(0);
1455 }
1456 
1457 /* -------------------------------------------------------------------------- */
1458 /*
1459    PCApplyTranspose_BDDC - Applies the transpose of the BDDC operator to a vector.
1460 
1461    Input Parameters:
1462 +  pc - the preconditioner context
1463 -  r - input vector (global)
1464 
1465    Output Parameter:
1466 .  z - output vector (global)
1467 
1468    Application Interface Routine: PCApplyTranspose()
1469  */
1470 #undef __FUNCT__
1471 #define __FUNCT__ "PCApplyTranspose_BDDC"
1472 PetscErrorCode PCApplyTranspose_BDDC(PC pc,Vec r,Vec z)
1473 {
1474   PC_IS             *pcis = (PC_IS*)(pc->data);
1475   PC_BDDC           *pcbddc = (PC_BDDC*)(pc->data);
1476   PetscInt          n_B = pcis->n_B, n_D = pcis->n - n_B;
1477   PetscErrorCode    ierr;
1478   const PetscScalar one = 1.0;
1479   const PetscScalar m_one = -1.0;
1480   const PetscScalar zero = 0.0;
1481 
1482   PetscFunctionBegin;
1483   if (!pcbddc->use_exact_dirichlet_trick) {
1484     ierr = VecCopy(r,z);CHKERRQ(ierr);
1485     /* First Dirichlet solve */
1486     ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1487     ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1488     /*
1489       Assembling right hand side for BDDC operator
1490       - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE)
1491       - pcis->vec1_B the interface part of the global vector z
1492     */
1493     if (n_D) {
1494       ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr);
1495       ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr);
1496       if (pcbddc->switch_static) { ierr = MatMultTransposeAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); }
1497       ierr = MatMultTranspose(pcis->A_IB,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr);
1498     } else {
1499       ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr);
1500     }
1501     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1502     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1503     ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr);
1504   } else {
1505     if (pcbddc->switch_static) {
1506       ierr = VecSet(pcis->vec1_D,zero);CHKERRQ(ierr);
1507     }
1508     ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr);
1509   }
1510 
1511   /* Apply interface preconditioner
1512      input/output vecs: pcis->vec1_B and pcis->vec1_D */
1513   ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_TRUE);CHKERRQ(ierr);
1514 
1515   /* Apply transpose of partition of unity operator */
1516   ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr);
1517 
1518   /* Second Dirichlet solve and assembling of output */
1519   ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1520   ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1521   if (n_B) {
1522     ierr = MatMultTranspose(pcis->A_BI,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr);
1523     if (pcbddc->switch_static) { ierr = MatMultTransposeAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D);CHKERRQ(ierr); }
1524   } else if (pcbddc->switch_static) {
1525     ierr = MatMultTranspose(pcis->A_II,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr);
1526   }
1527   ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr);
1528   if (!pcbddc->use_exact_dirichlet_trick) {
1529     if (pcbddc->switch_static) {
1530       ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr);
1531     } else {
1532       ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr);
1533     }
1534     ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1535     ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1536   } else {
1537     if (pcbddc->switch_static) {
1538       ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr);
1539     } else {
1540       ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr);
1541     }
1542     ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1543     ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1544   }
1545   PetscFunctionReturn(0);
1546 }
1547 /* -------------------------------------------------------------------------- */
1548 
1549 #undef __FUNCT__
1550 #define __FUNCT__ "PCDestroy_BDDC"
1551 PetscErrorCode PCDestroy_BDDC(PC pc)
1552 {
1553   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
1554   PetscErrorCode ierr;
1555 
1556   PetscFunctionBegin;
1557   /* free data created by PCIS */
1558   ierr = PCISDestroy(pc);CHKERRQ(ierr);
1559   /* free BDDC custom data  */
1560   ierr = PCBDDCResetCustomization(pc);CHKERRQ(ierr);
1561   /* destroy objects related to topography */
1562   ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr);
1563   /* free allocated graph structure */
1564   ierr = PetscFree(pcbddc->mat_graph);CHKERRQ(ierr);
1565   /* free allocated sub schurs structure */
1566   ierr = PetscFree(pcbddc->sub_schurs);CHKERRQ(ierr);
1567   /* destroy objects for scaling operator */
1568   ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr);
1569   ierr = PetscFree(pcbddc->deluxe_ctx);CHKERRQ(ierr);
1570   /* free solvers stuff */
1571   ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr);
1572   /* free global vectors needed in presolve */
1573   ierr = VecDestroy(&pcbddc->temp_solution);CHKERRQ(ierr);
1574   ierr = VecDestroy(&pcbddc->original_rhs);CHKERRQ(ierr);
1575   /* free stuff for change of basis hooks */
1576   if (pcbddc->new_global_mat) {
1577     PCBDDCChange_ctx change_ctx;
1578     ierr = MatShellGetContext(pcbddc->new_global_mat,&change_ctx);CHKERRQ(ierr);
1579     ierr = MatDestroy(&change_ctx->original_mat);CHKERRQ(ierr);
1580     ierr = MatDestroy(&change_ctx->global_change);CHKERRQ(ierr);
1581     ierr = VecDestroyVecs(2,&change_ctx->work);CHKERRQ(ierr);
1582     ierr = PetscFree(change_ctx);CHKERRQ(ierr);
1583   }
1584   ierr = MatDestroy(&pcbddc->new_global_mat);CHKERRQ(ierr);
1585   /* remove functions */
1586   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",NULL);CHKERRQ(ierr);
1587   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);CHKERRQ(ierr);
1588   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);CHKERRQ(ierr);
1589   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",NULL);CHKERRQ(ierr);
1590   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",NULL);CHKERRQ(ierr);
1591   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",NULL);CHKERRQ(ierr);
1592   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",NULL);CHKERRQ(ierr);
1593   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",NULL);CHKERRQ(ierr);
1594   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr);
1595   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",NULL);CHKERRQ(ierr);
1596   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr);
1597   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",NULL);CHKERRQ(ierr);
1598   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr);
1599   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",NULL);CHKERRQ(ierr);
1600   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr);
1601   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);CHKERRQ(ierr);
1602   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",NULL);CHKERRQ(ierr);
1603   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);CHKERRQ(ierr);
1604   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);CHKERRQ(ierr);
1605   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);CHKERRQ(ierr);
1606   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);CHKERRQ(ierr);
1607   /* Free the private data structure */
1608   ierr = PetscFree(pc->data);CHKERRQ(ierr);
1609   PetscFunctionReturn(0);
1610 }
1611 /* -------------------------------------------------------------------------- */
1612 
1613 #undef __FUNCT__
1614 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS_BDDC"
1615 static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
1616 {
1617   FETIDPMat_ctx  mat_ctx;
1618   Vec            copy_standard_rhs;
1619   PC_IS*         pcis;
1620   PC_BDDC*       pcbddc;
1621   PetscErrorCode ierr;
1622 
1623   PetscFunctionBegin;
1624   ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr);
1625   pcis = (PC_IS*)mat_ctx->pc->data;
1626   pcbddc = (PC_BDDC*)mat_ctx->pc->data;
1627 
1628   /*
1629      change of basis for physical rhs if needed
1630      It also changes the rhs in case of dirichlet boundaries
1631      TODO: better management when FETIDP will have its own class
1632   */
1633   ierr = VecDuplicate(standard_rhs,&copy_standard_rhs);CHKERRQ(ierr);
1634   ierr = VecCopy(standard_rhs,copy_standard_rhs);CHKERRQ(ierr);
1635   ierr = PCPreSolve_BDDC(mat_ctx->pc,NULL,copy_standard_rhs,NULL);CHKERRQ(ierr);
1636   /* store vectors for computation of fetidp final solution */
1637   ierr = VecScatterBegin(pcis->global_to_D,copy_standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1638   ierr = VecScatterEnd(pcis->global_to_D,copy_standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1639   /* scale rhs since it should be unassembled */
1640   /* TODO use counter scaling? (also below) */
1641   ierr = VecScatterBegin(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1642   ierr = VecScatterEnd(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1643   /* Apply partition of unity */
1644   ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr);
1645   /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,copy_standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */
1646   if (!pcbddc->switch_static) {
1647     /* compute partially subassembled Schur complement right-hand side */
1648     ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr);
1649     ierr = MatMult(pcis->A_BI,pcis->vec1_D,pcis->vec1_B);CHKERRQ(ierr);
1650     ierr = VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);CHKERRQ(ierr);
1651     ierr = VecSet(copy_standard_rhs,0.0);CHKERRQ(ierr);
1652     ierr = VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,copy_standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1653     ierr = VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,copy_standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1654     /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,copy_standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */
1655     ierr = VecScatterBegin(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1656     ierr = VecScatterEnd(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1657     ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr);
1658   }
1659   ierr = VecDestroy(&copy_standard_rhs);CHKERRQ(ierr);
1660   /* BDDC rhs */
1661   ierr = VecCopy(mat_ctx->temp_solution_B,pcis->vec1_B);CHKERRQ(ierr);
1662   if (pcbddc->switch_static) {
1663     ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr);
1664   }
1665   /* apply BDDC */
1666   ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr);
1667   /* Application of B_delta and assembling of rhs for fetidp fluxes */
1668   ierr = VecSet(fetidp_flux_rhs,0.0);CHKERRQ(ierr);
1669   ierr = MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);CHKERRQ(ierr);
1670   ierr = VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1671   ierr = VecScatterEnd(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1672   PetscFunctionReturn(0);
1673 }
1674 
1675 #undef __FUNCT__
1676 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS"
1677 /*@
1678  PCBDDCMatFETIDPGetRHS - Compute the right-hand side for FETI-DP linear system using the physical right-hand side
1679 
1680    Collective
1681 
1682    Input Parameters:
1683 +  fetidp_mat      - the FETI-DP matrix object obtained by a call to PCBDDCCreateFETIDPOperators
1684 -  standard_rhs    - the right-hand side of the original linear system
1685 
1686    Output Parameters:
1687 .  fetidp_flux_rhs - the right-hand side for the FETI-DP linear system
1688 
1689    Level: developer
1690 
1691    Notes:
1692 
1693 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetSolution
1694 @*/
1695 PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
1696 {
1697   FETIDPMat_ctx  mat_ctx;
1698   PetscErrorCode ierr;
1699 
1700   PetscFunctionBegin;
1701   ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr);
1702   ierr = PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetRHS_C",(Mat,Vec,Vec),(fetidp_mat,standard_rhs,fetidp_flux_rhs));CHKERRQ(ierr);
1703   PetscFunctionReturn(0);
1704 }
1705 /* -------------------------------------------------------------------------- */
1706 
1707 #undef __FUNCT__
1708 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution_BDDC"
1709 static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
1710 {
1711   FETIDPMat_ctx  mat_ctx;
1712   PC_IS*         pcis;
1713   PC_BDDC*       pcbddc;
1714   PetscErrorCode ierr;
1715 
1716   PetscFunctionBegin;
1717   ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr);
1718   pcis = (PC_IS*)mat_ctx->pc->data;
1719   pcbddc = (PC_BDDC*)mat_ctx->pc->data;
1720 
1721   /* apply B_delta^T */
1722   ierr = VecScatterBegin(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1723   ierr = VecScatterEnd  (mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1724   ierr = MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr);
1725   /* compute rhs for BDDC application */
1726   ierr = VecAYPX(pcis->vec1_B,-1.0,mat_ctx->temp_solution_B);CHKERRQ(ierr);
1727   if (pcbddc->switch_static) {
1728     ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr);
1729   }
1730   /* apply BDDC */
1731   ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr);
1732   /* put values into standard global vector */
1733   ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1734   ierr = VecScatterEnd  (pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1735   if (!pcbddc->switch_static) {
1736     /* compute values into the interior if solved for the partially subassembled Schur complement */
1737     ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);CHKERRQ(ierr);
1738     ierr = VecAXPY(mat_ctx->temp_solution_D,-1.0,pcis->vec1_D);CHKERRQ(ierr);
1739     ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr);
1740   }
1741   ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1742   ierr = VecScatterEnd  (pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1743   /* final change of basis if needed
1744      Is also sums the dirichlet part removed during RHS assembling */
1745   ierr = PCPostSolve_BDDC(mat_ctx->pc,NULL,NULL,standard_sol);CHKERRQ(ierr);
1746   PetscFunctionReturn(0);
1747 }
1748 
1749 #undef __FUNCT__
1750 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution"
1751 /*@
1752  PCBDDCMatFETIDPGetSolution - Compute the physical solution using the solution of the FETI-DP linear system
1753 
1754    Collective
1755 
1756    Input Parameters:
1757 +  fetidp_mat      - the FETI-DP matrix obtained by a call to PCBDDCCreateFETIDPOperators
1758 -  fetidp_flux_sol - the solution of the FETI-DP linear system
1759 
1760    Output Parameters:
1761 .  standard_sol    - the solution defined on the physical domain
1762 
1763    Level: developer
1764 
1765    Notes:
1766 
1767 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetRHS
1768 @*/
1769 PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
1770 {
1771   FETIDPMat_ctx  mat_ctx;
1772   PetscErrorCode ierr;
1773 
1774   PetscFunctionBegin;
1775   ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr);
1776   ierr = PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetSolution_C",(Mat,Vec,Vec),(fetidp_mat,fetidp_flux_sol,standard_sol));CHKERRQ(ierr);
1777   PetscFunctionReturn(0);
1778 }
1779 /* -------------------------------------------------------------------------- */
1780 
1781 extern PetscErrorCode FETIDPMatMult(Mat,Vec,Vec);
1782 extern PetscErrorCode FETIDPMatMultTranspose(Mat,Vec,Vec);
1783 extern PetscErrorCode PCBDDCDestroyFETIDPMat(Mat);
1784 extern PetscErrorCode FETIDPPCApply(PC,Vec,Vec);
1785 extern PetscErrorCode FETIDPPCApplyTranspose(PC,Vec,Vec);
1786 extern PetscErrorCode PCBDDCDestroyFETIDPPC(PC);
1787 
1788 #undef __FUNCT__
1789 #define __FUNCT__ "PCBDDCCreateFETIDPOperators_BDDC"
1790 static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, Mat *fetidp_mat, PC *fetidp_pc)
1791 {
1792 
1793   FETIDPMat_ctx  fetidpmat_ctx;
1794   Mat            newmat;
1795   FETIDPPC_ctx   fetidppc_ctx;
1796   PC             newpc;
1797   MPI_Comm       comm;
1798   PetscErrorCode ierr;
1799 
1800   PetscFunctionBegin;
1801   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
1802   /* FETIDP linear matrix */
1803   ierr = PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);CHKERRQ(ierr);
1804   ierr = PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);CHKERRQ(ierr);
1805   ierr = MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,fetidpmat_ctx->n_lambda,fetidpmat_ctx->n_lambda,fetidpmat_ctx,&newmat);CHKERRQ(ierr);
1806   ierr = MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);CHKERRQ(ierr);
1807   ierr = MatShellSetOperation(newmat,MATOP_MULT_TRANSPOSE,(void (*)(void))FETIDPMatMultTranspose);CHKERRQ(ierr);
1808   ierr = MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);CHKERRQ(ierr);
1809   ierr = MatSetUp(newmat);CHKERRQ(ierr);
1810   /* FETIDP preconditioner */
1811   ierr = PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);CHKERRQ(ierr);
1812   ierr = PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);CHKERRQ(ierr);
1813   ierr = PCCreate(comm,&newpc);CHKERRQ(ierr);
1814   ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr);
1815   ierr = PCShellSetContext(newpc,fetidppc_ctx);CHKERRQ(ierr);
1816   ierr = PCShellSetApply(newpc,FETIDPPCApply);CHKERRQ(ierr);
1817   ierr = PCShellSetApplyTranspose(newpc,FETIDPPCApplyTranspose);CHKERRQ(ierr);
1818   ierr = PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);CHKERRQ(ierr);
1819   ierr = PCSetOperators(newpc,newmat,newmat);CHKERRQ(ierr);
1820   ierr = PCSetUp(newpc);CHKERRQ(ierr);
1821   /* return pointers for objects created */
1822   *fetidp_mat=newmat;
1823   *fetidp_pc=newpc;
1824   PetscFunctionReturn(0);
1825 }
1826 
1827 #undef __FUNCT__
1828 #define __FUNCT__ "PCBDDCCreateFETIDPOperators"
1829 /*@
1830  PCBDDCCreateFETIDPOperators - Create FETI-DP operators
1831 
1832    Collective
1833 
1834    Input Parameters:
1835 .  pc - the BDDC preconditioning context (setup should have been called before)
1836 
1837    Output Parameters:
1838 +  fetidp_mat - shell FETI-DP matrix object
1839 -  fetidp_pc  - shell Dirichlet preconditioner for FETI-DP matrix
1840 
1841    Options Database Keys:
1842 .    -fetidp_fullyredundant <false> - use or not a fully redundant set of Lagrange multipliers
1843 
1844    Level: developer
1845 
1846    Notes:
1847      Currently the only operations provided for FETI-DP matrix are MatMult and MatMultTranspose
1848 
1849 .seealso: PCBDDC, PCBDDCMatFETIDPGetRHS, PCBDDCMatFETIDPGetSolution
1850 @*/
1851 PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, Mat *fetidp_mat, PC *fetidp_pc)
1852 {
1853   PetscErrorCode ierr;
1854 
1855   PetscFunctionBegin;
1856   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1857   if (pc->setupcalled) {
1858     ierr = PetscUseMethod(pc,"PCBDDCCreateFETIDPOperators_C",(PC,Mat*,PC*),(pc,fetidp_mat,fetidp_pc));CHKERRQ(ierr);
1859   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"You must call PCSetup_BDDC() first \n");
1860   PetscFunctionReturn(0);
1861 }
1862 /* -------------------------------------------------------------------------- */
1863 /*MC
1864    PCBDDC - Balancing Domain Decomposition by Constraints.
1865 
1866    An implementation of the BDDC preconditioner based on
1867 
1868 .vb
1869    [1] C. R. Dohrmann. "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007
1870    [2] A. Klawonn and O. B. Widlund. "Dual-Primal FETI Methods for Linear Elasticity", http://cs.nyu.edu/csweb/Research/TechReports/TR2004-855/TR2004-855.pdf
1871    [3] J. Mandel, B. Sousedik, C. R. Dohrmann. "Multispace and Multilevel BDDC", http://arxiv.org/abs/0712.3977
1872    [4] C. Pechstein and C. R. Dohrmann. "Modern domain decomposition methods BDDC, deluxe scaling, and an algebraic approach", Seminar talk, Linz, December 2013, http://people.ricam.oeaw.ac.at/c.pechstein/pechstein-bddc2013.pdf
1873 .ve
1874 
1875    The matrix to be preconditioned (Pmat) must be of type MATIS.
1876 
1877    Currently works with MATIS matrices with local matrices of type MATSEQAIJ, MATSEQBAIJ or MATSEQSBAIJ, either with real or complex numbers.
1878 
1879    It also works with unsymmetric and indefinite problems.
1880 
1881    Unlike 'conventional' interface preconditioners, PCBDDC iterates over all degrees of freedom, not just those on the interface. This allows the use of approximate solvers on the subdomains.
1882 
1883    Approximate local solvers are automatically adapted for singular linear problems (see [1]) if the user has provided the nullspace using PCBDDCSetNullSpace()
1884 
1885    Boundary nodes are split in vertices, edges and faces classes using information from the local to global mapping of dofs and the local connectivity graph of nodes. The latter can be customized by using PCBDDCSetLocalAdjacencyGraph()
1886    Additional information on dofs can be provided by using PCBDDCSetDofsSplitting(), PCBDDCSetDirichletBoundaries(), PCBDDCSetNeumannBoundaries(), and PCBDDCSetPrimalVerticesLocalIS()
1887 
1888    Constraints can be customized by attaching a MatNullSpace object to the MATIS matrix via MatSetNearNullSpace(). Non-singular modes are retained via SVD.
1889 
1890    Change of basis is performed similarly to [2] when requested. When more than one constraint is present on a single connected component (i.e. an edge or a face), a robust method based on local QR factorizations is used.
1891    User defined change of basis can be passed to PCBDDC by using PCBDDCSetChangeOfBasisMat()
1892 
1893    The PETSc implementation also supports multilevel BDDC [3]. Coarse grids are partitioned using a MatPartitioning object.
1894 
1895    Adaptive selection of primal constraints [4] is supported for SPD systems with high-contrast in the coefficients if MUMPS is present. Future versions of the code will also consider using MKL_PARDISO or PASTIX.
1896 
1897    An experimental interface to the FETI-DP method is available. FETI-DP operators could be created using PCBDDCCreateFETIDPOperators(). A stand-alone class for the FETI-DP method will be provided in the next releases.
1898    Deluxe scaling is not supported yet for FETI-DP.
1899 
1900    Options Database Keys (some of them, run with -h for a complete list):
1901 
1902 .    -pc_bddc_use_vertices <true> - use or not vertices in primal space
1903 .    -pc_bddc_use_edges <true> - use or not edges in primal space
1904 .    -pc_bddc_use_faces <false> - use or not faces in primal space
1905 .    -pc_bddc_symmetric <true> - symmetric computation of primal basis functions. Specify false for unsymmetric problems
1906 .    -pc_bddc_use_change_of_basis <false> - use change of basis approach (on edges only)
1907 .    -pc_bddc_use_change_on_faces <false> - use change of basis approach on faces if change of basis has been requested
1908 .    -pc_bddc_switch_static <false> - switches from M_2 (default) to M_3 operator (see reference article [1])
1909 .    -pc_bddc_levels <0> - maximum number of levels for multilevel
1910 .    -pc_bddc_coarsening_ratio <8> - number of subdomains which will be aggregated together at the coarser level (e.g. H/h ratio at the coarser level, significative only in the multilevel case)
1911 .    -pc_bddc_redistribute <0> - size of a subset of processors where the coarse problem will be remapped (the value is ignored if not at the coarsest level)
1912 .    -pc_bddc_use_deluxe_scaling <false> - use deluxe scaling
1913 .    -pc_bddc_schur_layers <-1> - select the economic version of deluxe scaling by specifying the number of layers (-1 corresponds to the original deluxe scaling)
1914 .    -pc_bddc_adaptive_threshold <0.0> - when a value greater than one is specified, adaptive selection of constraints is performed on edges and faces (requires deluxe scaling and MUMPS installed)
1915 -    -pc_bddc_check_level <0> - set verbosity level of debugging output
1916 
1917    Options for Dirichlet, Neumann or coarse solver can be set with
1918 .vb
1919       -pc_bddc_dirichlet_
1920       -pc_bddc_neumann_
1921       -pc_bddc_coarse_
1922 .ve
1923    e.g -pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg. PCBDDC uses by default KPSPREONLY and PCLU.
1924 
1925    When using a multilevel approach, solvers' options at the N-th level (N > 1) can be specified as
1926 .vb
1927       -pc_bddc_dirichlet_lN_
1928       -pc_bddc_neumann_lN_
1929       -pc_bddc_coarse_lN_
1930 .ve
1931    Note that level number ranges from the finest (0) to the coarsest (N).
1932    In order to specify options for the BDDC operators at the coarser levels (and not for the solvers), prepend -pc_bddc_coarse_ or -pc_bddc_coarse_l to the option, e.g.
1933 .vb
1934      -pc_bddc_coarse_pc_bddc_adaptive_threshold 5 -pc_bddc_coarse_l1_pc_bddc_redistribute 3
1935 .ve
1936    will use a threshold of 5 for constraints' selection at the first coarse level and will redistribute the coarse problem of the first coarse level on 3 processors
1937 
1938    Level: intermediate
1939 
1940    Developer notes:
1941 
1942    Contributed by Stefano Zampini
1943 
1944 .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,  MATIS
1945 M*/
1946 
1947 #undef __FUNCT__
1948 #define __FUNCT__ "PCCreate_BDDC"
1949 PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc)
1950 {
1951   PetscErrorCode      ierr;
1952   PC_BDDC             *pcbddc;
1953 
1954   PetscFunctionBegin;
1955   /* Creates the private data structure for this preconditioner and attach it to the PC object. */
1956   ierr      = PetscNewLog(pc,&pcbddc);CHKERRQ(ierr);
1957   pc->data  = (void*)pcbddc;
1958 
1959   /* create PCIS data structure */
1960   ierr = PCISCreate(pc);CHKERRQ(ierr);
1961 
1962   /* BDDC customization */
1963   pcbddc->use_local_adj       = PETSC_TRUE;
1964   pcbddc->use_vertices        = PETSC_TRUE;
1965   pcbddc->use_edges           = PETSC_TRUE;
1966   pcbddc->use_faces           = PETSC_FALSE;
1967   pcbddc->use_change_of_basis = PETSC_FALSE;
1968   pcbddc->use_change_on_faces = PETSC_FALSE;
1969   pcbddc->switch_static       = PETSC_FALSE;
1970   pcbddc->use_nnsp_true       = PETSC_FALSE;
1971   pcbddc->use_qr_single       = PETSC_FALSE;
1972   pcbddc->symmetric_primal    = PETSC_TRUE;
1973   pcbddc->dbg_flag            = 0;
1974   /* private */
1975   pcbddc->local_primal_size          = 0;
1976   pcbddc->local_primal_size_cc       = 0;
1977   pcbddc->local_primal_ref_node      = 0;
1978   pcbddc->local_primal_ref_mult      = 0;
1979   pcbddc->n_vertices                 = 0;
1980   pcbddc->primal_indices_local_idxs  = 0;
1981   pcbddc->recompute_topography       = PETSC_FALSE;
1982   pcbddc->coarse_size                = -1;
1983   pcbddc->new_primal_space           = PETSC_FALSE;
1984   pcbddc->new_primal_space_local     = PETSC_FALSE;
1985   pcbddc->global_primal_indices      = 0;
1986   pcbddc->onearnullspace             = 0;
1987   pcbddc->onearnullvecs_state        = 0;
1988   pcbddc->user_primal_vertices       = 0;
1989   pcbddc->NullSpace                  = 0;
1990   pcbddc->temp_solution              = 0;
1991   pcbddc->original_rhs               = 0;
1992   pcbddc->local_mat                  = 0;
1993   pcbddc->ChangeOfBasisMatrix        = 0;
1994   pcbddc->user_ChangeOfBasisMatrix   = 0;
1995   pcbddc->new_global_mat             = 0;
1996   pcbddc->coarse_vec                 = 0;
1997   pcbddc->coarse_ksp                 = 0;
1998   pcbddc->coarse_phi_B               = 0;
1999   pcbddc->coarse_phi_D               = 0;
2000   pcbddc->coarse_psi_B               = 0;
2001   pcbddc->coarse_psi_D               = 0;
2002   pcbddc->vec1_P                     = 0;
2003   pcbddc->vec1_R                     = 0;
2004   pcbddc->vec2_R                     = 0;
2005   pcbddc->local_auxmat1              = 0;
2006   pcbddc->local_auxmat2              = 0;
2007   pcbddc->R_to_B                     = 0;
2008   pcbddc->R_to_D                     = 0;
2009   pcbddc->ksp_D                      = 0;
2010   pcbddc->ksp_R                      = 0;
2011   pcbddc->NeumannBoundaries          = 0;
2012   pcbddc->NeumannBoundariesLocal     = 0;
2013   pcbddc->DirichletBoundaries        = 0;
2014   pcbddc->DirichletBoundariesLocal   = 0;
2015   pcbddc->user_provided_isfordofs    = PETSC_FALSE;
2016   pcbddc->n_ISForDofs                = 0;
2017   pcbddc->n_ISForDofsLocal           = 0;
2018   pcbddc->ISForDofs                  = 0;
2019   pcbddc->ISForDofsLocal             = 0;
2020   pcbddc->ConstraintMatrix           = 0;
2021   pcbddc->use_exact_dirichlet_trick  = PETSC_TRUE;
2022   pcbddc->coarse_loc_to_glob         = 0;
2023   pcbddc->coarsening_ratio           = 8;
2024   pcbddc->coarse_adj_red             = 0;
2025   pcbddc->current_level              = 0;
2026   pcbddc->max_levels                 = 0;
2027   pcbddc->use_coarse_estimates       = PETSC_FALSE;
2028   pcbddc->redistribute_coarse        = 0;
2029   pcbddc->coarse_subassembling       = 0;
2030   pcbddc->coarse_subassembling_init  = 0;
2031 
2032   /* create local graph structure */
2033   ierr = PCBDDCGraphCreate(&pcbddc->mat_graph);CHKERRQ(ierr);
2034 
2035   /* scaling */
2036   pcbddc->work_scaling          = 0;
2037   pcbddc->use_deluxe_scaling    = PETSC_FALSE;
2038   pcbddc->faster_deluxe         = PETSC_FALSE;
2039 
2040   /* create sub schurs structure */
2041   ierr = PCBDDCSubSchursCreate(&pcbddc->sub_schurs);CHKERRQ(ierr);
2042   pcbddc->sub_schurs_rebuild     = PETSC_FALSE;
2043   pcbddc->sub_schurs_layers      = -1;
2044   pcbddc->sub_schurs_use_useradj = PETSC_FALSE;
2045 
2046   pcbddc->computed_rowadj = PETSC_FALSE;
2047 
2048   /* adaptivity */
2049   pcbddc->adaptive_threshold      = 0.0;
2050   pcbddc->adaptive_nmax           = 0;
2051   pcbddc->adaptive_nmin           = 0;
2052 
2053   /* function pointers */
2054   pc->ops->apply               = PCApply_BDDC;
2055   pc->ops->applytranspose      = PCApplyTranspose_BDDC;
2056   pc->ops->setup               = PCSetUp_BDDC;
2057   pc->ops->destroy             = PCDestroy_BDDC;
2058   pc->ops->setfromoptions      = PCSetFromOptions_BDDC;
2059   pc->ops->view                = 0;
2060   pc->ops->applyrichardson     = 0;
2061   pc->ops->applysymmetricleft  = 0;
2062   pc->ops->applysymmetricright = 0;
2063   pc->ops->presolve            = PCPreSolve_BDDC;
2064   pc->ops->postsolve           = PCPostSolve_BDDC;
2065 
2066   /* composing function */
2067   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",PCBDDCSetChangeOfBasisMat_BDDC);CHKERRQ(ierr);
2068   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);CHKERRQ(ierr);
2069   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);CHKERRQ(ierr);
2070   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",PCBDDCSetLevel_BDDC);CHKERRQ(ierr);
2071   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",PCBDDCSetUseExactDirichlet_BDDC);CHKERRQ(ierr);
2072   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",PCBDDCSetLevels_BDDC);CHKERRQ(ierr);
2073   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",PCBDDCSetNullSpace_BDDC);CHKERRQ(ierr);
2074   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",PCBDDCSetDirichletBoundaries_BDDC);CHKERRQ(ierr);
2075   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",PCBDDCSetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr);
2076   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",PCBDDCSetNeumannBoundaries_BDDC);CHKERRQ(ierr);
2077   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",PCBDDCSetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr);
2078   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",PCBDDCGetDirichletBoundaries_BDDC);CHKERRQ(ierr);
2079   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",PCBDDCGetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr);
2080   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",PCBDDCGetNeumannBoundaries_BDDC);CHKERRQ(ierr);
2081   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",PCBDDCGetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr);
2082   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);CHKERRQ(ierr);
2083   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",PCBDDCSetDofsSplittingLocal_BDDC);CHKERRQ(ierr);
2084   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);CHKERRQ(ierr);
2085   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);CHKERRQ(ierr);
2086   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);CHKERRQ(ierr);
2087   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);CHKERRQ(ierr);
2088   PetscFunctionReturn(0);
2089 }
2090 
2091