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