xref: /petsc/src/ksp/pc/impls/bddc/bddc.c (revision a7fac7c2b47dbcd8ece0cc2dfdfe0e63be1bb7b5)
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->ctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1002       ierr = VecScatterEnd(matis->ctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1003       ierr = VecScatterBegin(matis->ctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1004       ierr = VecScatterEnd(matis->ctx,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->ctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1014       ierr = VecScatterEnd(matis->ctx,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   PetscBool      computetopography,computesolvers,computesubschurs;
1176   PetscBool      computeconstraintsmatrix;
1177   PetscBool      new_nearnullspace_provided,ismatis;
1178 
1179   PetscFunctionBegin;
1180   ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATIS,&ismatis);CHKERRQ(ierr);
1181   if (!ismatis) {
1182     SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PCBDDC preconditioner requires matrix of type MATIS");
1183   }
1184   matis = (Mat_IS*)pc->pmat->data;
1185   /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other future nonoverlapping preconditioners */
1186   /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup
1187      Also, BDDC directly build the Dirichlet problem */
1188   /* split work */
1189   if (pc->setupcalled) {
1190     if (pc->flag == SAME_NONZERO_PATTERN) {
1191       computetopography = PETSC_FALSE;
1192       computesolvers = PETSC_TRUE;
1193     } else { /* DIFFERENT_NONZERO_PATTERN */
1194       computetopography = PETSC_TRUE;
1195       computesolvers = PETSC_TRUE;
1196     }
1197   } else {
1198     computetopography = PETSC_TRUE;
1199     computesolvers = PETSC_TRUE;
1200   }
1201   if (pcbddc->recompute_topography) {
1202     computetopography = PETSC_TRUE;
1203   }
1204   computeconstraintsmatrix = PETSC_FALSE;
1205   if (pcbddc->adaptive_threshold > 0.0 && !pcbddc->use_deluxe_scaling) {
1206     SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot compute adaptive constraints without deluxe scaling. Rerun with -pc_bddc_use_deluxe_scaling");
1207   }
1208   pcbddc->adaptive_selection = (PetscBool)(pcbddc->adaptive_threshold > 0.0 && pcbddc->use_deluxe_scaling);
1209   if (pcbddc->adaptive_selection) pcbddc->use_faces = PETSC_TRUE;
1210 
1211   computesubschurs = (PetscBool)(pcbddc->adaptive_selection || pcbddc->use_deluxe_scaling);
1212   if (pcbddc->faster_deluxe && pcbddc->adaptive_selection && pcbddc->use_change_of_basis) {
1213     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");
1214   }
1215   /* Get stdout for dbg */
1216   if (pcbddc->dbg_flag) {
1217     if (!pcbddc->dbg_viewer) {
1218       pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc));
1219       ierr = PetscViewerASCIISynchronizedAllow(pcbddc->dbg_viewer,PETSC_TRUE);CHKERRQ(ierr);
1220     }
1221     ierr = PetscViewerASCIIAddTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr);
1222   }
1223 
1224   if (pcbddc->user_ChangeOfBasisMatrix) {
1225     /* use_change_of_basis flag is used to automatically compute a change of basis from constraints */
1226     pcbddc->use_change_of_basis = PETSC_FALSE;
1227     ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr);
1228   } else {
1229     ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr);
1230     ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr);
1231     pcbddc->local_mat = matis->A;
1232   }
1233 
1234   /* workaround for reals */
1235 #if !defined(PETSC_USE_COMPLEX)
1236   if (matis->A->symmetric_set) {
1237     ierr = MatSetOption(pcbddc->local_mat,MAT_HERMITIAN,matis->A->symmetric);CHKERRQ(ierr);
1238   }
1239 #endif
1240 
1241   /* Set up all the "iterative substructuring" common block without computing solvers */
1242   {
1243     Mat temp_mat;
1244 
1245     temp_mat = matis->A;
1246     matis->A = pcbddc->local_mat;
1247     ierr = PCISSetUp(pc,PETSC_FALSE);CHKERRQ(ierr);
1248     pcbddc->local_mat = matis->A;
1249     matis->A = temp_mat;
1250   }
1251 
1252   /* Analyze interface and setup sub_schurs data */
1253   if (computetopography) {
1254     ierr = PCBDDCAnalyzeInterface(pc);CHKERRQ(ierr);
1255     computeconstraintsmatrix = PETSC_TRUE;
1256   }
1257 
1258   /* Setup local dirichlet solver ksp_D and sub_schurs solvers */
1259   if (computesolvers) {
1260     PCBDDCSubSchurs sub_schurs=pcbddc->sub_schurs;
1261 
1262     if (computesubschurs && computetopography) {
1263       ierr = PCBDDCInitSubSchurs(pc);CHKERRQ(ierr);
1264     }
1265     if (sub_schurs->use_mumps) {
1266       if (computesubschurs) {
1267         ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr);
1268       }
1269       ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr);
1270     } else {
1271       ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr);
1272       if (computesubschurs) {
1273         ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr);
1274       }
1275     }
1276     if (pcbddc->adaptive_selection) {
1277       ierr = PCBDDCAdaptiveSelection(pc);CHKERRQ(ierr);
1278       computeconstraintsmatrix = PETSC_TRUE;
1279     }
1280   }
1281 
1282   /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */
1283   new_nearnullspace_provided = PETSC_FALSE;
1284   ierr = MatGetNearNullSpace(pc->pmat,&nearnullspace);CHKERRQ(ierr);
1285   if (pcbddc->onearnullspace) { /* already used nearnullspace */
1286     if (!nearnullspace) { /* near null space attached to mat has been destroyed */
1287       new_nearnullspace_provided = PETSC_TRUE;
1288     } else {
1289       /* determine if the two nullspaces are different (should be lightweight) */
1290       if (nearnullspace != pcbddc->onearnullspace) {
1291         new_nearnullspace_provided = PETSC_TRUE;
1292       } else { /* maybe the user has changed the content of the nearnullspace so check vectors ObjectStateId */
1293         PetscInt         i;
1294         const Vec        *nearnullvecs;
1295         PetscObjectState state;
1296         PetscInt         nnsp_size;
1297         ierr = MatNullSpaceGetVecs(nearnullspace,NULL,&nnsp_size,&nearnullvecs);CHKERRQ(ierr);
1298         for (i=0;i<nnsp_size;i++) {
1299           ierr = PetscObjectStateGet((PetscObject)nearnullvecs[i],&state);CHKERRQ(ierr);
1300           if (pcbddc->onearnullvecs_state[i] != state) {
1301             new_nearnullspace_provided = PETSC_TRUE;
1302             break;
1303           }
1304         }
1305       }
1306     }
1307   } else {
1308     if (!nearnullspace) { /* both nearnullspaces are null */
1309       new_nearnullspace_provided = PETSC_FALSE;
1310     } else { /* nearnullspace attached later */
1311       new_nearnullspace_provided = PETSC_TRUE;
1312     }
1313   }
1314 
1315   /* Setup constraints and related work vectors */
1316   /* reset primal space flags */
1317   pcbddc->new_primal_space = PETSC_FALSE;
1318   pcbddc->new_primal_space_local = PETSC_FALSE;
1319   if (computeconstraintsmatrix || new_nearnullspace_provided) {
1320     /* It also sets the primal space flags */
1321     ierr = PCBDDCConstraintsSetUp(pc);CHKERRQ(ierr);
1322     /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */
1323     ierr = PCBDDCSetUpLocalWorkVectors(pc);CHKERRQ(ierr);
1324   }
1325 
1326   if (computesolvers || pcbddc->new_primal_space) {
1327     if (pcbddc->use_change_of_basis) {
1328       PC_IS *pcis = (PC_IS*)(pc->data);
1329 
1330       ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr);
1331       /* get submatrices */
1332       ierr = MatDestroy(&pcis->A_IB);CHKERRQ(ierr);
1333       ierr = MatDestroy(&pcis->A_BI);CHKERRQ(ierr);
1334       ierr = MatDestroy(&pcis->A_BB);CHKERRQ(ierr);
1335       ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_BB);CHKERRQ(ierr);
1336       ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_I_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_IB);CHKERRQ(ierr);
1337       ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_BI);CHKERRQ(ierr);
1338       /* set flag in pcis to not reuse submatrices during PCISCreate */
1339       pcis->reusesubmatrices = PETSC_FALSE;
1340     } else if (!pcbddc->user_ChangeOfBasisMatrix) {
1341       ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr);
1342       ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr);
1343       pcbddc->local_mat = matis->A;
1344     }
1345     /* SetUp coarse and local Neumann solvers */
1346     ierr = PCBDDCSetUpSolvers(pc);CHKERRQ(ierr);
1347     /* SetUp Scaling operator */
1348     ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr);
1349   }
1350 
1351   if (pcbddc->dbg_flag) {
1352     ierr = PetscViewerASCIISubtractTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr);
1353   }
1354   PetscFunctionReturn(0);
1355 }
1356 
1357 /* -------------------------------------------------------------------------- */
1358 /*
1359    PCApply_BDDC - Applies the BDDC operator to a vector.
1360 
1361    Input Parameters:
1362 +  pc - the preconditioner context
1363 -  r - input vector (global)
1364 
1365    Output Parameter:
1366 .  z - output vector (global)
1367 
1368    Application Interface Routine: PCApply()
1369  */
1370 #undef __FUNCT__
1371 #define __FUNCT__ "PCApply_BDDC"
1372 PetscErrorCode PCApply_BDDC(PC pc,Vec r,Vec z)
1373 {
1374   PC_IS             *pcis = (PC_IS*)(pc->data);
1375   PC_BDDC           *pcbddc = (PC_BDDC*)(pc->data);
1376   PetscInt          n_B = pcis->n_B, n_D = pcis->n - n_B;
1377   PetscErrorCode    ierr;
1378   const PetscScalar one = 1.0;
1379   const PetscScalar m_one = -1.0;
1380   const PetscScalar zero = 0.0;
1381 
1382 /* This code is similar to that provided in nn.c for PCNN
1383    NN interface preconditioner changed to BDDC
1384    Added support for M_3 preconditioner in the reference article (code is active if pcbddc->switch_static == PETSC_TRUE) */
1385 
1386   PetscFunctionBegin;
1387   if (!pcbddc->use_exact_dirichlet_trick) {
1388     ierr = VecCopy(r,z);CHKERRQ(ierr);
1389     /* First Dirichlet solve */
1390     ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1391     ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1392     /*
1393       Assembling right hand side for BDDC operator
1394       - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE)
1395       - pcis->vec1_B the interface part of the global vector z
1396     */
1397     if (n_D) {
1398       ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr);
1399       ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr);
1400       if (pcbddc->switch_static) { ierr = MatMultAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); }
1401       ierr = MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr);
1402     } else {
1403       ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr);
1404     }
1405     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1406     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1407     ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr);
1408   } else {
1409     if (pcbddc->switch_static) {
1410       ierr = VecSet(pcis->vec1_D,zero);CHKERRQ(ierr);
1411     }
1412     ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr);
1413   }
1414 
1415   /* Apply interface preconditioner
1416      input/output vecs: pcis->vec1_B and pcis->vec1_D */
1417   ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_FALSE);CHKERRQ(ierr);
1418 
1419   /* Apply transpose of partition of unity operator */
1420   ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr);
1421 
1422   /* Second Dirichlet solve and assembling of output */
1423   ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1424   ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1425   if (n_B) {
1426     ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr);
1427     if (pcbddc->switch_static) { ierr = MatMultAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D);CHKERRQ(ierr); }
1428   } else if (pcbddc->switch_static) {
1429     ierr = MatMult(pcis->A_II,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr);
1430   }
1431   ierr = KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr);
1432   if (!pcbddc->use_exact_dirichlet_trick) {
1433     if (pcbddc->switch_static) {
1434       ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr);
1435     } else {
1436       ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr);
1437     }
1438     ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1439     ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1440   } else {
1441     if (pcbddc->switch_static) {
1442       ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr);
1443     } else {
1444       ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr);
1445     }
1446     ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1447     ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1448   }
1449   PetscFunctionReturn(0);
1450 }
1451 
1452 /* -------------------------------------------------------------------------- */
1453 /*
1454    PCApplyTranspose_BDDC - Applies the transpose of the BDDC operator to a vector.
1455 
1456    Input Parameters:
1457 +  pc - the preconditioner context
1458 -  r - input vector (global)
1459 
1460    Output Parameter:
1461 .  z - output vector (global)
1462 
1463    Application Interface Routine: PCApplyTranspose()
1464  */
1465 #undef __FUNCT__
1466 #define __FUNCT__ "PCApplyTranspose_BDDC"
1467 PetscErrorCode PCApplyTranspose_BDDC(PC pc,Vec r,Vec z)
1468 {
1469   PC_IS             *pcis = (PC_IS*)(pc->data);
1470   PC_BDDC           *pcbddc = (PC_BDDC*)(pc->data);
1471   PetscInt          n_B = pcis->n_B, n_D = pcis->n - n_B;
1472   PetscErrorCode    ierr;
1473   const PetscScalar one = 1.0;
1474   const PetscScalar m_one = -1.0;
1475   const PetscScalar zero = 0.0;
1476 
1477   PetscFunctionBegin;
1478   if (!pcbddc->use_exact_dirichlet_trick) {
1479     ierr = VecCopy(r,z);CHKERRQ(ierr);
1480     /* First Dirichlet solve */
1481     ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1482     ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1483     /*
1484       Assembling right hand side for BDDC operator
1485       - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE)
1486       - pcis->vec1_B the interface part of the global vector z
1487     */
1488     if (n_D) {
1489       ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr);
1490       ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr);
1491       if (pcbddc->switch_static) { ierr = MatMultTransposeAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); }
1492       ierr = MatMultTranspose(pcis->A_IB,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr);
1493     } else {
1494       ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr);
1495     }
1496     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1497     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1498     ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr);
1499   } else {
1500     if (pcbddc->switch_static) {
1501       ierr = VecSet(pcis->vec1_D,zero);CHKERRQ(ierr);
1502     }
1503     ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr);
1504   }
1505 
1506   /* Apply interface preconditioner
1507      input/output vecs: pcis->vec1_B and pcis->vec1_D */
1508   ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_TRUE);CHKERRQ(ierr);
1509 
1510   /* Apply transpose of partition of unity operator */
1511   ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr);
1512 
1513   /* Second Dirichlet solve and assembling of output */
1514   ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1515   ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1516   if (n_B) {
1517     ierr = MatMultTranspose(pcis->A_BI,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr);
1518     if (pcbddc->switch_static) { ierr = MatMultTransposeAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D);CHKERRQ(ierr); }
1519   } else if (pcbddc->switch_static) {
1520     ierr = MatMultTranspose(pcis->A_II,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr);
1521   }
1522   ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr);
1523   if (!pcbddc->use_exact_dirichlet_trick) {
1524     if (pcbddc->switch_static) {
1525       ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr);
1526     } else {
1527       ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr);
1528     }
1529     ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1530     ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1531   } else {
1532     if (pcbddc->switch_static) {
1533       ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr);
1534     } else {
1535       ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr);
1536     }
1537     ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1538     ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1539   }
1540   PetscFunctionReturn(0);
1541 }
1542 /* -------------------------------------------------------------------------- */
1543 
1544 #undef __FUNCT__
1545 #define __FUNCT__ "PCDestroy_BDDC"
1546 PetscErrorCode PCDestroy_BDDC(PC pc)
1547 {
1548   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
1549   PetscErrorCode ierr;
1550 
1551   PetscFunctionBegin;
1552   /* free data created by PCIS */
1553   ierr = PCISDestroy(pc);CHKERRQ(ierr);
1554   /* free BDDC custom data  */
1555   ierr = PCBDDCResetCustomization(pc);CHKERRQ(ierr);
1556   /* destroy objects related to topography */
1557   ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr);
1558   /* free allocated graph structure */
1559   ierr = PetscFree(pcbddc->mat_graph);CHKERRQ(ierr);
1560   /* free allocated sub schurs structure */
1561   ierr = PetscFree(pcbddc->sub_schurs);CHKERRQ(ierr);
1562   /* destroy objects for scaling operator */
1563   ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr);
1564   ierr = PetscFree(pcbddc->deluxe_ctx);CHKERRQ(ierr);
1565   /* free solvers stuff */
1566   ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr);
1567   /* free global vectors needed in presolve */
1568   ierr = VecDestroy(&pcbddc->temp_solution);CHKERRQ(ierr);
1569   ierr = VecDestroy(&pcbddc->original_rhs);CHKERRQ(ierr);
1570   /* free stuff for change of basis hooks */
1571   if (pcbddc->new_global_mat) {
1572     PCBDDCChange_ctx change_ctx;
1573     ierr = MatShellGetContext(pcbddc->new_global_mat,&change_ctx);CHKERRQ(ierr);
1574     ierr = MatDestroy(&change_ctx->original_mat);CHKERRQ(ierr);
1575     ierr = MatDestroy(&change_ctx->global_change);CHKERRQ(ierr);
1576     ierr = VecDestroyVecs(2,&change_ctx->work);CHKERRQ(ierr);
1577     ierr = PetscFree(change_ctx);CHKERRQ(ierr);
1578   }
1579   ierr = MatDestroy(&pcbddc->new_global_mat);CHKERRQ(ierr);
1580   /* remove functions */
1581   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",NULL);CHKERRQ(ierr);
1582   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);CHKERRQ(ierr);
1583   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);CHKERRQ(ierr);
1584   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",NULL);CHKERRQ(ierr);
1585   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",NULL);CHKERRQ(ierr);
1586   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",NULL);CHKERRQ(ierr);
1587   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",NULL);CHKERRQ(ierr);
1588   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",NULL);CHKERRQ(ierr);
1589   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr);
1590   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",NULL);CHKERRQ(ierr);
1591   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr);
1592   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",NULL);CHKERRQ(ierr);
1593   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr);
1594   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",NULL);CHKERRQ(ierr);
1595   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr);
1596   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);CHKERRQ(ierr);
1597   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",NULL);CHKERRQ(ierr);
1598   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);CHKERRQ(ierr);
1599   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);CHKERRQ(ierr);
1600   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);CHKERRQ(ierr);
1601   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);CHKERRQ(ierr);
1602   /* Free the private data structure */
1603   ierr = PetscFree(pc->data);CHKERRQ(ierr);
1604   PetscFunctionReturn(0);
1605 }
1606 /* -------------------------------------------------------------------------- */
1607 
1608 #undef __FUNCT__
1609 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS_BDDC"
1610 static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
1611 {
1612   FETIDPMat_ctx  mat_ctx;
1613   Vec            copy_standard_rhs;
1614   PC_IS*         pcis;
1615   PC_BDDC*       pcbddc;
1616   PetscErrorCode ierr;
1617 
1618   PetscFunctionBegin;
1619   ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr);
1620   pcis = (PC_IS*)mat_ctx->pc->data;
1621   pcbddc = (PC_BDDC*)mat_ctx->pc->data;
1622 
1623   /*
1624      change of basis for physical rhs if needed
1625      It also changes the rhs in case of dirichlet boundaries
1626      TODO: better management when FETIDP will have its own class
1627   */
1628   ierr = VecDuplicate(standard_rhs,&copy_standard_rhs);CHKERRQ(ierr);
1629   ierr = VecCopy(standard_rhs,copy_standard_rhs);CHKERRQ(ierr);
1630   ierr = PCPreSolve_BDDC(mat_ctx->pc,NULL,copy_standard_rhs,NULL);CHKERRQ(ierr);
1631   /* store vectors for computation of fetidp final solution */
1632   ierr = VecScatterBegin(pcis->global_to_D,copy_standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1633   ierr = VecScatterEnd(pcis->global_to_D,copy_standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1634   /* scale rhs since it should be unassembled */
1635   /* TODO use counter scaling? (also below) */
1636   ierr = VecScatterBegin(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1637   ierr = VecScatterEnd(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1638   /* Apply partition of unity */
1639   ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr);
1640   /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,copy_standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */
1641   if (!pcbddc->switch_static) {
1642     /* compute partially subassembled Schur complement right-hand side */
1643     ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr);
1644     ierr = MatMult(pcis->A_BI,pcis->vec1_D,pcis->vec1_B);CHKERRQ(ierr);
1645     ierr = VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);CHKERRQ(ierr);
1646     ierr = VecSet(copy_standard_rhs,0.0);CHKERRQ(ierr);
1647     ierr = VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,copy_standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1648     ierr = VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,copy_standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1649     /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,copy_standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */
1650     ierr = VecScatterBegin(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1651     ierr = VecScatterEnd(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1652     ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr);
1653   }
1654   ierr = VecDestroy(&copy_standard_rhs);CHKERRQ(ierr);
1655   /* BDDC rhs */
1656   ierr = VecCopy(mat_ctx->temp_solution_B,pcis->vec1_B);CHKERRQ(ierr);
1657   if (pcbddc->switch_static) {
1658     ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr);
1659   }
1660   /* apply BDDC */
1661   ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr);
1662   /* Application of B_delta and assembling of rhs for fetidp fluxes */
1663   ierr = VecSet(fetidp_flux_rhs,0.0);CHKERRQ(ierr);
1664   ierr = MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);CHKERRQ(ierr);
1665   ierr = VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1666   ierr = VecScatterEnd(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1667   PetscFunctionReturn(0);
1668 }
1669 
1670 #undef __FUNCT__
1671 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS"
1672 /*@
1673  PCBDDCMatFETIDPGetRHS - Compute the right-hand side for FETI-DP linear system using the physical right-hand side
1674 
1675    Collective
1676 
1677    Input Parameters:
1678 +  fetidp_mat      - the FETI-DP matrix object obtained by a call to PCBDDCCreateFETIDPOperators
1679 -  standard_rhs    - the right-hand side of the original linear system
1680 
1681    Output Parameters:
1682 .  fetidp_flux_rhs - the right-hand side for the FETI-DP linear system
1683 
1684    Level: developer
1685 
1686    Notes:
1687 
1688 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetSolution
1689 @*/
1690 PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
1691 {
1692   FETIDPMat_ctx  mat_ctx;
1693   PetscErrorCode ierr;
1694 
1695   PetscFunctionBegin;
1696   ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr);
1697   ierr = PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetRHS_C",(Mat,Vec,Vec),(fetidp_mat,standard_rhs,fetidp_flux_rhs));CHKERRQ(ierr);
1698   PetscFunctionReturn(0);
1699 }
1700 /* -------------------------------------------------------------------------- */
1701 
1702 #undef __FUNCT__
1703 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution_BDDC"
1704 static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
1705 {
1706   FETIDPMat_ctx  mat_ctx;
1707   PC_IS*         pcis;
1708   PC_BDDC*       pcbddc;
1709   PetscErrorCode ierr;
1710 
1711   PetscFunctionBegin;
1712   ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr);
1713   pcis = (PC_IS*)mat_ctx->pc->data;
1714   pcbddc = (PC_BDDC*)mat_ctx->pc->data;
1715 
1716   /* apply B_delta^T */
1717   ierr = VecScatterBegin(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1718   ierr = VecScatterEnd  (mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1719   ierr = MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr);
1720   /* compute rhs for BDDC application */
1721   ierr = VecAYPX(pcis->vec1_B,-1.0,mat_ctx->temp_solution_B);CHKERRQ(ierr);
1722   if (pcbddc->switch_static) {
1723     ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr);
1724   }
1725   /* apply BDDC */
1726   ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr);
1727   /* put values into standard global vector */
1728   ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1729   ierr = VecScatterEnd  (pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1730   if (!pcbddc->switch_static) {
1731     /* compute values into the interior if solved for the partially subassembled Schur complement */
1732     ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);CHKERRQ(ierr);
1733     ierr = VecAXPY(mat_ctx->temp_solution_D,-1.0,pcis->vec1_D);CHKERRQ(ierr);
1734     ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr);
1735   }
1736   ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1737   ierr = VecScatterEnd  (pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1738   /* final change of basis if needed
1739      Is also sums the dirichlet part removed during RHS assembling */
1740   ierr = PCPostSolve_BDDC(mat_ctx->pc,NULL,NULL,standard_sol);CHKERRQ(ierr);
1741   PetscFunctionReturn(0);
1742 }
1743 
1744 #undef __FUNCT__
1745 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution"
1746 /*@
1747  PCBDDCMatFETIDPGetSolution - Compute the physical solution using the solution of the FETI-DP linear system
1748 
1749    Collective
1750 
1751    Input Parameters:
1752 +  fetidp_mat      - the FETI-DP matrix obtained by a call to PCBDDCCreateFETIDPOperators
1753 -  fetidp_flux_sol - the solution of the FETI-DP linear system
1754 
1755    Output Parameters:
1756 .  standard_sol    - the solution defined on the physical domain
1757 
1758    Level: developer
1759 
1760    Notes:
1761 
1762 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetRHS
1763 @*/
1764 PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
1765 {
1766   FETIDPMat_ctx  mat_ctx;
1767   PetscErrorCode ierr;
1768 
1769   PetscFunctionBegin;
1770   ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr);
1771   ierr = PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetSolution_C",(Mat,Vec,Vec),(fetidp_mat,fetidp_flux_sol,standard_sol));CHKERRQ(ierr);
1772   PetscFunctionReturn(0);
1773 }
1774 /* -------------------------------------------------------------------------- */
1775 
1776 extern PetscErrorCode FETIDPMatMult(Mat,Vec,Vec);
1777 extern PetscErrorCode FETIDPMatMultTranspose(Mat,Vec,Vec);
1778 extern PetscErrorCode PCBDDCDestroyFETIDPMat(Mat);
1779 extern PetscErrorCode FETIDPPCApply(PC,Vec,Vec);
1780 extern PetscErrorCode FETIDPPCApplyTranspose(PC,Vec,Vec);
1781 extern PetscErrorCode PCBDDCDestroyFETIDPPC(PC);
1782 
1783 #undef __FUNCT__
1784 #define __FUNCT__ "PCBDDCCreateFETIDPOperators_BDDC"
1785 static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, Mat *fetidp_mat, PC *fetidp_pc)
1786 {
1787 
1788   FETIDPMat_ctx  fetidpmat_ctx;
1789   Mat            newmat;
1790   FETIDPPC_ctx   fetidppc_ctx;
1791   PC             newpc;
1792   MPI_Comm       comm;
1793   PetscErrorCode ierr;
1794 
1795   PetscFunctionBegin;
1796   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
1797   /* FETIDP linear matrix */
1798   ierr = PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);CHKERRQ(ierr);
1799   ierr = PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);CHKERRQ(ierr);
1800   ierr = MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,fetidpmat_ctx->n_lambda,fetidpmat_ctx->n_lambda,fetidpmat_ctx,&newmat);CHKERRQ(ierr);
1801   ierr = MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);CHKERRQ(ierr);
1802   ierr = MatShellSetOperation(newmat,MATOP_MULT_TRANSPOSE,(void (*)(void))FETIDPMatMultTranspose);CHKERRQ(ierr);
1803   ierr = MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);CHKERRQ(ierr);
1804   ierr = MatSetUp(newmat);CHKERRQ(ierr);
1805   /* FETIDP preconditioner */
1806   ierr = PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);CHKERRQ(ierr);
1807   ierr = PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);CHKERRQ(ierr);
1808   ierr = PCCreate(comm,&newpc);CHKERRQ(ierr);
1809   ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr);
1810   ierr = PCShellSetContext(newpc,fetidppc_ctx);CHKERRQ(ierr);
1811   ierr = PCShellSetApply(newpc,FETIDPPCApply);CHKERRQ(ierr);
1812   ierr = PCShellSetApplyTranspose(newpc,FETIDPPCApplyTranspose);CHKERRQ(ierr);
1813   ierr = PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);CHKERRQ(ierr);
1814   ierr = PCSetOperators(newpc,newmat,newmat);CHKERRQ(ierr);
1815   ierr = PCSetUp(newpc);CHKERRQ(ierr);
1816   /* return pointers for objects created */
1817   *fetidp_mat=newmat;
1818   *fetidp_pc=newpc;
1819   PetscFunctionReturn(0);
1820 }
1821 
1822 #undef __FUNCT__
1823 #define __FUNCT__ "PCBDDCCreateFETIDPOperators"
1824 /*@
1825  PCBDDCCreateFETIDPOperators - Create FETI-DP operators
1826 
1827    Collective
1828 
1829    Input Parameters:
1830 .  pc - the BDDC preconditioning context (setup should have been called before)
1831 
1832    Output Parameters:
1833 +  fetidp_mat - shell FETI-DP matrix object
1834 -  fetidp_pc  - shell Dirichlet preconditioner for FETI-DP matrix
1835 
1836    Options Database Keys:
1837 .    -fetidp_fullyredundant <false> - use or not a fully redundant set of Lagrange multipliers
1838 
1839    Level: developer
1840 
1841    Notes:
1842      Currently the only operations provided for FETI-DP matrix are MatMult and MatMultTranspose
1843 
1844 .seealso: PCBDDC, PCBDDCMatFETIDPGetRHS, PCBDDCMatFETIDPGetSolution
1845 @*/
1846 PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, Mat *fetidp_mat, PC *fetidp_pc)
1847 {
1848   PetscErrorCode ierr;
1849 
1850   PetscFunctionBegin;
1851   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1852   if (pc->setupcalled) {
1853     ierr = PetscUseMethod(pc,"PCBDDCCreateFETIDPOperators_C",(PC,Mat*,PC*),(pc,fetidp_mat,fetidp_pc));CHKERRQ(ierr);
1854   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"You must call PCSetup_BDDC() first \n");
1855   PetscFunctionReturn(0);
1856 }
1857 /* -------------------------------------------------------------------------- */
1858 /*MC
1859    PCBDDC - Balancing Domain Decomposition by Constraints.
1860 
1861    An implementation of the BDDC preconditioner based on
1862 
1863 .vb
1864    [1] C. R. Dohrmann. "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007
1865    [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
1866    [3] J. Mandel, B. Sousedik, C. R. Dohrmann. "Multispace and Multilevel BDDC", http://arxiv.org/abs/0712.3977
1867    [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
1868 .ve
1869 
1870    The matrix to be preconditioned (Pmat) must be of type MATIS.
1871 
1872    Currently works with MATIS matrices with local matrices of type MATSEQAIJ, MATSEQBAIJ or MATSEQSBAIJ, either with real or complex numbers.
1873 
1874    It also works with unsymmetric and indefinite problems.
1875 
1876    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.
1877 
1878    Approximate local solvers are automatically adapted for singular linear problems (see [1]) if the user has provided the nullspace using PCBDDCSetNullSpace()
1879 
1880    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()
1881    Additional information on dofs can be provided by using PCBDDCSetDofsSplitting(), PCBDDCSetDirichletBoundaries(), PCBDDCSetNeumannBoundaries(), and PCBDDCSetPrimalVerticesLocalIS()
1882 
1883    Constraints can be customized by attaching a MatNullSpace object to the MATIS matrix via MatSetNearNullSpace(). Non-singular modes are retained via SVD.
1884 
1885    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.
1886    User defined change of basis can be passed to PCBDDC by using PCBDDCSetChangeOfBasisMat()
1887 
1888    The PETSc implementation also supports multilevel BDDC [3]. Coarse grids are partitioned using a MatPartitioning object.
1889 
1890    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.
1891 
1892    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.
1893    Deluxe scaling is not supported yet for FETI-DP.
1894 
1895    Options Database Keys (some of them, run with -h for a complete list):
1896 
1897 .    -pc_bddc_use_vertices <true> - use or not vertices in primal space
1898 .    -pc_bddc_use_edges <true> - use or not edges in primal space
1899 .    -pc_bddc_use_faces <false> - use or not faces in primal space
1900 .    -pc_bddc_symmetric <true> - symmetric computation of primal basis functions. Specify false for unsymmetric problems
1901 .    -pc_bddc_use_change_of_basis <false> - use change of basis approach (on edges only)
1902 .    -pc_bddc_use_change_on_faces <false> - use change of basis approach on faces if change of basis has been requested
1903 .    -pc_bddc_switch_static <false> - switches from M_2 (default) to M_3 operator (see reference article [1])
1904 .    -pc_bddc_levels <0> - maximum number of levels for multilevel
1905 .    -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)
1906 .    -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)
1907 .    -pc_bddc_use_deluxe_scaling <false> - use deluxe scaling
1908 .    -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)
1909 .    -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)
1910 -    -pc_bddc_check_level <0> - set verbosity level of debugging output
1911 
1912    Options for Dirichlet, Neumann or coarse solver can be set with
1913 .vb
1914       -pc_bddc_dirichlet_
1915       -pc_bddc_neumann_
1916       -pc_bddc_coarse_
1917 .ve
1918    e.g -pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg. PCBDDC uses by default KPSPREONLY and PCLU.
1919 
1920    When using a multilevel approach, solvers' options at the N-th level (N > 1) can be specified as
1921 .vb
1922       -pc_bddc_dirichlet_lN_
1923       -pc_bddc_neumann_lN_
1924       -pc_bddc_coarse_lN_
1925 .ve
1926    Note that level number ranges from the finest (0) to the coarsest (N).
1927    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.
1928 .vb
1929      -pc_bddc_coarse_pc_bddc_adaptive_threshold 5 -pc_bddc_coarse_l1_pc_bddc_redistribute 3
1930 .ve
1931    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
1932 
1933    Level: intermediate
1934 
1935    Developer notes:
1936 
1937    Contributed by Stefano Zampini
1938 
1939 .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,  MATIS
1940 M*/
1941 
1942 #undef __FUNCT__
1943 #define __FUNCT__ "PCCreate_BDDC"
1944 PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc)
1945 {
1946   PetscErrorCode      ierr;
1947   PC_BDDC             *pcbddc;
1948 
1949   PetscFunctionBegin;
1950   /* Creates the private data structure for this preconditioner and attach it to the PC object. */
1951   ierr      = PetscNewLog(pc,&pcbddc);CHKERRQ(ierr);
1952   pc->data  = (void*)pcbddc;
1953 
1954   /* create PCIS data structure */
1955   ierr = PCISCreate(pc);CHKERRQ(ierr);
1956 
1957   /* BDDC customization */
1958   pcbddc->use_local_adj       = PETSC_TRUE;
1959   pcbddc->use_vertices        = PETSC_TRUE;
1960   pcbddc->use_edges           = PETSC_TRUE;
1961   pcbddc->use_faces           = PETSC_FALSE;
1962   pcbddc->use_change_of_basis = PETSC_FALSE;
1963   pcbddc->use_change_on_faces = PETSC_FALSE;
1964   pcbddc->switch_static       = PETSC_FALSE;
1965   pcbddc->use_nnsp_true       = PETSC_FALSE;
1966   pcbddc->use_qr_single       = PETSC_FALSE;
1967   pcbddc->symmetric_primal    = PETSC_TRUE;
1968   pcbddc->dbg_flag            = 0;
1969   /* private */
1970   pcbddc->local_primal_size          = 0;
1971   pcbddc->local_primal_size_cc       = 0;
1972   pcbddc->local_primal_ref_node      = 0;
1973   pcbddc->local_primal_ref_mult      = 0;
1974   pcbddc->n_vertices                 = 0;
1975   pcbddc->primal_indices_local_idxs  = 0;
1976   pcbddc->recompute_topography       = PETSC_FALSE;
1977   pcbddc->coarse_size                = -1;
1978   pcbddc->new_primal_space           = PETSC_FALSE;
1979   pcbddc->new_primal_space_local     = PETSC_FALSE;
1980   pcbddc->global_primal_indices      = 0;
1981   pcbddc->onearnullspace             = 0;
1982   pcbddc->onearnullvecs_state        = 0;
1983   pcbddc->user_primal_vertices       = 0;
1984   pcbddc->NullSpace                  = 0;
1985   pcbddc->temp_solution              = 0;
1986   pcbddc->original_rhs               = 0;
1987   pcbddc->local_mat                  = 0;
1988   pcbddc->ChangeOfBasisMatrix        = 0;
1989   pcbddc->user_ChangeOfBasisMatrix   = 0;
1990   pcbddc->new_global_mat             = 0;
1991   pcbddc->coarse_vec                 = 0;
1992   pcbddc->coarse_ksp                 = 0;
1993   pcbddc->coarse_phi_B               = 0;
1994   pcbddc->coarse_phi_D               = 0;
1995   pcbddc->coarse_psi_B               = 0;
1996   pcbddc->coarse_psi_D               = 0;
1997   pcbddc->vec1_P                     = 0;
1998   pcbddc->vec1_R                     = 0;
1999   pcbddc->vec2_R                     = 0;
2000   pcbddc->local_auxmat1              = 0;
2001   pcbddc->local_auxmat2              = 0;
2002   pcbddc->R_to_B                     = 0;
2003   pcbddc->R_to_D                     = 0;
2004   pcbddc->ksp_D                      = 0;
2005   pcbddc->ksp_R                      = 0;
2006   pcbddc->NeumannBoundaries          = 0;
2007   pcbddc->NeumannBoundariesLocal     = 0;
2008   pcbddc->DirichletBoundaries        = 0;
2009   pcbddc->DirichletBoundariesLocal   = 0;
2010   pcbddc->user_provided_isfordofs    = PETSC_FALSE;
2011   pcbddc->n_ISForDofs                = 0;
2012   pcbddc->n_ISForDofsLocal           = 0;
2013   pcbddc->ISForDofs                  = 0;
2014   pcbddc->ISForDofsLocal             = 0;
2015   pcbddc->ConstraintMatrix           = 0;
2016   pcbddc->use_exact_dirichlet_trick  = PETSC_TRUE;
2017   pcbddc->coarse_loc_to_glob         = 0;
2018   pcbddc->coarsening_ratio           = 8;
2019   pcbddc->coarse_adj_red             = 0;
2020   pcbddc->current_level              = 0;
2021   pcbddc->max_levels                 = 0;
2022   pcbddc->use_coarse_estimates       = PETSC_FALSE;
2023   pcbddc->redistribute_coarse        = 0;
2024   pcbddc->coarse_subassembling       = 0;
2025   pcbddc->coarse_subassembling_init  = 0;
2026 
2027   /* create local graph structure */
2028   ierr = PCBDDCGraphCreate(&pcbddc->mat_graph);CHKERRQ(ierr);
2029 
2030   /* scaling */
2031   pcbddc->work_scaling          = 0;
2032   pcbddc->use_deluxe_scaling    = PETSC_FALSE;
2033   pcbddc->faster_deluxe         = PETSC_FALSE;
2034 
2035   /* create sub schurs structure */
2036   ierr = PCBDDCSubSchursCreate(&pcbddc->sub_schurs);CHKERRQ(ierr);
2037   pcbddc->sub_schurs_rebuild     = PETSC_FALSE;
2038   pcbddc->sub_schurs_layers      = -1;
2039   pcbddc->sub_schurs_use_useradj = PETSC_FALSE;
2040 
2041   pcbddc->computed_rowadj = PETSC_FALSE;
2042 
2043   /* adaptivity */
2044   pcbddc->adaptive_threshold      = 0.0;
2045   pcbddc->adaptive_nmax           = 0;
2046   pcbddc->adaptive_nmin           = 0;
2047 
2048   /* function pointers */
2049   pc->ops->apply               = PCApply_BDDC;
2050   pc->ops->applytranspose      = PCApplyTranspose_BDDC;
2051   pc->ops->setup               = PCSetUp_BDDC;
2052   pc->ops->destroy             = PCDestroy_BDDC;
2053   pc->ops->setfromoptions      = PCSetFromOptions_BDDC;
2054   pc->ops->view                = 0;
2055   pc->ops->applyrichardson     = 0;
2056   pc->ops->applysymmetricleft  = 0;
2057   pc->ops->applysymmetricright = 0;
2058   pc->ops->presolve            = PCPreSolve_BDDC;
2059   pc->ops->postsolve           = PCPostSolve_BDDC;
2060 
2061   /* composing function */
2062   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",PCBDDCSetChangeOfBasisMat_BDDC);CHKERRQ(ierr);
2063   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);CHKERRQ(ierr);
2064   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);CHKERRQ(ierr);
2065   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",PCBDDCSetLevel_BDDC);CHKERRQ(ierr);
2066   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",PCBDDCSetUseExactDirichlet_BDDC);CHKERRQ(ierr);
2067   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",PCBDDCSetLevels_BDDC);CHKERRQ(ierr);
2068   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",PCBDDCSetNullSpace_BDDC);CHKERRQ(ierr);
2069   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",PCBDDCSetDirichletBoundaries_BDDC);CHKERRQ(ierr);
2070   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",PCBDDCSetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr);
2071   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",PCBDDCSetNeumannBoundaries_BDDC);CHKERRQ(ierr);
2072   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",PCBDDCSetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr);
2073   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",PCBDDCGetDirichletBoundaries_BDDC);CHKERRQ(ierr);
2074   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",PCBDDCGetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr);
2075   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",PCBDDCGetNeumannBoundaries_BDDC);CHKERRQ(ierr);
2076   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",PCBDDCGetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr);
2077   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);CHKERRQ(ierr);
2078   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",PCBDDCSetDofsSplittingLocal_BDDC);CHKERRQ(ierr);
2079   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);CHKERRQ(ierr);
2080   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);CHKERRQ(ierr);
2081   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);CHKERRQ(ierr);
2082   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);CHKERRQ(ierr);
2083   PetscFunctionReturn(0);
2084 }
2085 
2086