xref: /petsc/src/ksp/pc/impls/bddc/bddcprivate.c (revision aed5ffcbc31e41b96407a92d0c19d6464f10f960) 
1 #include "bddc.h"
2 #include "bddcprivate.h"
3 #include <petscblaslapack.h>
4 
5 #undef __FUNCT__
6 #define __FUNCT__ "PCBDDCResetCustomization"
7 PetscErrorCode PCBDDCResetCustomization(PC pc)
8 {
9   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
10   PetscInt       i;
11   PetscErrorCode ierr;
12 
13   PetscFunctionBegin;
14   ierr = PCBDDCGraphResetCSR(pcbddc->mat_graph);CHKERRQ(ierr);
15   ierr = ISDestroy(&pcbddc->user_primal_vertices);CHKERRQ(ierr);
16   ierr = MatNullSpaceDestroy(&pcbddc->NullSpace);CHKERRQ(ierr);
17   ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr);
18   ierr = ISDestroy(&pcbddc->DirichletBoundaries);CHKERRQ(ierr);
19   for (i=0;i<pcbddc->n_ISForDofs;i++) {
20     ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr);
21   }
22   ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr);
23   PetscFunctionReturn(0);
24 }
25 
26 #undef __FUNCT__
27 #define __FUNCT__ "PCBDDCResetTopography"
28 PetscErrorCode PCBDDCResetTopography(PC pc)
29 {
30   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
31   PetscErrorCode ierr;
32 
33   PetscFunctionBegin;
34   ierr = MatDestroy(&pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr);
35   ierr = MatDestroy(&pcbddc->ConstraintMatrix);CHKERRQ(ierr);
36   ierr = PCBDDCGraphReset(pcbddc->mat_graph);CHKERRQ(ierr);
37   PetscFunctionReturn(0);
38 }
39 
40 #undef __FUNCT__
41 #define __FUNCT__ "PCBDDCResetSolvers"
42 PetscErrorCode PCBDDCResetSolvers(PC pc)
43 {
44   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
45   PetscErrorCode ierr;
46 
47   PetscFunctionBegin;
48   ierr = VecDestroy(&pcbddc->temp_solution);CHKERRQ(ierr);
49   ierr = VecDestroy(&pcbddc->original_rhs);CHKERRQ(ierr);
50   ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr);
51   ierr = VecDestroy(&pcbddc->coarse_vec);CHKERRQ(ierr);
52   ierr = VecDestroy(&pcbddc->coarse_rhs);CHKERRQ(ierr);
53   ierr = KSPDestroy(&pcbddc->coarse_ksp);CHKERRQ(ierr);
54   ierr = MatDestroy(&pcbddc->coarse_mat);CHKERRQ(ierr);
55   ierr = MatDestroy(&pcbddc->coarse_phi_B);CHKERRQ(ierr);
56   ierr = MatDestroy(&pcbddc->coarse_phi_D);CHKERRQ(ierr);
57   ierr = MatDestroy(&pcbddc->coarse_psi_B);CHKERRQ(ierr);
58   ierr = MatDestroy(&pcbddc->coarse_psi_D);CHKERRQ(ierr);
59   ierr = VecDestroy(&pcbddc->vec1_P);CHKERRQ(ierr);
60   ierr = VecDestroy(&pcbddc->vec1_C);CHKERRQ(ierr);
61   ierr = MatDestroy(&pcbddc->local_auxmat1);CHKERRQ(ierr);
62   ierr = MatDestroy(&pcbddc->local_auxmat2);CHKERRQ(ierr);
63   ierr = VecDestroy(&pcbddc->vec1_R);CHKERRQ(ierr);
64   ierr = VecDestroy(&pcbddc->vec2_R);CHKERRQ(ierr);
65   ierr = VecDestroy(&pcbddc->vec4_D);CHKERRQ(ierr);
66   ierr = VecScatterDestroy(&pcbddc->R_to_B);CHKERRQ(ierr);
67   ierr = VecScatterDestroy(&pcbddc->R_to_D);CHKERRQ(ierr);
68   ierr = VecScatterDestroy(&pcbddc->coarse_loc_to_glob);CHKERRQ(ierr);
69   ierr = KSPDestroy(&pcbddc->ksp_D);CHKERRQ(ierr);
70   ierr = KSPDestroy(&pcbddc->ksp_R);CHKERRQ(ierr);
71   ierr = PetscFree(pcbddc->local_primal_indices);CHKERRQ(ierr);
72   ierr = PetscFree(pcbddc->replicated_local_primal_indices);CHKERRQ(ierr);
73   ierr = PetscFree(pcbddc->replicated_local_primal_values);CHKERRQ(ierr);
74   ierr = PetscFree(pcbddc->local_primal_displacements);CHKERRQ(ierr);
75   ierr = PetscFree(pcbddc->local_primal_sizes);CHKERRQ(ierr);
76   PetscFunctionReturn(0);
77 }
78 
79 #undef __FUNCT__
80 #define __FUNCT__ "PCBDDCSolveSaddlePoint"
81 static PetscErrorCode  PCBDDCSolveSaddlePoint(PC pc)
82 {
83   PetscErrorCode ierr;
84   PC_BDDC*       pcbddc = (PC_BDDC*)(pc->data);
85 
86   PetscFunctionBegin;
87   ierr = KSPSolve(pcbddc->ksp_R,pcbddc->vec1_R,pcbddc->vec2_R);CHKERRQ(ierr);
88   if (pcbddc->local_auxmat1) {
89     ierr = MatMult(pcbddc->local_auxmat1,pcbddc->vec2_R,pcbddc->vec1_C);CHKERRQ(ierr);
90     ierr = MatMultAdd(pcbddc->local_auxmat2,pcbddc->vec1_C,pcbddc->vec2_R,pcbddc->vec2_R);CHKERRQ(ierr);
91   }
92   PetscFunctionReturn(0);
93 }
94 
95 #undef __FUNCT__
96 #define __FUNCT__ "PCBDDCApplyInterfacePreconditioner"
97 PetscErrorCode  PCBDDCApplyInterfacePreconditioner(PC pc)
98 {
99   PetscErrorCode ierr;
100   PC_BDDC*        pcbddc = (PC_BDDC*)(pc->data);
101   PC_IS*            pcis = (PC_IS*)  (pc->data);
102   const PetscScalar zero = 0.0;
103 
104   PetscFunctionBegin;
105   /* Application of PHI^T (or PSI^T)  */
106   if (pcbddc->coarse_psi_B) {
107     ierr = MatMultTranspose(pcbddc->coarse_psi_B,pcis->vec1_B,pcbddc->vec1_P);CHKERRQ(ierr);
108     if (pcbddc->inexact_prec_type) { ierr = MatMultTransposeAdd(pcbddc->coarse_psi_D,pcis->vec1_D,pcbddc->vec1_P,pcbddc->vec1_P);CHKERRQ(ierr); }
109   } else {
110     ierr = MatMultTranspose(pcbddc->coarse_phi_B,pcis->vec1_B,pcbddc->vec1_P);CHKERRQ(ierr);
111     if (pcbddc->inexact_prec_type) { ierr = MatMultTransposeAdd(pcbddc->coarse_phi_D,pcis->vec1_D,pcbddc->vec1_P,pcbddc->vec1_P);CHKERRQ(ierr); }
112   }
113   /* Scatter data of coarse_rhs */
114   if (pcbddc->coarse_rhs) { ierr = VecSet(pcbddc->coarse_rhs,zero);CHKERRQ(ierr); }
115   ierr = PCBDDCScatterCoarseDataBegin(pc,pcbddc->vec1_P,pcbddc->coarse_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
116 
117   /* Local solution on R nodes */
118   ierr = VecSet(pcbddc->vec1_R,zero);CHKERRQ(ierr);
119   ierr = VecScatterBegin(pcbddc->R_to_B,pcis->vec1_B,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
120   ierr = VecScatterEnd  (pcbddc->R_to_B,pcis->vec1_B,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
121   if (pcbddc->inexact_prec_type) {
122     ierr = VecScatterBegin(pcbddc->R_to_D,pcis->vec1_D,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
123     ierr = VecScatterEnd  (pcbddc->R_to_D,pcis->vec1_D,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
124   }
125   ierr = PCBDDCSolveSaddlePoint(pc);CHKERRQ(ierr);
126   ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr);
127   ierr = VecScatterBegin(pcbddc->R_to_B,pcbddc->vec2_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
128   ierr = VecScatterEnd  (pcbddc->R_to_B,pcbddc->vec2_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
129   if (pcbddc->inexact_prec_type) {
130     ierr = VecScatterBegin(pcbddc->R_to_D,pcbddc->vec2_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
131     ierr = VecScatterEnd  (pcbddc->R_to_D,pcbddc->vec2_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
132   }
133 
134   /* Coarse solution */
135   ierr = PCBDDCScatterCoarseDataEnd(pc,pcbddc->vec1_P,pcbddc->coarse_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
136   if (pcbddc->coarse_rhs) { /* TODO remove null space when doing multilevel */
137     ierr = KSPSolve(pcbddc->coarse_ksp,pcbddc->coarse_rhs,pcbddc->coarse_vec);CHKERRQ(ierr);
138   }
139   ierr = PCBDDCScatterCoarseDataBegin(pc,pcbddc->coarse_vec,pcbddc->vec1_P,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
140   ierr = PCBDDCScatterCoarseDataEnd  (pc,pcbddc->coarse_vec,pcbddc->vec1_P,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
141 
142   /* Sum contributions from two levels */
143   ierr = MatMultAdd(pcbddc->coarse_phi_B,pcbddc->vec1_P,pcis->vec1_B,pcis->vec1_B);CHKERRQ(ierr);
144   if (pcbddc->inexact_prec_type) { ierr = MatMultAdd(pcbddc->coarse_phi_D,pcbddc->vec1_P,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); }
145   PetscFunctionReturn(0);
146 }
147 
148 #undef __FUNCT__
149 #define __FUNCT__ "PCBDDCScatterCoarseDataBegin"
150 PetscErrorCode PCBDDCScatterCoarseDataBegin(PC pc,Vec vec_from, Vec vec_to, InsertMode imode, ScatterMode smode)
151 {
152   PetscErrorCode ierr;
153   PC_BDDC*       pcbddc = (PC_BDDC*)(pc->data);
154 
155   PetscFunctionBegin;
156   switch (pcbddc->coarse_communications_type) {
157     case SCATTERS_BDDC:
158       ierr = VecScatterBegin(pcbddc->coarse_loc_to_glob,vec_from,vec_to,imode,smode);CHKERRQ(ierr);
159       break;
160     case GATHERS_BDDC:
161       break;
162   }
163   PetscFunctionReturn(0);
164 }
165 
166 #undef __FUNCT__
167 #define __FUNCT__ "PCBDDCScatterCoarseDataEnd"
168 PetscErrorCode PCBDDCScatterCoarseDataEnd(PC pc,Vec vec_from, Vec vec_to, InsertMode imode, ScatterMode smode)
169 {
170   PetscErrorCode ierr;
171   PC_BDDC*       pcbddc = (PC_BDDC*)(pc->data);
172   PetscScalar*   array_to;
173   PetscScalar*   array_from;
174   MPI_Comm       comm;
175   PetscInt       i;
176 
177   PetscFunctionBegin;
178   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
179   switch (pcbddc->coarse_communications_type) {
180     case SCATTERS_BDDC:
181       ierr = VecScatterEnd(pcbddc->coarse_loc_to_glob,vec_from,vec_to,imode,smode);CHKERRQ(ierr);
182       break;
183     case GATHERS_BDDC:
184       if (vec_from) {
185         ierr = VecGetArray(vec_from,&array_from);CHKERRQ(ierr);
186       }
187       if (vec_to) {
188         ierr = VecGetArray(vec_to,&array_to);CHKERRQ(ierr);
189       }
190       switch(pcbddc->coarse_problem_type){
191         case SEQUENTIAL_BDDC:
192           if (smode == SCATTER_FORWARD) {
193             ierr = MPI_Gatherv(&array_from[0],pcbddc->local_primal_size,MPIU_SCALAR,&pcbddc->replicated_local_primal_values[0],pcbddc->local_primal_sizes,pcbddc->local_primal_displacements,MPIU_SCALAR,0,comm);CHKERRQ(ierr);
194             if (vec_to) {
195               if (imode == ADD_VALUES) {
196                 for (i=0;i<pcbddc->replicated_primal_size;i++) {
197                   array_to[pcbddc->replicated_local_primal_indices[i]]+=pcbddc->replicated_local_primal_values[i];
198                 }
199               } else {
200                 for (i=0;i<pcbddc->replicated_primal_size;i++) {
201                   array_to[pcbddc->replicated_local_primal_indices[i]]=pcbddc->replicated_local_primal_values[i];
202                 }
203               }
204             }
205           } else {
206             if (vec_from) {
207               if (imode == ADD_VALUES) {
208                 MPI_Comm vec_from_comm;
209                 ierr = PetscObjectGetComm((PetscObject)(vec_from),&vec_from_comm);CHKERRQ(ierr);
210                 SETERRQ2(vec_from_comm,PETSC_ERR_SUP,"Unsupported insert mode ADD_VALUES for SCATTER_REVERSE in %s for case %d\n",__FUNCT__,pcbddc->coarse_problem_type);
211               }
212               for (i=0;i<pcbddc->replicated_primal_size;i++) {
213                 pcbddc->replicated_local_primal_values[i]=array_from[pcbddc->replicated_local_primal_indices[i]];
214               }
215             }
216             ierr = MPI_Scatterv(&pcbddc->replicated_local_primal_values[0],pcbddc->local_primal_sizes,pcbddc->local_primal_displacements,MPIU_SCALAR,&array_to[0],pcbddc->local_primal_size,MPIU_SCALAR,0,comm);CHKERRQ(ierr);
217           }
218           break;
219         case REPLICATED_BDDC:
220           if (smode == SCATTER_FORWARD) {
221             ierr = MPI_Allgatherv(&array_from[0],pcbddc->local_primal_size,MPIU_SCALAR,&pcbddc->replicated_local_primal_values[0],pcbddc->local_primal_sizes,pcbddc->local_primal_displacements,MPIU_SCALAR,comm);CHKERRQ(ierr);
222             if (imode == ADD_VALUES) {
223               for (i=0;i<pcbddc->replicated_primal_size;i++) {
224                 array_to[pcbddc->replicated_local_primal_indices[i]]+=pcbddc->replicated_local_primal_values[i];
225               }
226             } else {
227               for (i=0;i<pcbddc->replicated_primal_size;i++) {
228                 array_to[pcbddc->replicated_local_primal_indices[i]]=pcbddc->replicated_local_primal_values[i];
229               }
230             }
231           } else { /* no communications needed for SCATTER_REVERSE since needed data is already present */
232             if (imode == ADD_VALUES) {
233               for (i=0;i<pcbddc->local_primal_size;i++) {
234                 array_to[i]+=array_from[pcbddc->local_primal_indices[i]];
235               }
236             } else {
237               for (i=0;i<pcbddc->local_primal_size;i++) {
238                 array_to[i]=array_from[pcbddc->local_primal_indices[i]];
239               }
240             }
241           }
242           break;
243         case MULTILEVEL_BDDC:
244           break;
245         case PARALLEL_BDDC:
246           break;
247       }
248       if (vec_from) {
249         ierr = VecRestoreArray(vec_from,&array_from);CHKERRQ(ierr);
250       }
251       if (vec_to) {
252         ierr = VecRestoreArray(vec_to,&array_to);CHKERRQ(ierr);
253       }
254       break;
255   }
256   PetscFunctionReturn(0);
257 }
258 
259 #undef __FUNCT__
260 #define __FUNCT__ "PCBDDCConstraintsSetUp"
261 PetscErrorCode PCBDDCConstraintsSetUp(PC pc)
262 {
263   PetscErrorCode ierr;
264   PC_IS*         pcis = (PC_IS*)(pc->data);
265   PC_BDDC*       pcbddc = (PC_BDDC*)pc->data;
266   Mat_IS         *matis = (Mat_IS*)pc->pmat->data;
267   PetscInt       *nnz,*is_indices;
268   PetscScalar    *temp_quadrature_constraint;
269   PetscInt       *temp_indices,*temp_indices_to_constraint,*temp_indices_to_constraint_B,*local_to_B;
270   PetscInt       local_primal_size,i,j,k,total_counts,max_size_of_constraint;
271   PetscInt       n_vertices,size_of_constraint;
272   PetscReal      real_value;
273   PetscBool      nnsp_has_cnst=PETSC_FALSE,use_nnsp_true=pcbddc->use_nnsp_true;
274   PetscInt       nnsp_size=0,nnsp_addone=0,temp_constraints,temp_start_ptr,n_ISForFaces,n_ISForEdges;
275   IS             *used_IS,ISForVertices,*ISForFaces,*ISForEdges;
276   MatType        impMatType=MATSEQAIJ;
277   PetscBLASInt   Bs,Bt,lwork,lierr;
278   PetscReal      tol=1.0e-8;
279   MatNullSpace   nearnullsp;
280   const Vec      *nearnullvecs;
281   Vec            *localnearnullsp;
282   PetscScalar    *work,*temp_basis,*array_vector,*correlation_mat;
283   PetscReal      *rwork,*singular_vals;
284   PetscBLASInt   Bone=1,*ipiv;
285   Vec            temp_vec;
286   Mat            temp_mat;
287   KSP            temp_ksp;
288   PC             temp_pc;
289   PetscInt       s,start_constraint,dual_dofs;
290   PetscBool      compute_submatrix,useksp=PETSC_FALSE;
291   PetscInt       *aux_primal_permutation,*aux_primal_numbering;
292   PetscBool      boolforchange,*change_basis;
293 /* some ugly conditional declarations */
294 #if defined(PETSC_MISSING_LAPACK_GESVD)
295   PetscScalar    one=1.0,zero=0.0;
296   PetscInt       ii;
297   PetscScalar    *singular_vectors;
298   PetscBLASInt   *iwork,*ifail;
299   PetscReal      dummy_real,abs_tol;
300   PetscBLASInt   eigs_found;
301 #endif
302   PetscBLASInt   dummy_int;
303   PetscScalar    dummy_scalar;
304   PetscBool      used_vertex,get_faces,get_edges,get_vertices;
305 
306   PetscFunctionBegin;
307   /* Get index sets for faces, edges and vertices from graph */
308   get_faces = PETSC_TRUE;
309   get_edges = PETSC_TRUE;
310   get_vertices = PETSC_TRUE;
311   if (pcbddc->vertices_flag) {
312     get_faces = PETSC_FALSE;
313     get_edges = PETSC_FALSE;
314   }
315   if (pcbddc->constraints_flag) {
316     get_vertices = PETSC_FALSE;
317   }
318   if (pcbddc->faces_flag) {
319     get_edges = PETSC_FALSE;
320   }
321   if (pcbddc->edges_flag) {
322     get_faces = PETSC_FALSE;
323   }
324   /* default */
325   if (!get_faces && !get_edges && !get_vertices) {
326     get_vertices = PETSC_TRUE;
327   }
328   ierr = PCBDDCGraphGetCandidatesIS(pcbddc->mat_graph,get_faces,get_edges,get_vertices,&n_ISForFaces,&ISForFaces,&n_ISForEdges,&ISForEdges,&ISForVertices);
329   if (pcbddc->dbg_flag) {
330     ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"--------------------------------------------------------------\n");CHKERRQ(ierr);
331     i = 0;
332     if (ISForVertices) {
333       ierr = ISGetSize(ISForVertices,&i);CHKERRQ(ierr);
334     }
335     ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d got %02d local candidate vertices\n",PetscGlobalRank,i);CHKERRQ(ierr);
336     ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d got %02d local candidate edges\n",PetscGlobalRank,n_ISForEdges);CHKERRQ(ierr);
337     ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d got %02d local candidate faces\n",PetscGlobalRank,n_ISForFaces);CHKERRQ(ierr);
338     ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr);
339   }
340   /* check if near null space is attached to global mat */
341   ierr = MatGetNearNullSpace(pc->pmat,&nearnullsp);CHKERRQ(ierr);
342   if (nearnullsp) {
343     ierr = MatNullSpaceGetVecs(nearnullsp,&nnsp_has_cnst,&nnsp_size,&nearnullvecs);CHKERRQ(ierr);
344   } else { /* if near null space is not provided it uses constants */
345     nnsp_has_cnst = PETSC_TRUE;
346     use_nnsp_true = PETSC_TRUE;
347   }
348   if (nnsp_has_cnst) {
349     nnsp_addone = 1;
350   }
351   /*
352        Evaluate maximum storage size needed by the procedure
353        - temp_indices will contain start index of each constraint stored as follows
354        - temp_indices_to_constraint  [temp_indices[i],...,temp[indices[i+1]-1] will contain the indices (in local numbering) on which the constraint acts
355        - temp_indices_to_constraint_B[temp_indices[i],...,temp[indices[i+1]-1] will contain the indices (in boundary numbering) on which the constraint acts
356        - temp_quadrature_constraint  [temp_indices[i],...,temp[indices[i+1]-1] will contain the scalars representing the constraint itself
357                                                                                                                                                          */
358   total_counts = n_ISForFaces+n_ISForEdges;
359   total_counts *= (nnsp_addone+nnsp_size);
360   n_vertices = 0;
361   if (ISForVertices) {
362     ierr = ISGetSize(ISForVertices,&n_vertices);CHKERRQ(ierr);
363   }
364   total_counts += n_vertices;
365   ierr = PetscMalloc((total_counts+1)*sizeof(PetscInt),&temp_indices);CHKERRQ(ierr);
366   ierr = PetscMalloc((total_counts+1)*sizeof(PetscBool),&change_basis);CHKERRQ(ierr);
367   total_counts = 0;
368   max_size_of_constraint = 0;
369   for (i=0;i<n_ISForEdges+n_ISForFaces;i++) {
370     if (i<n_ISForEdges) {
371       used_IS = &ISForEdges[i];
372     } else {
373       used_IS = &ISForFaces[i-n_ISForEdges];
374     }
375     ierr = ISGetSize(*used_IS,&j);CHKERRQ(ierr);
376     total_counts += j;
377     max_size_of_constraint = PetscMax(j,max_size_of_constraint);
378   }
379   total_counts *= (nnsp_addone+nnsp_size);
380   total_counts += n_vertices;
381   ierr = PetscMalloc(total_counts*sizeof(PetscScalar),&temp_quadrature_constraint);CHKERRQ(ierr);
382   ierr = PetscMalloc(total_counts*sizeof(PetscInt),&temp_indices_to_constraint);CHKERRQ(ierr);
383   ierr = PetscMalloc(total_counts*sizeof(PetscInt),&temp_indices_to_constraint_B);CHKERRQ(ierr);
384   ierr = PetscMalloc(pcis->n*sizeof(PetscInt),&local_to_B);CHKERRQ(ierr);
385   ierr = ISGetIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr);
386   for (i=0;i<pcis->n;i++) {
387     local_to_B[i]=-1;
388   }
389   for (i=0;i<pcis->n_B;i++) {
390     local_to_B[is_indices[i]]=i;
391   }
392   ierr = ISRestoreIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr);
393 
394   /* First we issue queries to allocate optimal workspace for LAPACKgesvd or LAPACKsyev/LAPACKheev */
395   rwork = 0;
396   work = 0;
397   singular_vals = 0;
398   temp_basis = 0;
399   correlation_mat = 0;
400   if (!pcbddc->use_nnsp_true) {
401     PetscScalar temp_work;
402 #if defined(PETSC_MISSING_LAPACK_GESVD)
403     /* POD */
404     PetscInt max_n;
405     max_n = nnsp_addone+nnsp_size;
406     /* using some techniques borrowed from Proper Orthogonal Decomposition */
407     ierr = PetscMalloc(max_n*max_n*sizeof(PetscScalar),&correlation_mat);CHKERRQ(ierr);
408     ierr = PetscMalloc(max_n*max_n*sizeof(PetscScalar),&singular_vectors);CHKERRQ(ierr);
409     ierr = PetscMalloc(max_n*sizeof(PetscReal),&singular_vals);CHKERRQ(ierr);
410     ierr = PetscMalloc(max_size_of_constraint*(nnsp_addone+nnsp_size)*sizeof(PetscScalar),&temp_basis);CHKERRQ(ierr);
411 #if defined(PETSC_USE_COMPLEX)
412     ierr = PetscMalloc(3*max_n*sizeof(PetscReal),&rwork);CHKERRQ(ierr);
413 #endif
414     ierr = PetscMalloc(5*max_n*sizeof(PetscBLASInt),&iwork);CHKERRQ(ierr);
415     ierr = PetscMalloc(max_n*sizeof(PetscBLASInt),&ifail);CHKERRQ(ierr);
416     /* now we evaluate the optimal workspace using query with lwork=-1 */
417     ierr = PetscBLASIntCast(max_n,&Bt);CHKERRQ(ierr);
418     lwork=-1;
419     ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
420 #if !defined(PETSC_USE_COMPLEX)
421     abs_tol=1.e-8;
422 /*    LAPACKsyev_("V","U",&Bt,correlation_mat,&Bt,singular_vals,&temp_work,&lwork,&lierr); */
423     PetscStackCallBLAS("LAPACKsyevx",LAPACKsyevx_("V","A","U",&Bt,correlation_mat,&Bt,&dummy_real,&dummy_real,&dummy_int,&dummy_int,&abs_tol,&eigs_found,singular_vals,singular_vectors,&Bt,&temp_work,&lwork,iwork,ifail,&lierr));
424 #else
425 /*    LAPACKsyev_("V","U",&Bt,correlation_mat,&Bt,singular_vals,&temp_work,&lwork,rwork,&lierr); */
426 /*  LAPACK call is missing here! TODO */
427     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not yet implemented for complexes when PETSC_MISSING_GESVD = 1");
428 #endif
429     if ( lierr ) {
430       SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in query to SYEVX Lapack routine %d",(int)lierr);
431     }
432     ierr = PetscFPTrapPop();CHKERRQ(ierr);
433 #else /* on missing GESVD */
434     /* SVD */
435     PetscInt max_n,min_n;
436     max_n = max_size_of_constraint;
437     min_n = nnsp_addone+nnsp_size;
438     if (max_size_of_constraint < ( nnsp_addone+nnsp_size ) ) {
439       min_n = max_size_of_constraint;
440       max_n = nnsp_addone+nnsp_size;
441     }
442     ierr = PetscMalloc(min_n*sizeof(PetscReal),&singular_vals);CHKERRQ(ierr);
443 #if defined(PETSC_USE_COMPLEX)
444     ierr = PetscMalloc(5*min_n*sizeof(PetscReal),&rwork);CHKERRQ(ierr);
445 #endif
446     /* now we evaluate the optimal workspace using query with lwork=-1 */
447     lwork=-1;
448     ierr = PetscBLASIntCast(max_n,&Bs);CHKERRQ(ierr);
449     ierr = PetscBLASIntCast(min_n,&Bt);CHKERRQ(ierr);
450     dummy_int = Bs;
451     ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
452 #if !defined(PETSC_USE_COMPLEX)
453     PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&Bs,&Bt,&temp_quadrature_constraint[0],&Bs,singular_vals,&dummy_scalar,&dummy_int,&dummy_scalar,&dummy_int,&temp_work,&lwork,&lierr));
454 #else
455     PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&Bs,&Bt,&temp_quadrature_constraint[0],&Bs,singular_vals,&dummy_scalar,&dummy_int,&dummy_scalar,&dummy_int,&temp_work,&lwork,rwork,&lierr));
456 #endif
457     if ( lierr ) {
458       SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in query to SVD Lapack routine %d",(int)lierr);
459     }
460     ierr = PetscFPTrapPop();CHKERRQ(ierr);
461 #endif
462     /* Allocate optimal workspace */
463     ierr = PetscBLASIntCast((PetscInt)PetscRealPart(temp_work),&lwork);CHKERRQ(ierr);
464     total_counts = (PetscInt)lwork;
465     ierr = PetscMalloc(total_counts*sizeof(PetscScalar),&work);CHKERRQ(ierr);
466   }
467   /* get local part of global near null space vectors */
468   ierr = PetscMalloc(nnsp_size*sizeof(Vec),&localnearnullsp);CHKERRQ(ierr);
469   for (k=0;k<nnsp_size;k++) {
470     ierr = VecDuplicate(pcis->vec1_N,&localnearnullsp[k]);CHKERRQ(ierr);
471     ierr = VecScatterBegin(matis->ctx,nearnullvecs[k],localnearnullsp[k],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
472     ierr = VecScatterEnd(matis->ctx,nearnullvecs[k],localnearnullsp[k],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
473   }
474   /* Now we can loop on constraining sets */
475   total_counts = 0;
476   temp_indices[0] = 0;
477   /* vertices */
478   if (ISForVertices) {
479     ierr = ISGetIndices(ISForVertices,(const PetscInt**)&is_indices);CHKERRQ(ierr);
480     if (nnsp_has_cnst) { /* consider all vertices */
481       for (i=0;i<n_vertices;i++) {
482         temp_indices_to_constraint[temp_indices[total_counts]]=is_indices[i];
483         temp_indices_to_constraint_B[temp_indices[total_counts]]=local_to_B[is_indices[i]];
484         temp_quadrature_constraint[temp_indices[total_counts]]=1.0;
485         temp_indices[total_counts+1]=temp_indices[total_counts]+1;
486         change_basis[total_counts]=PETSC_FALSE;
487         total_counts++;
488       }
489     } else { /* consider vertices for which exist at least a localnearnullsp which is not null there */
490       for (i=0;i<n_vertices;i++) {
491         used_vertex=PETSC_FALSE;
492         k=0;
493         while (!used_vertex && k<nnsp_size) {
494           ierr = VecGetArrayRead(localnearnullsp[k],(const PetscScalar**)&array_vector);CHKERRQ(ierr);
495           if (PetscAbsScalar(array_vector[is_indices[i]])>0.0) {
496             temp_indices_to_constraint[temp_indices[total_counts]]=is_indices[i];
497             temp_indices_to_constraint_B[temp_indices[total_counts]]=local_to_B[is_indices[i]];
498             temp_quadrature_constraint[temp_indices[total_counts]]=1.0;
499             temp_indices[total_counts+1]=temp_indices[total_counts]+1;
500             change_basis[total_counts]=PETSC_FALSE;
501             total_counts++;
502             used_vertex=PETSC_TRUE;
503           }
504           ierr = VecRestoreArrayRead(localnearnullsp[k],(const PetscScalar**)&array_vector);CHKERRQ(ierr);
505           k++;
506         }
507       }
508     }
509     ierr = ISRestoreIndices(ISForVertices,(const PetscInt**)&is_indices);CHKERRQ(ierr);
510     n_vertices = total_counts;
511   }
512   /* edges and faces */
513   for (i=0;i<n_ISForEdges+n_ISForFaces;i++) {
514     if (i<n_ISForEdges) {
515       used_IS = &ISForEdges[i];
516       boolforchange = pcbddc->use_change_of_basis;
517     } else {
518       used_IS = &ISForFaces[i-n_ISForEdges];
519       boolforchange = pcbddc->use_change_on_faces;
520     }
521     temp_constraints = 0;          /* zero the number of constraints I have on this conn comp */
522     temp_start_ptr = total_counts; /* need to know the starting index of constraints stored */
523     ierr = ISGetSize(*used_IS,&size_of_constraint);CHKERRQ(ierr);
524     ierr = ISGetIndices(*used_IS,(const PetscInt**)&is_indices);CHKERRQ(ierr);
525     /* HACK: change of basis should not performed on local periodic nodes */
526     if (pcbddc->mat_graph->mirrors && pcbddc->mat_graph->mirrors[is_indices[0]]) {
527       boolforchange = PETSC_FALSE;
528     }
529     if (nnsp_has_cnst) {
530       PetscScalar quad_value;
531       temp_constraints++;
532       quad_value = (PetscScalar)(1.0/PetscSqrtReal((PetscReal)size_of_constraint));
533       for (j=0;j<size_of_constraint;j++) {
534         temp_indices_to_constraint[temp_indices[total_counts]+j]=is_indices[j];
535         temp_indices_to_constraint_B[temp_indices[total_counts]+j]=local_to_B[is_indices[j]];
536         temp_quadrature_constraint[temp_indices[total_counts]+j]=quad_value;
537       }
538       temp_indices[total_counts+1]=temp_indices[total_counts]+size_of_constraint;  /* store new starting point */
539       change_basis[total_counts]=boolforchange;
540       total_counts++;
541     }
542     for (k=0;k<nnsp_size;k++) {
543       ierr = VecGetArrayRead(localnearnullsp[k],(const PetscScalar**)&array_vector);CHKERRQ(ierr);
544       for (j=0;j<size_of_constraint;j++) {
545         temp_indices_to_constraint[temp_indices[total_counts]+j]=is_indices[j];
546         temp_indices_to_constraint_B[temp_indices[total_counts]+j]=local_to_B[is_indices[j]];
547         temp_quadrature_constraint[temp_indices[total_counts]+j]=array_vector[is_indices[j]];
548       }
549       ierr = VecRestoreArrayRead(localnearnullsp[k],(const PetscScalar**)&array_vector);CHKERRQ(ierr);
550       real_value = 1.0;
551       if (use_nnsp_true) { /* check if array is null on the connected component in case use_nnsp_true has been requested */
552         ierr = PetscBLASIntCast(size_of_constraint,&Bs);CHKERRQ(ierr);
553         PetscStackCallBLAS("BLASasum",real_value = BLASasum_(&Bs,&temp_quadrature_constraint[temp_indices[total_counts]],&Bone));
554       }
555       if (real_value > 0.0) { /* keep indices and values */
556         temp_constraints++;
557         temp_indices[total_counts+1]=temp_indices[total_counts]+size_of_constraint;  /* store new starting point */
558         change_basis[total_counts]=boolforchange;
559         total_counts++;
560       }
561     }
562     ierr = ISRestoreIndices(*used_IS,(const PetscInt**)&is_indices);CHKERRQ(ierr);
563     /* perform SVD on the constraint if use_nnsp_true has not be requested by the user */
564     if (!use_nnsp_true) {
565       ierr = PetscBLASIntCast(size_of_constraint,&Bs);CHKERRQ(ierr);
566       ierr = PetscBLASIntCast(temp_constraints,&Bt);CHKERRQ(ierr);
567 
568 #if defined(PETSC_MISSING_LAPACK_GESVD)
569       ierr = PetscMemzero(correlation_mat,Bt*Bt*sizeof(PetscScalar));CHKERRQ(ierr);
570       /* Store upper triangular part of correlation matrix */
571       for (j=0;j<temp_constraints;j++) {
572         for (k=0;k<j+1;k++) {
573           PetscStackCallBLAS("BLASdot",correlation_mat[j*temp_constraints+k]=BLASdot_(&Bs,&temp_quadrature_constraint[temp_indices[temp_start_ptr+j]],&Bone,&temp_quadrature_constraint[temp_indices[temp_start_ptr+k]],&Bone));
574 
575         }
576       }
577       ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
578 #if !defined(PETSC_USE_COMPLEX)
579 /*      LAPACKsyev_("V","U",&Bt,correlation_mat,&Bt,singular_vals,work,&lwork,&lierr); */
580       PetscStackCallBLAS("LAPACKsyevx",LAPACKsyevx_("V","A","U",&Bt,correlation_mat,&Bt,&dummy_real,&dummy_real,&dummy_int,&dummy_int,&abs_tol,&eigs_found,singular_vals,singular_vectors,&Bt,work,&lwork,iwork,ifail,&lierr));
581 #else
582 /*  LAPACK call is missing here! TODO */
583       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not yet implemented for complexes when PETSC_MISSING_GESVD = 1");
584 #endif
585       if (lierr) {
586         SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SYEVX Lapack routine %d",(int)lierr);
587       }
588       ierr = PetscFPTrapPop();CHKERRQ(ierr);
589       /* retain eigenvalues greater than tol: note that lapack SYEV gives eigs in ascending order */
590       j=0;
591       while (j < Bt && singular_vals[j] < tol) j++;
592       total_counts=total_counts-j;
593       if (j<temp_constraints) {
594         for (k=j;k<Bt;k++) {
595           singular_vals[k]=1.0/PetscSqrtReal(singular_vals[k]);
596         }
597         ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
598         PetscStackCallBLAS("BLASgemm_",BLASgemm_("N","N",&Bs,&Bt,&Bt,&one,&temp_quadrature_constraint[temp_indices[temp_start_ptr]],&Bs,correlation_mat,&Bt,&zero,temp_basis,&Bs));
599         ierr = PetscFPTrapPop();CHKERRQ(ierr);
600         /* copy POD basis into used quadrature memory */
601         for (k=0;k<Bt-j;k++) {
602           for (ii=0;ii<size_of_constraint;ii++) {
603             temp_quadrature_constraint[temp_indices[temp_start_ptr+k]+ii]=singular_vals[Bt-1-k]*temp_basis[(Bt-1-k)*size_of_constraint+ii];
604           }
605         }
606       }
607 
608 #else  /* on missing GESVD */
609       PetscInt min_n = temp_constraints;
610       if (min_n > size_of_constraint) min_n = size_of_constraint;
611       dummy_int = Bs;
612       ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
613 #if !defined(PETSC_USE_COMPLEX)
614       PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&Bs,&Bt,&temp_quadrature_constraint[temp_indices[temp_start_ptr]],&Bs,singular_vals,&dummy_scalar,&dummy_int,&dummy_scalar,&dummy_int,work,&lwork,&lierr));
615 #else
616       PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&Bs,&Bt,&temp_quadrature_constraint[temp_indices[temp_start_ptr]],&Bs,singular_vals,&dummy_scalar,&dummy_int,&dummy_scalar,&dummy_int,work,&lwork,rwork,&lierr));
617 #endif
618       if (lierr) {
619         SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SVD Lapack routine %d",(int)lierr);
620       }
621       ierr = PetscFPTrapPop();CHKERRQ(ierr);
622       /* retain eigenvalues greater than tol: note that lapack SVD gives eigs in descending order */
623       j = 0;
624       while (j < min_n && singular_vals[min_n-j-1] < tol) j++;
625       total_counts = total_counts-(PetscInt)Bt+(min_n-j);
626 #endif
627     }
628   }
629   /* free index sets of faces, edges and vertices */
630   for (i=0;i<n_ISForFaces;i++) {
631     ierr = ISDestroy(&ISForFaces[i]);CHKERRQ(ierr);
632   }
633   ierr = PetscFree(ISForFaces);CHKERRQ(ierr);
634   for (i=0;i<n_ISForEdges;i++) {
635     ierr = ISDestroy(&ISForEdges[i]);CHKERRQ(ierr);
636   }
637   ierr = PetscFree(ISForEdges);CHKERRQ(ierr);
638   ierr = ISDestroy(&ISForVertices);CHKERRQ(ierr);
639 
640   /* set quantities in pcbddc data structure */
641   /* n_vertices defines the number of point primal dofs */
642   /* n_constraints defines the number of averages (they can be point primal dofs if change of basis is requested) */
643   local_primal_size = total_counts;
644   pcbddc->n_vertices = n_vertices;
645   pcbddc->n_constraints = total_counts-n_vertices;
646   pcbddc->local_primal_size = local_primal_size;
647 
648   /* Create constraint matrix */
649   /* The constraint matrix is used to compute the l2g map of primal dofs */
650   /* so we need to set it up properly either with or without change of basis */
651   ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->ConstraintMatrix);CHKERRQ(ierr);
652   ierr = MatSetType(pcbddc->ConstraintMatrix,impMatType);CHKERRQ(ierr);
653   ierr = MatSetSizes(pcbddc->ConstraintMatrix,local_primal_size,pcis->n,local_primal_size,pcis->n);CHKERRQ(ierr);
654   /* compute a local numbering of constraints : vertices first then constraints */
655   ierr = VecSet(pcis->vec1_N,0.0);CHKERRQ(ierr);
656   ierr = VecGetArray(pcis->vec1_N,&array_vector);CHKERRQ(ierr);
657   ierr = PetscMalloc(local_primal_size*sizeof(PetscInt),&aux_primal_numbering);CHKERRQ(ierr);
658   ierr = PetscMalloc(local_primal_size*sizeof(PetscInt),&aux_primal_permutation);CHKERRQ(ierr);
659   total_counts=0;
660   /* find vertices: subdomain corners plus dofs with basis changed */
661   for (i=0;i<local_primal_size;i++) {
662     size_of_constraint=temp_indices[i+1]-temp_indices[i];
663     if (change_basis[i] || size_of_constraint == 1) {
664       k=0;
665       while(k < size_of_constraint && array_vector[temp_indices_to_constraint[temp_indices[i]+size_of_constraint-k-1]] != 0.0) {
666         k=k+1;
667       }
668       j=temp_indices_to_constraint[temp_indices[i]+size_of_constraint-k-1];
669       array_vector[j] = 1.0;
670       aux_primal_numbering[total_counts]=j;
671       aux_primal_permutation[total_counts]=total_counts;
672       total_counts++;
673     }
674   }
675   ierr = VecRestoreArray(pcis->vec1_N,&array_vector);CHKERRQ(ierr);
676   /* permute indices in order to have a sorted set of vertices */
677   ierr = PetscSortIntWithPermutation(total_counts,aux_primal_numbering,aux_primal_permutation);
678   /* nonzero structure */
679   ierr = PetscMalloc(local_primal_size*sizeof(PetscInt),&nnz);CHKERRQ(ierr);
680   for (i=0;i<total_counts;i++) {
681     nnz[i]=1;
682   }
683   j=total_counts;
684   for (i=n_vertices;i<local_primal_size;i++) {
685     if (!change_basis[i]) {
686       nnz[j]=temp_indices[i+1]-temp_indices[i];
687       j++;
688     }
689   }
690   ierr = MatSeqAIJSetPreallocation(pcbddc->ConstraintMatrix,0,nnz);CHKERRQ(ierr);
691   ierr = PetscFree(nnz);CHKERRQ(ierr);
692   /* set values in constraint matrix */
693   for (i=0;i<total_counts;i++) {
694     j = aux_primal_permutation[i];
695     k = aux_primal_numbering[j];
696     ierr = MatSetValue(pcbddc->ConstraintMatrix,i,k,1.0,INSERT_VALUES);CHKERRQ(ierr);
697   }
698   for (i=n_vertices;i<local_primal_size;i++) {
699     if (!change_basis[i]) {
700       size_of_constraint=temp_indices[i+1]-temp_indices[i];
701       ierr = MatSetValues(pcbddc->ConstraintMatrix,1,&total_counts,size_of_constraint,&temp_indices_to_constraint[temp_indices[i]],&temp_quadrature_constraint[temp_indices[i]],INSERT_VALUES);CHKERRQ(ierr);
702       total_counts++;
703     }
704   }
705   ierr = PetscFree(aux_primal_numbering);CHKERRQ(ierr);
706   ierr = PetscFree(aux_primal_permutation);CHKERRQ(ierr);
707   /* assembling */
708   ierr = MatAssemblyBegin(pcbddc->ConstraintMatrix,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
709   ierr = MatAssemblyEnd(pcbddc->ConstraintMatrix,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
710 
711   /* Create matrix for change of basis. We don't need it in case pcbddc->use_change_of_basis is FALSE */
712   if (pcbddc->use_change_of_basis) {
713     ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr);
714     ierr = MatSetType(pcbddc->ChangeOfBasisMatrix,impMatType);CHKERRQ(ierr);
715     ierr = MatSetSizes(pcbddc->ChangeOfBasisMatrix,pcis->n_B,pcis->n_B,pcis->n_B,pcis->n_B);CHKERRQ(ierr);
716     /* work arrays */
717     /* we need to reuse these arrays, so we free them */
718     ierr = PetscFree(temp_basis);CHKERRQ(ierr);
719     ierr = PetscFree(work);CHKERRQ(ierr);
720     ierr = PetscMalloc(pcis->n_B*sizeof(PetscInt),&nnz);CHKERRQ(ierr);
721     ierr = PetscMalloc((nnsp_addone+nnsp_size)*(nnsp_addone+nnsp_size)*sizeof(PetscScalar),&temp_basis);CHKERRQ(ierr);
722     ierr = PetscMalloc((nnsp_addone+nnsp_size)*sizeof(PetscScalar),&work);CHKERRQ(ierr);
723     ierr = PetscMalloc((nnsp_addone+nnsp_size)*sizeof(PetscBLASInt),&ipiv);CHKERRQ(ierr);
724     for (i=0;i<pcis->n_B;i++) {
725       nnz[i]=1;
726     }
727     /* Overestimated nonzeros per row */
728     k=1;
729     for (i=pcbddc->n_vertices;i<local_primal_size;i++) {
730       if (change_basis[i]) {
731         size_of_constraint = temp_indices[i+1]-temp_indices[i];
732         if (k < size_of_constraint) {
733           k = size_of_constraint;
734         }
735         for (j=0;j<size_of_constraint;j++) {
736           nnz[temp_indices_to_constraint_B[temp_indices[i]+j]] = size_of_constraint;
737         }
738       }
739     }
740     ierr = MatSeqAIJSetPreallocation(pcbddc->ChangeOfBasisMatrix,0,nnz);CHKERRQ(ierr);
741     ierr = PetscFree(nnz);CHKERRQ(ierr);
742     /* Temporary array to store indices */
743     ierr = PetscMalloc(k*sizeof(PetscInt),&is_indices);CHKERRQ(ierr);
744     /* Set initial identity in the matrix */
745     for (i=0;i<pcis->n_B;i++) {
746       ierr = MatSetValue(pcbddc->ChangeOfBasisMatrix,i,i,1.0,INSERT_VALUES);CHKERRQ(ierr);
747     }
748     /* Now we loop on the constraints which need a change of basis */
749     /* Change of basis matrix is evaluated as the FIRST APPROACH in */
750     /* Klawonn and Widlund, Dual-primal FETI-DP methods for linear elasticity, (6.2.1) */
751     temp_constraints = 0;
752     if (pcbddc->n_vertices < local_primal_size) {
753       temp_start_ptr = temp_indices_to_constraint_B[temp_indices[pcbddc->n_vertices]];
754     }
755     for (i=pcbddc->n_vertices;i<local_primal_size;i++) {
756       if (change_basis[i]) {
757         compute_submatrix = PETSC_FALSE;
758         useksp = PETSC_FALSE;
759         if (temp_start_ptr == temp_indices_to_constraint_B[temp_indices[i]]) {
760           temp_constraints++;
761           if (i == local_primal_size -1 ||  temp_start_ptr != temp_indices_to_constraint_B[temp_indices[i+1]]) {
762             compute_submatrix = PETSC_TRUE;
763           }
764         }
765         if (compute_submatrix) {
766           if (temp_constraints > 1 || pcbddc->use_nnsp_true) {
767             useksp = PETSC_TRUE;
768           }
769           size_of_constraint = temp_indices[i+1]-temp_indices[i];
770           if (useksp) { /* experimental TODO: reuse KSP and MAT instead of creating them each time */
771             ierr = MatCreate(PETSC_COMM_SELF,&temp_mat);CHKERRQ(ierr);
772             ierr = MatSetType(temp_mat,impMatType);CHKERRQ(ierr);
773             ierr = MatSetSizes(temp_mat,size_of_constraint,size_of_constraint,size_of_constraint,size_of_constraint);CHKERRQ(ierr);
774             ierr = MatSeqAIJSetPreallocation(temp_mat,size_of_constraint,NULL);CHKERRQ(ierr);
775           }
776           /* First _size_of_constraint-temp_constraints_ columns */
777           dual_dofs = size_of_constraint-temp_constraints;
778           start_constraint = i+1-temp_constraints;
779           for (s=0;s<dual_dofs;s++) {
780             is_indices[0] = s;
781             for (j=0;j<temp_constraints;j++) {
782               for (k=0;k<temp_constraints;k++) {
783                 temp_basis[j*temp_constraints+k]=temp_quadrature_constraint[temp_indices[start_constraint+k]+s+j+1];
784               }
785               work[j]=-temp_quadrature_constraint[temp_indices[start_constraint+j]+s];
786               is_indices[j+1]=s+j+1;
787             }
788             Bt = temp_constraints;
789             ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
790             PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&Bt,&Bone,temp_basis,&Bt,ipiv,work,&Bt,&lierr));
791             if ( lierr ) {
792               SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GESV Lapack routine %d",(int)lierr);
793             }
794             ierr = PetscFPTrapPop();CHKERRQ(ierr);
795             j = temp_indices_to_constraint_B[temp_indices[start_constraint]+s];
796             ierr = MatSetValues(pcbddc->ChangeOfBasisMatrix,temp_constraints,&temp_indices_to_constraint_B[temp_indices[start_constraint]+s+1],1,&j,work,INSERT_VALUES);CHKERRQ(ierr);
797             if (useksp) {
798               /* temp mat with transposed rows and columns */
799               ierr = MatSetValues(temp_mat,1,&s,temp_constraints,&is_indices[1],work,INSERT_VALUES);CHKERRQ(ierr);
800               ierr = MatSetValue(temp_mat,is_indices[0],is_indices[0],1.0,INSERT_VALUES);CHKERRQ(ierr);
801             }
802           }
803           if (useksp) {
804             /* last rows of temp_mat */
805             for (j=0;j<size_of_constraint;j++) {
806               is_indices[j] = j;
807             }
808             for (s=0;s<temp_constraints;s++) {
809               k = s + dual_dofs;
810               ierr = MatSetValues(temp_mat,1,&k,size_of_constraint,is_indices,&temp_quadrature_constraint[temp_indices[start_constraint+s]],INSERT_VALUES);CHKERRQ(ierr);
811             }
812             ierr = MatAssemblyBegin(temp_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
813             ierr = MatAssemblyEnd(temp_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
814             ierr = MatGetVecs(temp_mat,&temp_vec,NULL);CHKERRQ(ierr);
815             ierr = KSPCreate(PETSC_COMM_SELF,&temp_ksp);CHKERRQ(ierr);
816             ierr = KSPSetOperators(temp_ksp,temp_mat,temp_mat,SAME_PRECONDITIONER);CHKERRQ(ierr);
817             ierr = KSPSetType(temp_ksp,KSPPREONLY);CHKERRQ(ierr);
818             ierr = KSPGetPC(temp_ksp,&temp_pc);CHKERRQ(ierr);
819             ierr = PCSetType(temp_pc,PCLU);CHKERRQ(ierr);
820             ierr = KSPSetUp(temp_ksp);CHKERRQ(ierr);
821             for (s=0;s<temp_constraints;s++) {
822               ierr = VecSet(temp_vec,0.0);CHKERRQ(ierr);
823               ierr = VecSetValue(temp_vec,s+dual_dofs,1.0,INSERT_VALUES);CHKERRQ(ierr);
824               ierr = VecAssemblyBegin(temp_vec);CHKERRQ(ierr);
825               ierr = VecAssemblyEnd(temp_vec);CHKERRQ(ierr);
826               ierr = KSPSolve(temp_ksp,temp_vec,temp_vec);CHKERRQ(ierr);
827               ierr = VecGetArray(temp_vec,&array_vector);CHKERRQ(ierr);
828               j = temp_indices_to_constraint_B[temp_indices[start_constraint+s]+size_of_constraint-s-1];
829               /* last columns of change of basis matrix associated to new primal dofs */
830               ierr = MatSetValues(pcbddc->ChangeOfBasisMatrix,size_of_constraint,&temp_indices_to_constraint_B[temp_indices[start_constraint+s]],1,&j,array_vector,INSERT_VALUES);CHKERRQ(ierr);
831               ierr = VecRestoreArray(temp_vec,&array_vector);CHKERRQ(ierr);
832             }
833             ierr = MatDestroy(&temp_mat);CHKERRQ(ierr);
834             ierr = KSPDestroy(&temp_ksp);CHKERRQ(ierr);
835             ierr = VecDestroy(&temp_vec);CHKERRQ(ierr);
836           } else {
837             /* last columns of change of basis matrix associated to new primal dofs */
838             for (s=0;s<temp_constraints;s++) {
839               j = temp_indices_to_constraint_B[temp_indices[start_constraint+s]+size_of_constraint-s-1];
840               ierr = MatSetValues(pcbddc->ChangeOfBasisMatrix,size_of_constraint,&temp_indices_to_constraint_B[temp_indices[start_constraint+s]],1,&j,&temp_quadrature_constraint[temp_indices[start_constraint+s]],INSERT_VALUES);CHKERRQ(ierr);
841             }
842           }
843           /* prepare for the next cycle */
844           temp_constraints = 0;
845           if (i != local_primal_size -1 ) {
846             temp_start_ptr = temp_indices_to_constraint_B[temp_indices[i+1]];
847           }
848         }
849       }
850     }
851     /* assembling */
852     ierr = MatAssemblyBegin(pcbddc->ChangeOfBasisMatrix,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
853     ierr = MatAssemblyEnd(pcbddc->ChangeOfBasisMatrix,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
854     ierr = PetscFree(ipiv);CHKERRQ(ierr);
855     ierr = PetscFree(is_indices);CHKERRQ(ierr);
856   }
857   /* free workspace no longer needed */
858   ierr = PetscFree(rwork);CHKERRQ(ierr);
859   ierr = PetscFree(work);CHKERRQ(ierr);
860   ierr = PetscFree(temp_basis);CHKERRQ(ierr);
861   ierr = PetscFree(singular_vals);CHKERRQ(ierr);
862   ierr = PetscFree(correlation_mat);CHKERRQ(ierr);
863   ierr = PetscFree(temp_indices);CHKERRQ(ierr);
864   ierr = PetscFree(change_basis);CHKERRQ(ierr);
865   ierr = PetscFree(temp_indices_to_constraint);CHKERRQ(ierr);
866   ierr = PetscFree(temp_indices_to_constraint_B);CHKERRQ(ierr);
867   ierr = PetscFree(local_to_B);CHKERRQ(ierr);
868   ierr = PetscFree(temp_quadrature_constraint);CHKERRQ(ierr);
869 #if defined(PETSC_MISSING_LAPACK_GESVD)
870   ierr = PetscFree(iwork);CHKERRQ(ierr);
871   ierr = PetscFree(ifail);CHKERRQ(ierr);
872   ierr = PetscFree(singular_vectors);CHKERRQ(ierr);
873 #endif
874   for (k=0;k<nnsp_size;k++) {
875     ierr = VecDestroy(&localnearnullsp[k]);CHKERRQ(ierr);
876   }
877   ierr = PetscFree(localnearnullsp);CHKERRQ(ierr);
878   PetscFunctionReturn(0);
879 }
880 
881 #undef __FUNCT__
882 #define __FUNCT__ "PCBDDCAnalyzeInterface"
883 PetscErrorCode PCBDDCAnalyzeInterface(PC pc)
884 {
885   PC_BDDC     *pcbddc = (PC_BDDC*)pc->data;
886   PC_IS       *pcis = (PC_IS*)pc->data;
887   Mat_IS      *matis  = (Mat_IS*)pc->pmat->data;
888   PetscInt    bs,ierr,i,vertex_size;
889   PetscViewer viewer=pcbddc->dbg_viewer;
890 
891   PetscFunctionBegin;
892   /* Init local Graph struct */
893   ierr = PCBDDCGraphInit(pcbddc->mat_graph,matis->mapping);CHKERRQ(ierr);
894 
895   /* Check validity of the csr graph passed in by the user */
896   if (pcbddc->mat_graph->nvtxs_csr != pcbddc->mat_graph->nvtxs) {
897     ierr = PCBDDCGraphResetCSR(pcbddc->mat_graph);CHKERRQ(ierr);
898   }
899   /* Set default CSR adjacency of local dofs if not provided by the user with PCBDDCSetLocalAdjacencyGraph */
900   if (!pcbddc->mat_graph->xadj || !pcbddc->mat_graph->adjncy) {
901     Mat mat_adj;
902     const PetscInt *xadj,*adjncy;
903     PetscBool flg_row=PETSC_TRUE;
904 
905     ierr = MatConvert(matis->A,MATMPIADJ,MAT_INITIAL_MATRIX,&mat_adj);CHKERRQ(ierr);
906     ierr = MatGetRowIJ(mat_adj,0,PETSC_TRUE,PETSC_FALSE,&i,&xadj,&adjncy,&flg_row);CHKERRQ(ierr);
907     if (!flg_row) {
908       SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in MatGetRowIJ called in %s\n",__FUNCT__);
909     }
910     ierr = PCBDDCSetLocalAdjacencyGraph(pc,i,xadj,adjncy,PETSC_COPY_VALUES);CHKERRQ(ierr);
911     ierr = MatRestoreRowIJ(mat_adj,0,PETSC_TRUE,PETSC_FALSE,&i,&xadj,&adjncy,&flg_row);CHKERRQ(ierr);
912     if (!flg_row) {
913       SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called in %s\n",__FUNCT__);
914     }
915     ierr = MatDestroy(&mat_adj);CHKERRQ(ierr);
916   }
917 
918   /* Set default dofs' splitting if no information has been provided by the user with PCBDDCSetDofsSplitting */
919   vertex_size = 1;
920   if (!pcbddc->n_ISForDofs) {
921     IS *custom_ISForDofs;
922 
923     ierr = MatGetBlockSize(matis->A,&bs);CHKERRQ(ierr);
924     ierr = PetscMalloc(bs*sizeof(IS),&custom_ISForDofs);CHKERRQ(ierr);
925     for (i=0;i<bs;i++) {
926       ierr = ISCreateStride(PETSC_COMM_SELF,pcis->n/bs,i,bs,&custom_ISForDofs[i]);CHKERRQ(ierr);
927     }
928     ierr = PCBDDCSetDofsSplitting(pc,bs,custom_ISForDofs);CHKERRQ(ierr);
929     /* remove my references to IS objects */
930     for (i=0;i<bs;i++) {
931       ierr = ISDestroy(&custom_ISForDofs[i]);CHKERRQ(ierr);
932     }
933     ierr = PetscFree(custom_ISForDofs);CHKERRQ(ierr);
934   } else { /* mat block size as vertex size (used for elasticity) */
935     ierr = MatGetBlockSize(matis->A,&vertex_size);CHKERRQ(ierr);
936   }
937 
938   /* Setup of Graph */
939   ierr = PCBDDCGraphSetUp(pcbddc->mat_graph,vertex_size,pcbddc->NeumannBoundaries,pcbddc->DirichletBoundaries,pcbddc->n_ISForDofs,pcbddc->ISForDofs,pcbddc->user_primal_vertices);
940 
941   /* Graph's connected components analysis */
942   ierr = PCBDDCGraphComputeConnectedComponents(pcbddc->mat_graph);CHKERRQ(ierr);
943 
944   /* print some info to stdout */
945   if (pcbddc->dbg_flag) {
946     ierr = PCBDDCGraphASCIIView(pcbddc->mat_graph,pcbddc->dbg_flag,viewer);
947   }
948   PetscFunctionReturn(0);
949 }
950 
951 #undef __FUNCT__
952 #define __FUNCT__ "PCBDDCGetPrimalVerticesLocalIdx"
953 PetscErrorCode  PCBDDCGetPrimalVerticesLocalIdx(PC pc, PetscInt *n_vertices, PetscInt *vertices_idx[])
954 {
955   PC_BDDC        *pcbddc = (PC_BDDC*)(pc->data);
956   PetscInt       *vertices,*row_cmat_indices,n,i,size_of_constraint,local_primal_size;
957   PetscErrorCode ierr;
958 
959   PetscFunctionBegin;
960   n = 0;
961   vertices = 0;
962   if (pcbddc->ConstraintMatrix) {
963     ierr = MatGetSize(pcbddc->ConstraintMatrix,&local_primal_size,&i);CHKERRQ(ierr);
964     for (i=0;i<local_primal_size;i++) {
965       ierr = MatGetRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,NULL,NULL);CHKERRQ(ierr);
966       if (size_of_constraint == 1) n++;
967       ierr = MatRestoreRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,NULL,NULL);CHKERRQ(ierr);
968     }
969     ierr = PetscMalloc(n*sizeof(PetscInt),&vertices);CHKERRQ(ierr);
970     n = 0;
971     for (i=0;i<local_primal_size;i++) {
972       ierr = MatGetRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,NULL);CHKERRQ(ierr);
973       if (size_of_constraint == 1) {
974         vertices[n++]=row_cmat_indices[0];
975       }
976       ierr = MatRestoreRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,NULL);CHKERRQ(ierr);
977     }
978   }
979   *n_vertices = n;
980   *vertices_idx = vertices;
981   PetscFunctionReturn(0);
982 }
983 
984 #undef __FUNCT__
985 #define __FUNCT__ "PCBDDCGetPrimalConstraintsLocalIdx"
986 PetscErrorCode  PCBDDCGetPrimalConstraintsLocalIdx(PC pc, PetscInt *n_constraints, PetscInt *constraints_idx[])
987 {
988   PC_BDDC        *pcbddc = (PC_BDDC*)(pc->data);
989   PetscInt       *constraints_index,*row_cmat_indices,*row_cmat_global_indices;
990   PetscInt       n,i,j,size_of_constraint,local_primal_size,local_size,max_size_of_constraint,min_index,min_loc;
991   PetscBool      *touched;
992   PetscErrorCode ierr;
993 
994   PetscFunctionBegin;
995   n = 0;
996   constraints_index = 0;
997   if (pcbddc->ConstraintMatrix) {
998     ierr = MatGetSize(pcbddc->ConstraintMatrix,&local_primal_size,&local_size);CHKERRQ(ierr);
999     max_size_of_constraint = 0;
1000     for (i=0;i<local_primal_size;i++) {
1001       ierr = MatGetRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,NULL,NULL);CHKERRQ(ierr);
1002       if (size_of_constraint > 1) {
1003         n++;
1004       }
1005       max_size_of_constraint = PetscMax(size_of_constraint,max_size_of_constraint);
1006       ierr = MatRestoreRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,NULL,NULL);CHKERRQ(ierr);
1007     }
1008     ierr = PetscMalloc(n*sizeof(PetscInt),&constraints_index);CHKERRQ(ierr);
1009     ierr = PetscMalloc(max_size_of_constraint*sizeof(PetscInt),&row_cmat_global_indices);CHKERRQ(ierr);
1010     ierr = PetscMalloc(local_size*sizeof(PetscBool),&touched);CHKERRQ(ierr);
1011     ierr = PetscMemzero(touched,local_size*sizeof(PetscBool));CHKERRQ(ierr);
1012     n = 0;
1013     for (i=0;i<local_primal_size;i++) {
1014       ierr = MatGetRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,NULL);CHKERRQ(ierr);
1015       if (size_of_constraint > 1) {
1016         ierr = ISLocalToGlobalMappingApply(pcbddc->mat_graph->l2gmap,size_of_constraint,row_cmat_indices,row_cmat_global_indices);CHKERRQ(ierr);
1017         /* find first untouched local node */
1018         j = 0;
1019         while(touched[row_cmat_indices[j]]) j++;
1020         min_index = row_cmat_global_indices[j];
1021         min_loc = j;
1022         /* search the minimum among nodes not yet touched on the connected component
1023            since there can be more than one constraint on a single cc */
1024         for (j=1;j<size_of_constraint;j++) {
1025           if (min_index > row_cmat_global_indices[j] && !touched[row_cmat_indices[j]]) {
1026             min_index = row_cmat_global_indices[j];
1027             min_loc = j;
1028           }
1029         }
1030         touched[row_cmat_indices[min_loc]] = PETSC_TRUE;
1031         constraints_index[n++] = row_cmat_indices[min_loc];
1032       }
1033       ierr = MatRestoreRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,NULL);CHKERRQ(ierr);
1034     }
1035   }
1036   ierr = PetscFree(touched);CHKERRQ(ierr);
1037   ierr = PetscFree(row_cmat_global_indices);CHKERRQ(ierr);
1038   *n_constraints = n;
1039   *constraints_idx = constraints_index;
1040   PetscFunctionReturn(0);
1041 }
1042 
1043 /* the next two functions has been adapted from pcis.c */
1044 #undef __FUNCT__
1045 #define __FUNCT__ "PCBDDCApplySchur"
1046 PetscErrorCode  PCBDDCApplySchur(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D)
1047 {
1048   PetscErrorCode ierr;
1049   PC_IS          *pcis = (PC_IS*)(pc->data);
1050 
1051   PetscFunctionBegin;
1052   if (!vec2_B) { vec2_B = v; }
1053   ierr = MatMult(pcis->A_BB,v,vec1_B);CHKERRQ(ierr);
1054   ierr = MatMult(pcis->A_IB,v,vec1_D);CHKERRQ(ierr);
1055   ierr = KSPSolve(pcis->ksp_D,vec1_D,vec2_D);CHKERRQ(ierr);
1056   ierr = MatMult(pcis->A_BI,vec2_D,vec2_B);CHKERRQ(ierr);
1057   ierr = VecAXPY(vec1_B,-1.0,vec2_B);CHKERRQ(ierr);
1058   PetscFunctionReturn(0);
1059 }
1060 
1061 #undef __FUNCT__
1062 #define __FUNCT__ "PCBDDCApplySchurTranspose"
1063 PetscErrorCode  PCBDDCApplySchurTranspose(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D)
1064 {
1065   PetscErrorCode ierr;
1066   PC_IS          *pcis = (PC_IS*)(pc->data);
1067 
1068   PetscFunctionBegin;
1069   if (!vec2_B) { vec2_B = v; }
1070   ierr = MatMultTranspose(pcis->A_BB,v,vec1_B);CHKERRQ(ierr);
1071   ierr = MatMultTranspose(pcis->A_BI,v,vec1_D);CHKERRQ(ierr);
1072   ierr = KSPSolveTranspose(pcis->ksp_D,vec1_D,vec2_D);CHKERRQ(ierr);
1073   ierr = MatMultTranspose(pcis->A_IB,vec2_D,vec2_B);CHKERRQ(ierr);
1074   ierr = VecAXPY(vec1_B,-1.0,vec2_B);CHKERRQ(ierr);
1075   PetscFunctionReturn(0);
1076 }
1077 
1078 #undef __FUNCT__
1079 #define __FUNCT__ "PCBDDCSubsetNumbering"
1080 PetscErrorCode PCBDDCSubsetNumbering(MPI_Comm comm,ISLocalToGlobalMapping l2gmap, PetscInt n_local_dofs, PetscInt local_dofs[], PetscInt local_dofs_mult[], PetscInt* n_global_subset, PetscInt* global_numbering_subset[])
1081 {
1082   Vec            local_vec,global_vec;
1083   IS             seqis,paris;
1084   VecScatter     scatter_ctx;
1085   PetscScalar    *array;
1086   PetscInt       *temp_global_dofs;
1087   PetscScalar    globalsum;
1088   PetscInt       i,j,s;
1089   PetscInt       nlocals,first_index,old_index,max_local;
1090   PetscMPIInt    rank_prec_comm,size_prec_comm,max_global;
1091   PetscMPIInt    *dof_sizes,*dof_displs;
1092   PetscBool      first_found;
1093   PetscErrorCode ierr;
1094 
1095   PetscFunctionBegin;
1096   /* mpi buffers */
1097   MPI_Comm_size(comm,&size_prec_comm);
1098   MPI_Comm_rank(comm,&rank_prec_comm);
1099   j = ( !rank_prec_comm ? size_prec_comm : 0);
1100   ierr = PetscMalloc(j*sizeof(*dof_sizes),&dof_sizes);CHKERRQ(ierr);
1101   ierr = PetscMalloc(j*sizeof(*dof_displs),&dof_displs);CHKERRQ(ierr);
1102   /* get maximum size of subset */
1103   ierr = PetscMalloc(n_local_dofs*sizeof(PetscInt),&temp_global_dofs);CHKERRQ(ierr);
1104   ierr = ISLocalToGlobalMappingApply(l2gmap,n_local_dofs,local_dofs,temp_global_dofs);CHKERRQ(ierr);
1105   max_local = 0;
1106   if (n_local_dofs) {
1107     max_local = temp_global_dofs[0];
1108     for (i=1;i<n_local_dofs;i++) {
1109       if (max_local < temp_global_dofs[i] ) {
1110         max_local = temp_global_dofs[i];
1111       }
1112     }
1113   }
1114   ierr = MPI_Allreduce(&max_local,&max_global,1,MPIU_INT,MPI_MAX,comm);
1115   max_global++;
1116   max_local = 0;
1117   if (n_local_dofs) {
1118     max_local = local_dofs[0];
1119     for (i=1;i<n_local_dofs;i++) {
1120       if (max_local < local_dofs[i] ) {
1121         max_local = local_dofs[i];
1122       }
1123     }
1124   }
1125   max_local++;
1126   /* allocate workspace */
1127   ierr = VecCreate(PETSC_COMM_SELF,&local_vec);CHKERRQ(ierr);
1128   ierr = VecSetSizes(local_vec,PETSC_DECIDE,max_local);CHKERRQ(ierr);
1129   ierr = VecSetType(local_vec,VECSEQ);CHKERRQ(ierr);
1130   ierr = VecCreate(comm,&global_vec);CHKERRQ(ierr);
1131   ierr = VecSetSizes(global_vec,PETSC_DECIDE,max_global);CHKERRQ(ierr);
1132   ierr = VecSetType(global_vec,VECMPI);CHKERRQ(ierr);
1133   /* create scatter */
1134   ierr = ISCreateGeneral(PETSC_COMM_SELF,n_local_dofs,local_dofs,PETSC_COPY_VALUES,&seqis);CHKERRQ(ierr);
1135   ierr = ISCreateGeneral(comm,n_local_dofs,temp_global_dofs,PETSC_COPY_VALUES,&paris);CHKERRQ(ierr);
1136   ierr = VecScatterCreate(local_vec,seqis,global_vec,paris,&scatter_ctx);CHKERRQ(ierr);
1137   ierr = ISDestroy(&seqis);CHKERRQ(ierr);
1138   ierr = ISDestroy(&paris);CHKERRQ(ierr);
1139   /* init array */
1140   ierr = VecSet(global_vec,0.0);CHKERRQ(ierr);
1141   ierr = VecSet(local_vec,0.0);CHKERRQ(ierr);
1142   ierr = VecGetArray(local_vec,&array);CHKERRQ(ierr);
1143   if (local_dofs_mult) {
1144     for (i=0;i<n_local_dofs;i++) {
1145       array[local_dofs[i]]=(PetscScalar)local_dofs_mult[i];
1146     }
1147   } else {
1148     for (i=0;i<n_local_dofs;i++) {
1149       array[local_dofs[i]]=1.0;
1150     }
1151   }
1152   ierr = VecRestoreArray(local_vec,&array);CHKERRQ(ierr);
1153   /* scatter into global vec and get total number of global dofs */
1154   ierr = VecScatterBegin(scatter_ctx,local_vec,global_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1155   ierr = VecScatterEnd(scatter_ctx,local_vec,global_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1156   ierr = VecSum(global_vec,&globalsum);CHKERRQ(ierr);
1157   *n_global_subset = (PetscInt)PetscRealPart(globalsum);
1158   /* Fill global_vec with cumulative function for global numbering */
1159   ierr = VecGetArray(global_vec,&array);CHKERRQ(ierr);
1160   ierr = VecGetLocalSize(global_vec,&s);CHKERRQ(ierr);
1161   nlocals = 0;
1162   first_index = -1;
1163   first_found = PETSC_FALSE;
1164   for (i=0;i<s;i++) {
1165     if (!first_found && PetscRealPart(array[i]) > 0.0) {
1166       first_found = PETSC_TRUE;
1167       first_index = i;
1168     }
1169     nlocals += (PetscInt)PetscRealPart(array[i]);
1170   }
1171   ierr = MPI_Gather(&nlocals,1,MPIU_INT,dof_sizes,1,MPIU_INT,0,comm);CHKERRQ(ierr);
1172   if (!rank_prec_comm) {
1173     dof_displs[0]=0;
1174     for (i=1;i<size_prec_comm;i++) {
1175       dof_displs[i] = dof_displs[i-1]+dof_sizes[i-1];
1176     }
1177   }
1178   ierr = MPI_Scatter(dof_displs,1,MPIU_INT,&nlocals,1,MPIU_INT,0,comm);CHKERRQ(ierr);
1179   if (first_found) {
1180     array[first_index] += (PetscScalar)nlocals;
1181     old_index = first_index;
1182     for (i=first_index+1;i<s;i++) {
1183       if (PetscRealPart(array[i]) > 0.0) {
1184         array[i] += array[old_index];
1185         old_index = i;
1186       }
1187     }
1188   }
1189   ierr = VecRestoreArray(global_vec,&array);CHKERRQ(ierr);
1190   ierr = VecSet(local_vec,0.0);CHKERRQ(ierr);
1191   ierr = VecScatterBegin(scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1192   ierr = VecScatterEnd  (scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
1193   /* get global ordering of local dofs */
1194   ierr = VecGetArray(local_vec,&array);CHKERRQ(ierr);
1195   if (local_dofs_mult) {
1196     for (i=0;i<n_local_dofs;i++) {
1197       temp_global_dofs[i] = (PetscInt)PetscRealPart(array[local_dofs[i]])-local_dofs_mult[i];
1198     }
1199   } else {
1200     for (i=0;i<n_local_dofs;i++) {
1201       temp_global_dofs[i] = (PetscInt)PetscRealPart(array[local_dofs[i]])-1;
1202     }
1203   }
1204   ierr = VecRestoreArray(local_vec,&array);CHKERRQ(ierr);
1205   /* free workspace */
1206   ierr = VecScatterDestroy(&scatter_ctx);CHKERRQ(ierr);
1207   ierr = VecDestroy(&local_vec);CHKERRQ(ierr);
1208   ierr = VecDestroy(&global_vec);CHKERRQ(ierr);
1209   ierr = PetscFree(dof_sizes);CHKERRQ(ierr);
1210   ierr = PetscFree(dof_displs);CHKERRQ(ierr);
1211   /* return pointer to global ordering of local dofs */
1212   *global_numbering_subset = temp_global_dofs;
1213   PetscFunctionReturn(0);
1214 }
1215 
1216 #undef __FUNCT__
1217 #define __FUNCT__ "PCBDDCOrthonormalizeVecs"
1218 PetscErrorCode PCBDDCOrthonormalizeVecs(PetscInt n, Vec vecs[])
1219 {
1220   PetscInt       i,j;
1221   PetscScalar    *alphas;
1222   PetscErrorCode ierr;
1223 
1224   PetscFunctionBegin;
1225   /* this implements stabilized Gram-Schmidt */
1226   ierr = PetscMalloc(n*sizeof(PetscScalar),&alphas);CHKERRQ(ierr);
1227   for (i=0;i<n;i++) {
1228     ierr = VecNormalize(vecs[i],NULL);CHKERRQ(ierr);
1229     if (i<n) { ierr = VecMDot(vecs[i],n-i-1,&vecs[i+1],&alphas[i+1]);CHKERRQ(ierr); }
1230     for (j=i+1;j<n;j++) { ierr = VecAXPY(vecs[j],PetscConj(-alphas[j]),vecs[i]);CHKERRQ(ierr); }
1231   }
1232   ierr = PetscFree(alphas);CHKERRQ(ierr);
1233   PetscFunctionReturn(0);
1234 }
1235 
1236 
1237