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