1 #include <petsc/private/petscimpl.h>
2 #include <../src/ksp/pc/impls/bddc/bddcprivate.h>
3 #include <../src/ksp/pc/impls/bddc/bddcstructs.h>
4 
PCBDDCGraphGetDirichletDofsB(PCBDDCGraph graph,IS * dirdofs)5 PetscErrorCode PCBDDCGraphGetDirichletDofsB(PCBDDCGraph graph, IS* dirdofs)
6 {
7   PetscErrorCode ierr;
8 
9   PetscFunctionBegin;
10   if (graph->dirdofsB) {
11     ierr = PetscObjectReference((PetscObject)graph->dirdofsB);CHKERRQ(ierr);
12   } else if (graph->has_dirichlet) {
13     PetscInt i,size;
14     PetscInt *dirdofs_idxs;
15 
16     size = 0;
17     for (i=0;i<graph->nvtxs;i++) {
18       if (graph->count[i] && graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK) size++;
19     }
20 
21     ierr = PetscMalloc1(size,&dirdofs_idxs);CHKERRQ(ierr);
22     size = 0;
23     for (i=0;i<graph->nvtxs;i++) {
24       if (graph->count[i] && graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK) dirdofs_idxs[size++] = i;
25     }
26     ierr = ISCreateGeneral(PETSC_COMM_SELF,size,dirdofs_idxs,PETSC_OWN_POINTER,&graph->dirdofsB);CHKERRQ(ierr);
27     ierr = PetscObjectReference((PetscObject)graph->dirdofsB);CHKERRQ(ierr);
28   }
29   *dirdofs = graph->dirdofsB;
30   PetscFunctionReturn(0);
31 }
32 
PCBDDCGraphGetDirichletDofs(PCBDDCGraph graph,IS * dirdofs)33 PetscErrorCode PCBDDCGraphGetDirichletDofs(PCBDDCGraph graph, IS* dirdofs)
34 {
35   PetscErrorCode ierr;
36 
37   PetscFunctionBegin;
38   if (graph->dirdofs) {
39     ierr = PetscObjectReference((PetscObject)graph->dirdofs);CHKERRQ(ierr);
40   } else if (graph->has_dirichlet) {
41     PetscInt i,size;
42     PetscInt *dirdofs_idxs;
43 
44     size = 0;
45     for (i=0;i<graph->nvtxs;i++) {
46       if (graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK) size++;
47     }
48 
49     ierr = PetscMalloc1(size,&dirdofs_idxs);CHKERRQ(ierr);
50     size = 0;
51     for (i=0;i<graph->nvtxs;i++) {
52       if (graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK) dirdofs_idxs[size++] = i;
53     }
54     ierr = ISCreateGeneral(PetscObjectComm((PetscObject)graph->l2gmap),size,dirdofs_idxs,PETSC_OWN_POINTER,&graph->dirdofs);CHKERRQ(ierr);
55     ierr = PetscObjectReference((PetscObject)graph->dirdofs);CHKERRQ(ierr);
56   }
57   *dirdofs = graph->dirdofs;
58   PetscFunctionReturn(0);
59 }
60 
PCBDDCGraphASCIIView(PCBDDCGraph graph,PetscInt verbosity_level,PetscViewer viewer)61 PetscErrorCode PCBDDCGraphASCIIView(PCBDDCGraph graph, PetscInt verbosity_level, PetscViewer viewer)
62 {
63   PetscInt       i,j,tabs;
64   PetscInt*      queue_in_global_numbering;
65   PetscErrorCode ierr;
66 
67   PetscFunctionBegin;
68   ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr);
69   ierr = PetscViewerASCIIGetTab(viewer,&tabs);CHKERRQ(ierr);
70   ierr = PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");CHKERRQ(ierr);
71   ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
72   ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Local BDDC graph for subdomain %04d\n",PetscGlobalRank);CHKERRQ(ierr);
73   ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Number of vertices %d\n",graph->nvtxs);CHKERRQ(ierr);
74   ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Custom minimal size %d\n",graph->custom_minimal_size);CHKERRQ(ierr);
75   if (graph->maxcount != PETSC_MAX_INT) {
76     ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Max count %d\n",graph->maxcount);CHKERRQ(ierr);
77   }
78   ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Topological two dim? %d (set %d)\n",graph->twodim,graph->twodimset);CHKERRQ(ierr);
79   if (verbosity_level > 2) {
80     for (i=0;i<graph->nvtxs;i++) {
81       ierr = PetscViewerASCIISynchronizedPrintf(viewer,"%d:\n",i);CHKERRQ(ierr);
82       ierr = PetscViewerASCIISynchronizedPrintf(viewer,"   which_dof: %d\n",graph->which_dof[i]);CHKERRQ(ierr);
83       ierr = PetscViewerASCIISynchronizedPrintf(viewer,"   special_dof: %d\n",graph->special_dof[i]);CHKERRQ(ierr);
84       ierr = PetscViewerASCIISynchronizedPrintf(viewer,"   neighbours: %d\n",graph->count[i]);CHKERRQ(ierr);
85       ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr);
86       if (graph->count[i]) {
87         ierr = PetscViewerASCIISynchronizedPrintf(viewer,"     set of neighbours:");CHKERRQ(ierr);
88         for (j=0;j<graph->count[i];j++) {
89           ierr = PetscViewerASCIISynchronizedPrintf(viewer," %d",graph->neighbours_set[i][j]);CHKERRQ(ierr);
90         }
91         ierr = PetscViewerASCIISynchronizedPrintf(viewer,"\n");CHKERRQ(ierr);
92       }
93       ierr = PetscViewerASCIISetTab(viewer,tabs);CHKERRQ(ierr);
94       ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr);
95       if (graph->mirrors) {
96         ierr = PetscViewerASCIISynchronizedPrintf(viewer,"   mirrors: %d\n",graph->mirrors[i]);CHKERRQ(ierr);
97         if (graph->mirrors[i]) {
98           ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr);
99           ierr = PetscViewerASCIISynchronizedPrintf(viewer,"     set of mirrors:");CHKERRQ(ierr);
100           for (j=0;j<graph->mirrors[i];j++) {
101             ierr = PetscViewerASCIISynchronizedPrintf(viewer," %d",graph->mirrors_set[i][j]);CHKERRQ(ierr);
102           }
103           ierr = PetscViewerASCIISynchronizedPrintf(viewer,"\n");CHKERRQ(ierr);
104           ierr = PetscViewerASCIISetTab(viewer,tabs);CHKERRQ(ierr);
105           ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr);
106         }
107       }
108       if (verbosity_level > 3) {
109         if (graph->xadj) {
110           ierr = PetscViewerASCIISynchronizedPrintf(viewer,"   local adj list:");CHKERRQ(ierr);
111           ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr);
112           for (j=graph->xadj[i];j<graph->xadj[i+1];j++) {
113             ierr = PetscViewerASCIISynchronizedPrintf(viewer," %d",graph->adjncy[j]);CHKERRQ(ierr);
114           }
115           ierr = PetscViewerASCIISynchronizedPrintf(viewer,"\n");CHKERRQ(ierr);
116           ierr = PetscViewerASCIISetTab(viewer,tabs);CHKERRQ(ierr);
117           ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr);
118         } else {
119           ierr = PetscViewerASCIISynchronizedPrintf(viewer,"   no adj info\n");CHKERRQ(ierr);
120         }
121       }
122       if (graph->n_local_subs) {
123         ierr = PetscViewerASCIISynchronizedPrintf(viewer,"   local sub id: %d\n",graph->local_subs[i]);CHKERRQ(ierr);
124       }
125       ierr = PetscViewerASCIISynchronizedPrintf(viewer,"   interface subset id: %d\n",graph->subset[i]);CHKERRQ(ierr);
126       if (graph->subset[i] && graph->subset_ncc) {
127         ierr = PetscViewerASCIISynchronizedPrintf(viewer,"   ncc for subset: %d\n",graph->subset_ncc[graph->subset[i]-1]);CHKERRQ(ierr);
128       }
129     }
130   }
131   ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Total number of connected components %d\n",graph->ncc);CHKERRQ(ierr);
132   ierr = PetscMalloc1(graph->cptr[graph->ncc],&queue_in_global_numbering);CHKERRQ(ierr);
133   ierr = ISLocalToGlobalMappingApply(graph->l2gmap,graph->cptr[graph->ncc],graph->queue,queue_in_global_numbering);CHKERRQ(ierr);
134   for (i=0;i<graph->ncc;i++) {
135     PetscInt node_num=graph->queue[graph->cptr[i]];
136     PetscBool printcc = PETSC_FALSE;
137     ierr = PetscViewerASCIISynchronizedPrintf(viewer,"  cc %d (size %d, fid %d, neighs:",i,graph->cptr[i+1]-graph->cptr[i],graph->which_dof[node_num]);CHKERRQ(ierr);
138     ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr);
139     for (j=0;j<graph->count[node_num];j++) {
140       ierr = PetscViewerASCIISynchronizedPrintf(viewer," %d",graph->neighbours_set[node_num][j]);CHKERRQ(ierr);
141     }
142     if (verbosity_level > 1) {
143       ierr = PetscViewerASCIISynchronizedPrintf(viewer,"):");CHKERRQ(ierr);
144       if (verbosity_level > 2 || graph->twodim || graph->count[node_num] > 1 || (graph->count[node_num] == 1 && graph->special_dof[node_num] == PCBDDCGRAPH_NEUMANN_MARK)) {
145         printcc = PETSC_TRUE;
146       }
147       if (printcc) {
148         for (j=graph->cptr[i];j<graph->cptr[i+1];j++) {
149           ierr = PetscViewerASCIISynchronizedPrintf(viewer," %d (%d)",graph->queue[j],queue_in_global_numbering[j]);CHKERRQ(ierr);
150         }
151       }
152     } else {
153       ierr = PetscViewerASCIISynchronizedPrintf(viewer,")");CHKERRQ(ierr);
154     }
155     ierr = PetscViewerASCIISynchronizedPrintf(viewer,"\n");CHKERRQ(ierr);
156     ierr = PetscViewerASCIISetTab(viewer,tabs);CHKERRQ(ierr);
157     ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr);
158   }
159   ierr = PetscFree(queue_in_global_numbering);CHKERRQ(ierr);
160   ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
161   PetscFunctionReturn(0);
162 }
163 
PCBDDCGraphRestoreCandidatesIS(PCBDDCGraph graph,PetscInt * n_faces,IS * FacesIS[],PetscInt * n_edges,IS * EdgesIS[],IS * VerticesIS)164 PetscErrorCode PCBDDCGraphRestoreCandidatesIS(PCBDDCGraph graph, PetscInt *n_faces, IS *FacesIS[], PetscInt *n_edges, IS *EdgesIS[], IS *VerticesIS)
165 {
166   PetscInt       i;
167   PetscErrorCode ierr;
168 
169   PetscFunctionBegin;
170   if (n_faces) {
171     if (FacesIS) {
172       for (i=0;i<*n_faces;i++) {
173         ierr = ISDestroy(&((*FacesIS)[i]));CHKERRQ(ierr);
174       }
175       ierr = PetscFree(*FacesIS);CHKERRQ(ierr);
176     }
177     *n_faces = 0;
178   }
179   if (n_edges) {
180     if (EdgesIS) {
181       for (i=0;i<*n_edges;i++) {
182         ierr = ISDestroy(&((*EdgesIS)[i]));CHKERRQ(ierr);
183       }
184       ierr = PetscFree(*EdgesIS);CHKERRQ(ierr);
185     }
186     *n_edges = 0;
187   }
188   if (VerticesIS) {
189     ierr = ISDestroy(VerticesIS);CHKERRQ(ierr);
190   }
191   PetscFunctionReturn(0);
192 }
193 
PCBDDCGraphGetCandidatesIS(PCBDDCGraph graph,PetscInt * n_faces,IS * FacesIS[],PetscInt * n_edges,IS * EdgesIS[],IS * VerticesIS)194 PetscErrorCode PCBDDCGraphGetCandidatesIS(PCBDDCGraph graph, PetscInt *n_faces, IS *FacesIS[], PetscInt *n_edges, IS *EdgesIS[], IS *VerticesIS)
195 {
196   IS             *ISForFaces,*ISForEdges,ISForVertices;
197   PetscInt       i,nfc,nec,nvc,*idx,*mark;
198   PetscErrorCode ierr;
199 
200   PetscFunctionBegin;
201   ierr = PetscCalloc1(graph->ncc,&mark);CHKERRQ(ierr);
202   /* loop on ccs to evalute number of faces, edges and vertices */
203   nfc = 0;
204   nec = 0;
205   nvc = 0;
206   for (i=0;i<graph->ncc;i++) {
207     PetscInt repdof = graph->queue[graph->cptr[i]];
208     if (graph->cptr[i+1]-graph->cptr[i] > graph->custom_minimal_size && graph->count[repdof] < graph->maxcount) {
209       if (!graph->twodim && graph->count[repdof] == 1 && graph->special_dof[repdof] != PCBDDCGRAPH_NEUMANN_MARK) {
210         nfc++;
211         mark[i] = 2;
212       } else {
213         nec++;
214         mark[i] = 1;
215       }
216     } else {
217       nvc += graph->cptr[i+1]-graph->cptr[i];
218     }
219   }
220 
221   /* allocate IS arrays for faces, edges. Vertices need a single index set. */
222   if (FacesIS) {
223     ierr = PetscMalloc1(nfc,&ISForFaces);CHKERRQ(ierr);
224   }
225   if (EdgesIS) {
226     ierr = PetscMalloc1(nec,&ISForEdges);CHKERRQ(ierr);
227   }
228   if (VerticesIS) {
229     ierr = PetscMalloc1(nvc,&idx);CHKERRQ(ierr);
230   }
231 
232   /* loop on ccs to compute index sets for faces and edges */
233   if (!graph->queue_sorted) {
234     PetscInt *queue_global;
235 
236     ierr = PetscMalloc1(graph->cptr[graph->ncc],&queue_global);CHKERRQ(ierr);
237     ierr = ISLocalToGlobalMappingApply(graph->l2gmap,graph->cptr[graph->ncc],graph->queue,queue_global);CHKERRQ(ierr);
238     for (i=0;i<graph->ncc;i++) {
239       ierr = PetscSortIntWithArray(graph->cptr[i+1]-graph->cptr[i],&queue_global[graph->cptr[i]],&graph->queue[graph->cptr[i]]);CHKERRQ(ierr);
240     }
241     ierr = PetscFree(queue_global);CHKERRQ(ierr);
242     graph->queue_sorted = PETSC_TRUE;
243   }
244   nfc = 0;
245   nec = 0;
246   for (i=0;i<graph->ncc;i++) {
247     if (mark[i] == 2) {
248       if (FacesIS) {
249         ierr = ISCreateGeneral(PETSC_COMM_SELF,graph->cptr[i+1]-graph->cptr[i],&graph->queue[graph->cptr[i]],PETSC_USE_POINTER,&ISForFaces[nfc]);CHKERRQ(ierr);
250       }
251       nfc++;
252     } else if (mark[i] == 1) {
253       if (EdgesIS) {
254         ierr = ISCreateGeneral(PETSC_COMM_SELF,graph->cptr[i+1]-graph->cptr[i],&graph->queue[graph->cptr[i]],PETSC_USE_POINTER,&ISForEdges[nec]);CHKERRQ(ierr);
255       }
256       nec++;
257     }
258   }
259 
260   /* index set for vertices */
261   if (VerticesIS) {
262     nvc = 0;
263     for (i=0;i<graph->ncc;i++) {
264       if (!mark[i]) {
265         PetscInt j;
266 
267         for (j=graph->cptr[i];j<graph->cptr[i+1];j++) {
268           idx[nvc]=graph->queue[j];
269           nvc++;
270         }
271       }
272     }
273     /* sort vertex set (by local ordering) */
274     ierr = PetscSortInt(nvc,idx);CHKERRQ(ierr);
275     ierr = ISCreateGeneral(PETSC_COMM_SELF,nvc,idx,PETSC_OWN_POINTER,&ISForVertices);CHKERRQ(ierr);
276   }
277   ierr = PetscFree(mark);CHKERRQ(ierr);
278 
279   /* get back info */
280   if (n_faces)       *n_faces = nfc;
281   if (FacesIS)       *FacesIS = ISForFaces;
282   if (n_edges)       *n_edges = nec;
283   if (EdgesIS)       *EdgesIS = ISForEdges;
284   if (VerticesIS) *VerticesIS = ISForVertices;
285   PetscFunctionReturn(0);
286 }
287 
PCBDDCGraphComputeConnectedComponents(PCBDDCGraph graph)288 PetscErrorCode PCBDDCGraphComputeConnectedComponents(PCBDDCGraph graph)
289 {
290   PetscBool      adapt_interface_reduced;
291   MPI_Comm       interface_comm;
292   PetscMPIInt    size;
293   PetscInt       i;
294   PetscBT        cornerp;
295   PetscErrorCode ierr;
296 
297   PetscFunctionBegin;
298   /* compute connected components locally */
299   ierr = PetscObjectGetComm((PetscObject)(graph->l2gmap),&interface_comm);CHKERRQ(ierr);
300   ierr = PCBDDCGraphComputeConnectedComponentsLocal(graph);CHKERRQ(ierr);
301 
302   cornerp = NULL;
303   if (graph->active_coords) { /* face based corner selection */
304     PetscBT   excluded;
305     PetscReal *wdist;
306     PetscInt  n_neigh,*neigh,*n_shared,**shared;
307     PetscInt  maxc, ns;
308 
309     ierr = PetscBTCreate(graph->nvtxs,&cornerp);CHKERRQ(ierr);
310     ierr = ISLocalToGlobalMappingGetInfo(graph->l2gmap,&n_neigh,&neigh,&n_shared,&shared);CHKERRQ(ierr);
311     for (ns = 1, maxc = 0; ns < n_neigh; ns++) maxc = PetscMax(maxc,n_shared[ns]);
312     ierr = PetscMalloc1(maxc*graph->cdim,&wdist);CHKERRQ(ierr);
313     ierr = PetscBTCreate(maxc,&excluded);CHKERRQ(ierr);
314 
315     for (ns = 1; ns < n_neigh; ns++) { /* first proc is self */
316       PetscReal *anchor,mdist;
317       PetscInt  fst,j,k,d,cdim = graph->cdim,n = n_shared[ns];
318       PetscInt  point1,point2,point3;
319 
320       /* import coordinates on shared interface */
321       ierr = PetscBTMemzero(n,excluded);CHKERRQ(ierr);
322       for (j=0,fst=-1,k=0;j<n;j++) {
323         PetscBool skip = PETSC_FALSE;
324         for (d=0;d<cdim;d++) {
325           PetscReal c = graph->coords[shared[ns][j]*cdim+d];
326           skip = (PetscBool)(skip || c == PETSC_MAX_REAL);
327           wdist[k++] = c;
328         }
329         if (skip) {
330           ierr = PetscBTSet(excluded,j);CHKERRQ(ierr);
331         } else if (fst == -1) fst = j;
332       }
333       if (fst == -1) continue;
334 
335       /* the dofs are sorted by global numbering, so each rank start from the same id and will detect the same corners from the given set */
336       anchor = wdist + fst*cdim;
337 
338       /* find the farthest point from the starting one */
339       mdist  = -1.0;
340       point1 = fst;
341       for (j=fst;j<n;j++) {
342         PetscReal dist = 0.0;
343 
344         if (PetscUnlikely(PetscBTLookup(excluded,j))) continue;
345         for (d=0;d<cdim;d++) dist += (wdist[j*cdim+d]-anchor[d])*(wdist[j*cdim+d]-anchor[d]);
346         if (dist > mdist) { mdist = dist; point1 = j; }
347       }
348 
349       /* find the farthest point from point1 */
350       anchor = wdist + point1*cdim;
351       mdist  = -1.0;
352       point2 = point1;
353       for (j=fst;j<n;j++) {
354         PetscReal dist = 0.0;
355 
356         if (PetscUnlikely(PetscBTLookup(excluded,j))) continue;
357         for (d=0;d<cdim;d++) dist += (wdist[j*cdim+d]-anchor[d])*(wdist[j*cdim+d]-anchor[d]);
358         if (dist > mdist) { mdist = dist; point2 = j; }
359       }
360 
361 
362       /* find the third point maximizing the triangle area */
363       point3 = point2;
364       if (cdim > 2) {
365         PetscReal a = 0.0;
366 
367         for (d=0;d<cdim;d++) a += (wdist[point1*cdim+d]-wdist[point2*cdim+d])*(wdist[point1*cdim+d]-wdist[point2*cdim+d]);
368         a = PetscSqrtReal(a);
369         mdist = -1.0;
370         for (j=fst;j<n;j++) {
371           PetscReal area,b = 0.0, c = 0.0,s;
372 
373           if (PetscUnlikely(PetscBTLookup(excluded,j))) continue;
374           for (d=0;d<cdim;d++) {
375             b += (wdist[point1*cdim+d]-wdist[j*cdim+d])*(wdist[point1*cdim+d]-wdist[j*cdim+d]);
376             c += (wdist[point2*cdim+d]-wdist[j*cdim+d])*(wdist[point2*cdim+d]-wdist[j*cdim+d]);
377           }
378           b = PetscSqrtReal(b);
379           c = PetscSqrtReal(c);
380           s = 0.5*(a+b+c);
381 
382           /* Heron's formula, area squared */
383           area = s*(s-a)*(s-b)*(s-c);
384           if (area > mdist) { mdist = area; point3 = j; }
385         }
386       }
387 
388       ierr = PetscBTSet(cornerp,shared[ns][point1]);CHKERRQ(ierr);
389       ierr = PetscBTSet(cornerp,shared[ns][point2]);CHKERRQ(ierr);
390       ierr = PetscBTSet(cornerp,shared[ns][point3]);CHKERRQ(ierr);
391 
392       /* all dofs having the same coordinates will be primal */
393       for (j=fst;j<n;j++) {
394         PetscBool same[3] = {PETSC_TRUE,PETSC_TRUE,PETSC_TRUE};
395 
396         if (PetscUnlikely(PetscBTLookup(excluded,j))) continue;
397         for (d=0;d<cdim;d++) {
398           same[0] = (PetscBool)(same[0] && (PetscAbsReal(wdist[j*cdim + d]-wdist[point1*cdim+d]) < PETSC_SMALL));
399           same[1] = (PetscBool)(same[1] && (PetscAbsReal(wdist[j*cdim + d]-wdist[point2*cdim+d]) < PETSC_SMALL));
400           same[2] = (PetscBool)(same[2] && (PetscAbsReal(wdist[j*cdim + d]-wdist[point3*cdim+d]) < PETSC_SMALL));
401         }
402         if (same[0] || same[1] || same[2]) {
403           ierr = PetscBTSet(cornerp,shared[ns][j]);CHKERRQ(ierr);
404         }
405       }
406     }
407     ierr = PetscBTDestroy(&excluded);CHKERRQ(ierr);
408     ierr = PetscFree(wdist);CHKERRQ(ierr);
409     ierr = ISLocalToGlobalMappingRestoreInfo(graph->l2gmap,&n_neigh,&neigh,&n_shared,&shared);CHKERRQ(ierr);
410   }
411 
412   /* check consistency of connected components among neighbouring subdomains -> it adapt them in case it is needed */
413   ierr = MPI_Comm_size(interface_comm,&size);CHKERRQ(ierr);
414   adapt_interface_reduced = PETSC_FALSE;
415   if (size > 1) {
416     PetscInt i;
417     PetscBool adapt_interface = cornerp ? PETSC_TRUE : PETSC_FALSE;
418     for (i=0;i<graph->n_subsets && !adapt_interface;i++) {
419       /* We are not sure that on a given subset of the local interface,
420          with two connected components, the latters be the same among sharing subdomains */
421       if (graph->subset_ncc[i] > 1) adapt_interface = PETSC_TRUE;
422     }
423     ierr = MPIU_Allreduce(&adapt_interface,&adapt_interface_reduced,1,MPIU_BOOL,MPI_LOR,interface_comm);CHKERRQ(ierr);
424   }
425 
426   if (graph->n_subsets && adapt_interface_reduced) {
427     PetscBT     subset_cc_adapt;
428     MPI_Request *send_requests,*recv_requests;
429     PetscInt    *send_buffer,*recv_buffer;
430     PetscInt    sum_requests,start_of_recv,start_of_send;
431     PetscInt    *cum_recv_counts;
432     PetscInt    *labels;
433     PetscInt    ncc,cum_queue,mss,mns,j,k,s;
434     PetscInt    **refine_buffer=NULL,*private_labels = NULL;
435     PetscBool   *subset_has_corn,*recv_buffer_bool,*send_buffer_bool;
436 
437     ierr = PetscCalloc1(graph->n_subsets,&subset_has_corn);CHKERRQ(ierr);
438     if (cornerp) {
439       for (i=0;i<graph->n_subsets;i++) {
440         for (j=0;j<graph->subset_size[i];j++) {
441           if (PetscBTLookup(cornerp,graph->subset_idxs[i][j])) {
442             subset_has_corn[i] = PETSC_TRUE;
443             break;
444           }
445         }
446       }
447     }
448     ierr = PetscMalloc1(graph->nvtxs,&labels);CHKERRQ(ierr);
449     ierr = PetscArrayzero(labels,graph->nvtxs);CHKERRQ(ierr);
450     for (i=0,k=0;i<graph->ncc;i++) {
451       PetscInt s = 1;
452       for (j=graph->cptr[i];j<graph->cptr[i+1];j++) {
453         if (cornerp && PetscBTLookup(cornerp,graph->queue[j])) {
454           labels[graph->queue[j]] = k+s;
455           s += 1;
456         } else {
457           labels[graph->queue[j]] = k;
458         }
459       }
460       k += s;
461     }
462 
463     /* allocate some space */
464     ierr = PetscMalloc1(graph->n_subsets+1,&cum_recv_counts);CHKERRQ(ierr);
465     ierr = PetscArrayzero(cum_recv_counts,graph->n_subsets+1);CHKERRQ(ierr);
466 
467     /* first count how many neighbours per connected component I will receive from */
468     cum_recv_counts[0] = 0;
469     for (i=0;i<graph->n_subsets;i++) cum_recv_counts[i+1] = cum_recv_counts[i]+graph->count[graph->subset_idxs[i][0]];
470     ierr = PetscMalloc1(graph->n_subsets,&send_buffer_bool);CHKERRQ(ierr);
471     ierr = PetscMalloc1(cum_recv_counts[graph->n_subsets],&recv_buffer_bool);CHKERRQ(ierr);
472     ierr = PetscMalloc2(cum_recv_counts[graph->n_subsets],&send_requests,cum_recv_counts[graph->n_subsets],&recv_requests);CHKERRQ(ierr);
473     for (i=0;i<cum_recv_counts[graph->n_subsets];i++) {
474       send_requests[i] = MPI_REQUEST_NULL;
475       recv_requests[i] = MPI_REQUEST_NULL;
476     }
477 
478     /* exchange with my neighbours the number of my connected components on the subset of interface */
479     sum_requests = 0;
480     for (i=0;i<graph->n_subsets;i++) {
481       send_buffer_bool[i] = (PetscBool)(graph->subset_ncc[i] > 1 || subset_has_corn[i]);
482     }
483     for (i=0;i<graph->n_subsets;i++) {
484       PetscMPIInt neigh,tag;
485       PetscInt    count,*neighs;
486 
487       count  = graph->count[graph->subset_idxs[i][0]];
488       neighs = graph->neighbours_set[graph->subset_idxs[i][0]];
489       ierr   = PetscMPIIntCast(2*graph->subset_ref_node[i],&tag);CHKERRQ(ierr);
490       for (k=0;k<count;k++) {
491 
492         ierr = PetscMPIIntCast(neighs[k],&neigh);CHKERRQ(ierr);
493         ierr = MPI_Isend(send_buffer_bool + i,           1,MPIU_BOOL,neigh,tag,interface_comm,&send_requests[sum_requests]);CHKERRQ(ierr);
494         ierr = MPI_Irecv(recv_buffer_bool + sum_requests,1,MPIU_BOOL,neigh,tag,interface_comm,&recv_requests[sum_requests]);CHKERRQ(ierr);
495         sum_requests++;
496       }
497     }
498     ierr = MPI_Waitall(sum_requests,recv_requests,MPI_STATUSES_IGNORE);CHKERRQ(ierr);
499     ierr = MPI_Waitall(sum_requests,send_requests,MPI_STATUSES_IGNORE);CHKERRQ(ierr);
500 
501     /* determine the subsets I have to adapt (those having more than 1 cc) */
502     ierr = PetscBTCreate(graph->n_subsets,&subset_cc_adapt);CHKERRQ(ierr);
503     ierr = PetscBTMemzero(graph->n_subsets,subset_cc_adapt);CHKERRQ(ierr);
504     for (i=0;i<graph->n_subsets;i++) {
505       if (graph->subset_ncc[i] > 1 || subset_has_corn[i]) {
506         ierr = PetscBTSet(subset_cc_adapt,i);CHKERRQ(ierr);
507         continue;
508       }
509       for (j=cum_recv_counts[i];j<cum_recv_counts[i+1];j++){
510          if (recv_buffer_bool[j]) {
511           ierr = PetscBTSet(subset_cc_adapt,i);CHKERRQ(ierr);
512           break;
513         }
514       }
515     }
516     ierr = PetscFree(send_buffer_bool);CHKERRQ(ierr);
517     ierr = PetscFree(recv_buffer_bool);CHKERRQ(ierr);
518     ierr = PetscFree(subset_has_corn);CHKERRQ(ierr);
519 
520     /* determine send/recv buffers sizes */
521     j = 0;
522     mss = 0;
523     for (i=0;i<graph->n_subsets;i++) {
524       if (PetscBTLookup(subset_cc_adapt,i)) {
525         j  += graph->subset_size[i];
526         mss = PetscMax(graph->subset_size[i],mss);
527       }
528     }
529     k = 0;
530     mns = 0;
531     for (i=0;i<graph->n_subsets;i++) {
532       if (PetscBTLookup(subset_cc_adapt,i)) {
533         k  += (cum_recv_counts[i+1]-cum_recv_counts[i])*graph->subset_size[i];
534         mns = PetscMax(cum_recv_counts[i+1]-cum_recv_counts[i],mns);
535       }
536     }
537     ierr = PetscMalloc2(j,&send_buffer,k,&recv_buffer);CHKERRQ(ierr);
538 
539     /* fill send buffer (order matters: subset_idxs ordered by global ordering) */
540     j = 0;
541     for (i=0;i<graph->n_subsets;i++)
542       if (PetscBTLookup(subset_cc_adapt,i))
543         for (k=0;k<graph->subset_size[i];k++)
544           send_buffer[j++] = labels[graph->subset_idxs[i][k]];
545 
546     /* now exchange the data */
547     start_of_recv = 0;
548     start_of_send = 0;
549     sum_requests  = 0;
550     for (i=0;i<graph->n_subsets;i++) {
551       if (PetscBTLookup(subset_cc_adapt,i)) {
552         PetscMPIInt neigh,tag;
553         PetscInt    size_of_send = graph->subset_size[i];
554 
555         j    = graph->subset_idxs[i][0];
556         ierr = PetscMPIIntCast(2*graph->subset_ref_node[i]+1,&tag);CHKERRQ(ierr);
557         for (k=0;k<graph->count[j];k++) {
558           ierr = PetscMPIIntCast(graph->neighbours_set[j][k],&neigh);CHKERRQ(ierr);
559           ierr = MPI_Isend(&send_buffer[start_of_send],size_of_send,MPIU_INT,neigh,tag,interface_comm,&send_requests[sum_requests]);CHKERRQ(ierr);
560           ierr = MPI_Irecv(&recv_buffer[start_of_recv],size_of_send,MPIU_INT,neigh,tag,interface_comm,&recv_requests[sum_requests]);CHKERRQ(ierr);
561           start_of_recv += size_of_send;
562           sum_requests++;
563         }
564         start_of_send += size_of_send;
565       }
566     }
567     ierr = MPI_Waitall(sum_requests,recv_requests,MPI_STATUSES_IGNORE);CHKERRQ(ierr);
568 
569     /* refine connected components */
570     start_of_recv = 0;
571     /* allocate some temporary space */
572     if (mss) {
573       ierr = PetscMalloc1(mss,&refine_buffer);CHKERRQ(ierr);
574       ierr = PetscMalloc2(mss*(mns+1),&refine_buffer[0],mss,&private_labels);CHKERRQ(ierr);
575     }
576     ncc = 0;
577     cum_queue = 0;
578     graph->cptr[0] = 0;
579     for (i=0;i<graph->n_subsets;i++) {
580       if (PetscBTLookup(subset_cc_adapt,i)) {
581         PetscInt subset_counter = 0;
582         PetscInt sharingprocs = cum_recv_counts[i+1]-cum_recv_counts[i]+1; /* count myself */
583         PetscInt buffer_size = graph->subset_size[i];
584 
585         /* compute pointers */
586         for (j=1;j<buffer_size;j++) refine_buffer[j] = refine_buffer[j-1] + sharingprocs;
587         /* analyze contributions from subdomains that share the i-th subset
588            The structure of refine_buffer is suitable to find intersections of ccs among sharingprocs.
589            supposing the current subset is shared by 3 processes and has dimension 5 with global dofs 0,1,2,3,4 (local 0,4,3,1,2)
590            sharing procs connected components:
591              neigh 0: [0 1 4], [2 3], labels [4,7]  (2 connected components)
592              neigh 1: [0 1], [2 3 4], labels [3 2]  (2 connected components)
593              neigh 2: [0 4], [1], [2 3], labels [1 5 6] (3 connected components)
594            refine_buffer will be filled as:
595              [ 4, 3, 1;
596                4, 2, 1;
597                7, 2, 6;
598                4, 3, 5;
599                7, 2, 6; ];
600            The connected components in local ordering are [0], [1], [2 3], [4] */
601         /* fill temp_buffer */
602         for (k=0;k<buffer_size;k++) refine_buffer[k][0] = labels[graph->subset_idxs[i][k]];
603         for (j=0;j<sharingprocs-1;j++) {
604           for (k=0;k<buffer_size;k++) refine_buffer[k][j+1] = recv_buffer[start_of_recv+k];
605           start_of_recv += buffer_size;
606         }
607         ierr = PetscArrayzero(private_labels,buffer_size);CHKERRQ(ierr);
608         for (j=0;j<buffer_size;j++) {
609           if (!private_labels[j]) { /* found a new cc  */
610             PetscBool same_set;
611 
612             graph->cptr[ncc] = cum_queue;
613             ncc++;
614             subset_counter++;
615             private_labels[j] = subset_counter;
616             graph->queue[cum_queue++] = graph->subset_idxs[i][j];
617             for (k=j+1;k<buffer_size;k++) { /* check for other nodes in new cc */
618               same_set = PETSC_TRUE;
619               for (s=0;s<sharingprocs;s++) {
620                 if (refine_buffer[j][s] != refine_buffer[k][s]) {
621                   same_set = PETSC_FALSE;
622                   break;
623                 }
624               }
625               if (same_set) {
626                 private_labels[k] = subset_counter;
627                 graph->queue[cum_queue++] = graph->subset_idxs[i][k];
628               }
629             }
630           }
631         }
632         graph->cptr[ncc]     = cum_queue;
633         graph->subset_ncc[i] = subset_counter;
634         graph->queue_sorted  = PETSC_FALSE;
635       } else { /* this subset does not need to be adapted */
636         ierr = PetscArraycpy(graph->queue+cum_queue,graph->subset_idxs[i],graph->subset_size[i]);CHKERRQ(ierr);
637         ncc++;
638         cum_queue += graph->subset_size[i];
639         graph->cptr[ncc] = cum_queue;
640       }
641     }
642     graph->cptr[ncc] = cum_queue;
643     graph->ncc       = ncc;
644     if (mss) {
645       ierr = PetscFree2(refine_buffer[0],private_labels);CHKERRQ(ierr);
646       ierr = PetscFree(refine_buffer);CHKERRQ(ierr);
647     }
648     ierr = PetscFree(labels);CHKERRQ(ierr);
649     ierr = MPI_Waitall(sum_requests,send_requests,MPI_STATUSES_IGNORE);CHKERRQ(ierr);
650     ierr = PetscFree2(send_requests,recv_requests);CHKERRQ(ierr);
651     ierr = PetscFree2(send_buffer,recv_buffer);CHKERRQ(ierr);
652     ierr = PetscFree(cum_recv_counts);CHKERRQ(ierr);
653     ierr = PetscBTDestroy(&subset_cc_adapt);CHKERRQ(ierr);
654   }
655   ierr = PetscBTDestroy(&cornerp);CHKERRQ(ierr);
656 
657   /* Determine if we are in 2D or 3D */
658   if (!graph->twodimset) {
659     PetscBool twodim = PETSC_TRUE;
660     for (i=0;i<graph->ncc;i++) {
661       PetscInt repdof = graph->queue[graph->cptr[i]];
662       PetscInt ccsize = graph->cptr[i+1]-graph->cptr[i];
663       if (graph->count[repdof] > 1 && ccsize > graph->custom_minimal_size) {
664         twodim = PETSC_FALSE;
665         break;
666       }
667     }
668     ierr = MPIU_Allreduce(&twodim,&graph->twodim,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)graph->l2gmap));CHKERRQ(ierr);
669     graph->twodimset = PETSC_TRUE;
670   }
671   PetscFunctionReturn(0);
672 }
673 
674 
PCBDDCGraphComputeCC_Private(PCBDDCGraph graph,PetscInt pid,PetscInt * queue_tip,PetscInt n_prev,PetscInt * n_added)675 PETSC_STATIC_INLINE PetscErrorCode PCBDDCGraphComputeCC_Private(PCBDDCGraph graph,PetscInt pid,PetscInt* queue_tip,PetscInt n_prev,PetscInt* n_added)
676 {
677   PetscInt       i,j,n;
678   PetscInt       *xadj = graph->xadj,*adjncy = graph->adjncy;
679   PetscBT        touched = graph->touched;
680   PetscBool      havecsr = (PetscBool)(!!xadj);
681   PetscBool      havesubs = (PetscBool)(!!graph->n_local_subs);
682   PetscErrorCode ierr;
683 
684   PetscFunctionBegin;
685   n = 0;
686   if (havecsr && !havesubs) {
687     for (i=-n_prev;i<0;i++) {
688       PetscInt start_dof = queue_tip[i];
689       /* we assume that if a dof has a size 1 adjacency list and the corresponding entry is negative, it is connected to all dofs */
690       if (xadj[start_dof+1]-xadj[start_dof] == 1 && adjncy[xadj[start_dof]] < 0) {
691         for (j=0;j<graph->subset_size[pid-1];j++) { /* pid \in [1,graph->n_subsets] */
692           PetscInt dof = graph->subset_idxs[pid-1][j];
693           if (!PetscBTLookup(touched,dof) && graph->subset[dof] == pid) {
694             ierr = PetscBTSet(touched,dof);CHKERRQ(ierr);
695             queue_tip[n] = dof;
696             n++;
697           }
698         }
699       } else {
700         for (j=xadj[start_dof];j<xadj[start_dof+1];j++) {
701           PetscInt dof = adjncy[j];
702           if (!PetscBTLookup(touched,dof) && graph->subset[dof] == pid) {
703             ierr = PetscBTSet(touched,dof);CHKERRQ(ierr);
704             queue_tip[n] = dof;
705             n++;
706           }
707         }
708       }
709     }
710   } else if (havecsr && havesubs) {
711     PetscInt sid = graph->local_subs[queue_tip[-n_prev]];
712     for (i=-n_prev;i<0;i++) {
713       PetscInt start_dof = queue_tip[i];
714       /* we assume that if a dof has a size 1 adjacency list and the corresponding entry is negative, it is connected to all dofs belonging to the local sub */
715       if (xadj[start_dof+1]-xadj[start_dof] == 1 && adjncy[xadj[start_dof]] < 0) {
716         for (j=0;j<graph->subset_size[pid-1];j++) { /* pid \in [1,graph->n_subsets] */
717           PetscInt dof = graph->subset_idxs[pid-1][j];
718           if (!PetscBTLookup(touched,dof) && graph->subset[dof] == pid && graph->local_subs[dof] == sid) {
719             ierr = PetscBTSet(touched,dof);CHKERRQ(ierr);
720             queue_tip[n] = dof;
721             n++;
722           }
723         }
724       } else {
725         for (j=xadj[start_dof];j<xadj[start_dof+1];j++) {
726           PetscInt dof = adjncy[j];
727           if (!PetscBTLookup(touched,dof) && graph->subset[dof] == pid && graph->local_subs[dof] == sid) {
728             ierr = PetscBTSet(touched,dof);CHKERRQ(ierr);
729             queue_tip[n] = dof;
730             n++;
731           }
732         }
733       }
734     }
735   } else if (havesubs) { /* sub info only */
736     PetscInt sid = graph->local_subs[queue_tip[-n_prev]];
737     for (j=0;j<graph->subset_size[pid-1];j++) { /* pid \in [1,graph->n_subsets] */
738       PetscInt dof = graph->subset_idxs[pid-1][j];
739       if (!PetscBTLookup(touched,dof) && graph->subset[dof] == pid && graph->local_subs[dof] == sid) {
740         ierr = PetscBTSet(touched,dof);CHKERRQ(ierr);
741         queue_tip[n] = dof;
742         n++;
743       }
744     }
745   } else {
746     for (j=0;j<graph->subset_size[pid-1];j++) { /* pid \in [1,graph->n_subsets] */
747       PetscInt dof = graph->subset_idxs[pid-1][j];
748       if (!PetscBTLookup(touched,dof) && graph->subset[dof] == pid) {
749         ierr = PetscBTSet(touched,dof);CHKERRQ(ierr);
750         queue_tip[n] = dof;
751         n++;
752       }
753     }
754   }
755   *n_added = n;
756   PetscFunctionReturn(0);
757 }
758 
PCBDDCGraphComputeConnectedComponentsLocal(PCBDDCGraph graph)759 PetscErrorCode PCBDDCGraphComputeConnectedComponentsLocal(PCBDDCGraph graph)
760 {
761   PetscInt       ncc,cum_queue,n;
762   PetscMPIInt    commsize;
763   PetscErrorCode ierr;
764 
765   PetscFunctionBegin;
766   if (!graph->setupcalled) SETERRQ(PetscObjectComm((PetscObject)graph->l2gmap),PETSC_ERR_ORDER,"PCBDDCGraphSetUp should be called first");
767   /* quiet return if there isn't any local info */
768   if (!graph->xadj && !graph->n_local_subs) {
769     PetscFunctionReturn(0);
770   }
771 
772   /* reset any previous search of connected components */
773   ierr = PetscBTMemzero(graph->nvtxs,graph->touched);CHKERRQ(ierr);
774   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)graph->l2gmap),&commsize);CHKERRQ(ierr);
775   if (commsize > graph->commsizelimit) {
776     PetscInt i;
777     for (i=0;i<graph->nvtxs;i++) {
778       if (graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK || !graph->count[i]) {
779         ierr = PetscBTSet(graph->touched,i);CHKERRQ(ierr);
780       }
781     }
782   }
783 
784   /* begin search for connected components */
785   cum_queue = 0;
786   ncc = 0;
787   for (n=0;n<graph->n_subsets;n++) {
788     PetscInt pid = n+1;  /* partition labeled by 0 is discarded */
789     PetscInt found = 0,prev = 0,first = 0,ncc_pid = 0;
790     while (found != graph->subset_size[n]) {
791       PetscInt added = 0;
792       if (!prev) { /* search for new starting dof */
793         while (PetscBTLookup(graph->touched,graph->subset_idxs[n][first])) first++;
794         ierr = PetscBTSet(graph->touched,graph->subset_idxs[n][first]);CHKERRQ(ierr);
795         graph->queue[cum_queue] = graph->subset_idxs[n][first];
796         graph->cptr[ncc] = cum_queue;
797         prev = 1;
798         cum_queue++;
799         found++;
800         ncc_pid++;
801         ncc++;
802       }
803       ierr = PCBDDCGraphComputeCC_Private(graph,pid,graph->queue + cum_queue,prev,&added);CHKERRQ(ierr);
804       if (!added) {
805         graph->subset_ncc[n] = ncc_pid;
806         graph->cptr[ncc] = cum_queue;
807       }
808       prev = added;
809       found += added;
810       cum_queue += added;
811       if (added && found == graph->subset_size[n]) {
812         graph->subset_ncc[n] = ncc_pid;
813         graph->cptr[ncc] = cum_queue;
814       }
815     }
816   }
817   graph->ncc = ncc;
818   graph->queue_sorted = PETSC_FALSE;
819   PetscFunctionReturn(0);
820 }
821 
PCBDDCGraphSetUp(PCBDDCGraph graph,PetscInt custom_minimal_size,IS neumann_is,IS dirichlet_is,PetscInt n_ISForDofs,IS ISForDofs[],IS custom_primal_vertices)822 PetscErrorCode PCBDDCGraphSetUp(PCBDDCGraph graph, PetscInt custom_minimal_size, IS neumann_is, IS dirichlet_is, PetscInt n_ISForDofs, IS ISForDofs[], IS custom_primal_vertices)
823 {
824   IS             subset,subset_n;
825   MPI_Comm       comm;
826   const PetscInt *is_indices;
827   PetscInt       n_neigh,*neigh,*n_shared,**shared,*queue_global;
828   PetscInt       i,j,k,s,total_counts,nodes_touched,is_size;
829   PetscMPIInt    commsize;
830   PetscBool      same_set,mirrors_found;
831   PetscErrorCode ierr;
832 
833   PetscFunctionBegin;
834   PetscValidLogicalCollectiveInt(graph->l2gmap,custom_minimal_size,2);
835   if (neumann_is) {
836     PetscValidHeaderSpecific(neumann_is,IS_CLASSID,3);
837     PetscCheckSameComm(graph->l2gmap,1,neumann_is,3);
838   }
839   graph->has_dirichlet = PETSC_FALSE;
840   if (dirichlet_is) {
841     PetscValidHeaderSpecific(dirichlet_is,IS_CLASSID,4);
842     PetscCheckSameComm(graph->l2gmap,1,dirichlet_is,4);
843     graph->has_dirichlet = PETSC_TRUE;
844   }
845   PetscValidLogicalCollectiveInt(graph->l2gmap,n_ISForDofs,5);
846   for (i=0;i<n_ISForDofs;i++) {
847     PetscValidHeaderSpecific(ISForDofs[i],IS_CLASSID,6);
848     PetscCheckSameComm(graph->l2gmap,1,ISForDofs[i],6);
849   }
850   if (custom_primal_vertices) {
851     PetscValidHeaderSpecific(custom_primal_vertices,IS_CLASSID,7);
852     PetscCheckSameComm(graph->l2gmap,1,custom_primal_vertices,7);
853   }
854   ierr = PetscObjectGetComm((PetscObject)(graph->l2gmap),&comm);CHKERRQ(ierr);
855   ierr = MPI_Comm_size(comm,&commsize);CHKERRQ(ierr);
856 
857   /* custom_minimal_size */
858   graph->custom_minimal_size = custom_minimal_size;
859   /* get info l2gmap and allocate work vectors  */
860   ierr = ISLocalToGlobalMappingGetInfo(graph->l2gmap,&n_neigh,&neigh,&n_shared,&shared);CHKERRQ(ierr);
861   /* check if we have any local periodic nodes (periodic BCs) */
862   mirrors_found = PETSC_FALSE;
863   if (graph->nvtxs && n_neigh) {
864     for (i=0; i<n_shared[0]; i++) graph->count[shared[0][i]] += 1;
865     for (i=0; i<n_shared[0]; i++) {
866       if (graph->count[shared[0][i]] > 1) {
867         mirrors_found = PETSC_TRUE;
868         break;
869       }
870     }
871   }
872   /* compute local mirrors (if any) */
873   if (mirrors_found) {
874     IS       to,from;
875     PetscInt *local_indices,*global_indices;
876 
877     ierr = ISCreateStride(PETSC_COMM_SELF,graph->nvtxs,0,1,&to);CHKERRQ(ierr);
878     ierr = ISLocalToGlobalMappingApplyIS(graph->l2gmap,to,&from);CHKERRQ(ierr);
879     /* get arrays of local and global indices */
880     ierr = PetscMalloc1(graph->nvtxs,&local_indices);CHKERRQ(ierr);
881     ierr = ISGetIndices(to,(const PetscInt**)&is_indices);CHKERRQ(ierr);
882     ierr = PetscArraycpy(local_indices,is_indices,graph->nvtxs);CHKERRQ(ierr);
883     ierr = ISRestoreIndices(to,(const PetscInt**)&is_indices);CHKERRQ(ierr);
884     ierr = PetscMalloc1(graph->nvtxs,&global_indices);CHKERRQ(ierr);
885     ierr = ISGetIndices(from,(const PetscInt**)&is_indices);CHKERRQ(ierr);
886     ierr = PetscArraycpy(global_indices,is_indices,graph->nvtxs);CHKERRQ(ierr);
887     ierr = ISRestoreIndices(from,(const PetscInt**)&is_indices);CHKERRQ(ierr);
888     /* allocate space for mirrors */
889     ierr = PetscMalloc2(graph->nvtxs,&graph->mirrors,graph->nvtxs,&graph->mirrors_set);CHKERRQ(ierr);
890     ierr = PetscArrayzero(graph->mirrors,graph->nvtxs);CHKERRQ(ierr);
891     graph->mirrors_set[0] = NULL;
892 
893     k=0;
894     for (i=0;i<n_shared[0];i++) {
895       j=shared[0][i];
896       if (graph->count[j] > 1) {
897         graph->mirrors[j]++;
898         k++;
899       }
900     }
901     /* allocate space for set of mirrors */
902     ierr = PetscMalloc1(k,&graph->mirrors_set[0]);CHKERRQ(ierr);
903     for (i=1;i<graph->nvtxs;i++)
904       graph->mirrors_set[i]=graph->mirrors_set[i-1]+graph->mirrors[i-1];
905 
906     /* fill arrays */
907     ierr = PetscArrayzero(graph->mirrors,graph->nvtxs);CHKERRQ(ierr);
908     for (j=0;j<n_shared[0];j++) {
909       i=shared[0][j];
910       if (graph->count[i] > 1)
911         graph->mirrors_set[i][graph->mirrors[i]++]=global_indices[i];
912     }
913     ierr = PetscSortIntWithArray(graph->nvtxs,global_indices,local_indices);CHKERRQ(ierr);
914     for (i=0;i<graph->nvtxs;i++) {
915       if (graph->mirrors[i] > 0) {
916         ierr = PetscFindInt(graph->mirrors_set[i][0],graph->nvtxs,global_indices,&k);CHKERRQ(ierr);
917         j = global_indices[k];
918         while (k > 0 && global_indices[k-1] == j) k--;
919         for (j=0;j<graph->mirrors[i];j++) {
920           graph->mirrors_set[i][j]=local_indices[k+j];
921         }
922         ierr = PetscSortInt(graph->mirrors[i],graph->mirrors_set[i]);CHKERRQ(ierr);
923       }
924     }
925     ierr = PetscFree(local_indices);CHKERRQ(ierr);
926     ierr = PetscFree(global_indices);CHKERRQ(ierr);
927     ierr = ISDestroy(&to);CHKERRQ(ierr);
928     ierr = ISDestroy(&from);CHKERRQ(ierr);
929   }
930   ierr = PetscArrayzero(graph->count,graph->nvtxs);CHKERRQ(ierr);
931 
932   /* Count total number of neigh per node */
933   k = 0;
934   for (i=1;i<n_neigh;i++) {
935     k += n_shared[i];
936     for (j=0;j<n_shared[i];j++) {
937       graph->count[shared[i][j]] += 1;
938     }
939   }
940   /* Allocate space for storing the set of neighbours for each node */
941   if (graph->nvtxs) {
942     ierr = PetscMalloc1(k,&graph->neighbours_set[0]);CHKERRQ(ierr);
943   }
944   for (i=1;i<graph->nvtxs;i++) { /* dont count myself */
945     graph->neighbours_set[i]=graph->neighbours_set[i-1]+graph->count[i-1];
946   }
947   /* Get information for sharing subdomains */
948   ierr = PetscArrayzero(graph->count,graph->nvtxs);CHKERRQ(ierr);
949   for (i=1;i<n_neigh;i++) { /* dont count myself */
950     s = n_shared[i];
951     for (j=0;j<s;j++) {
952       k = shared[i][j];
953       graph->neighbours_set[k][graph->count[k]] = neigh[i];
954       graph->count[k] += 1;
955     }
956   }
957   /* sort set of sharing subdomains */
958   for (i=0;i<graph->nvtxs;i++) {
959     ierr = PetscSortRemoveDupsInt(&graph->count[i],graph->neighbours_set[i]);CHKERRQ(ierr);
960   }
961   /* free memory allocated by ISLocalToGlobalMappingGetInfo */
962   ierr = ISLocalToGlobalMappingRestoreInfo(graph->l2gmap,&n_neigh,&neigh,&n_shared,&shared);CHKERRQ(ierr);
963 
964   /*
965      Get info for dofs splitting
966      User can specify just a subset; an additional field is considered as a complementary field
967   */
968   for (i=0,k=0;i<n_ISForDofs;i++) {
969     PetscInt bs;
970 
971     ierr = ISGetBlockSize(ISForDofs[i],&bs);CHKERRQ(ierr);
972     k   += bs;
973   }
974   for (i=0;i<graph->nvtxs;i++) graph->which_dof[i] = k; /* by default a dof belongs to the complement set */
975   for (i=0,k=0;i<n_ISForDofs;i++) {
976     PetscInt bs;
977 
978     ierr = ISGetLocalSize(ISForDofs[i],&is_size);CHKERRQ(ierr);
979     ierr = ISGetBlockSize(ISForDofs[i],&bs);CHKERRQ(ierr);
980     ierr = ISGetIndices(ISForDofs[i],(const PetscInt**)&is_indices);CHKERRQ(ierr);
981     for (j=0;j<is_size/bs;j++) {
982       PetscInt b;
983 
984       for (b=0;b<bs;b++) {
985         PetscInt jj = bs*j + b;
986 
987         if (is_indices[jj] > -1 && is_indices[jj] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */
988           graph->which_dof[is_indices[jj]] = k+b;
989         }
990       }
991     }
992     ierr = ISRestoreIndices(ISForDofs[i],(const PetscInt**)&is_indices);CHKERRQ(ierr);
993     k   += bs;
994   }
995 
996   /* Take into account Neumann nodes */
997   if (neumann_is) {
998     ierr = ISGetLocalSize(neumann_is,&is_size);CHKERRQ(ierr);
999     ierr = ISGetIndices(neumann_is,(const PetscInt**)&is_indices);CHKERRQ(ierr);
1000     for (i=0;i<is_size;i++) {
1001       if (is_indices[i] > -1 && is_indices[i] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */
1002         graph->special_dof[is_indices[i]] = PCBDDCGRAPH_NEUMANN_MARK;
1003       }
1004     }
1005     ierr = ISRestoreIndices(neumann_is,(const PetscInt**)&is_indices);CHKERRQ(ierr);
1006   }
1007   /* Take into account Dirichlet nodes (they overwrite any neumann boundary mark previously set) */
1008   if (dirichlet_is) {
1009     ierr = ISGetLocalSize(dirichlet_is,&is_size);CHKERRQ(ierr);
1010     ierr = ISGetIndices(dirichlet_is,(const PetscInt**)&is_indices);CHKERRQ(ierr);
1011     for (i=0;i<is_size;i++){
1012       if (is_indices[i] > -1 && is_indices[i] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */
1013         if (commsize > graph->commsizelimit) { /* dirichlet nodes treated as internal */
1014           ierr = PetscBTSet(graph->touched,is_indices[i]);CHKERRQ(ierr);
1015           graph->subset[is_indices[i]] = 0;
1016         }
1017         graph->special_dof[is_indices[i]] = PCBDDCGRAPH_DIRICHLET_MARK;
1018       }
1019     }
1020     ierr = ISRestoreIndices(dirichlet_is,(const PetscInt**)&is_indices);CHKERRQ(ierr);
1021   }
1022   /* mark local periodic nodes (if any) and adapt CSR graph (if any) */
1023   if (graph->mirrors) {
1024     for (i=0;i<graph->nvtxs;i++)
1025       if (graph->mirrors[i])
1026         graph->special_dof[i] = PCBDDCGRAPH_LOCAL_PERIODIC_MARK;
1027 
1028     if (graph->xadj) {
1029       PetscInt *new_xadj,*new_adjncy;
1030       /* sort CSR graph */
1031       for (i=0;i<graph->nvtxs;i++)
1032         ierr = PetscSortInt(graph->xadj[i+1]-graph->xadj[i],&graph->adjncy[graph->xadj[i]]);CHKERRQ(ierr);
1033 
1034       /* adapt local CSR graph in case of local periodicity */
1035       k = 0;
1036       for (i=0;i<graph->nvtxs;i++)
1037         for (j=graph->xadj[i];j<graph->xadj[i+1];j++)
1038           k += graph->mirrors[graph->adjncy[j]];
1039 
1040       ierr = PetscMalloc1(graph->nvtxs+1,&new_xadj);CHKERRQ(ierr);
1041       ierr = PetscMalloc1(k+graph->xadj[graph->nvtxs],&new_adjncy);CHKERRQ(ierr);
1042       new_xadj[0] = 0;
1043       for (i=0;i<graph->nvtxs;i++) {
1044         k = graph->xadj[i+1]-graph->xadj[i];
1045         ierr = PetscArraycpy(&new_adjncy[new_xadj[i]],&graph->adjncy[graph->xadj[i]],k);CHKERRQ(ierr);
1046         new_xadj[i+1] = new_xadj[i]+k;
1047         for (j=graph->xadj[i];j<graph->xadj[i+1];j++) {
1048           k = graph->mirrors[graph->adjncy[j]];
1049           ierr = PetscArraycpy(&new_adjncy[new_xadj[i+1]],graph->mirrors_set[graph->adjncy[j]],k);CHKERRQ(ierr);
1050           new_xadj[i+1] += k;
1051         }
1052         k = new_xadj[i+1]-new_xadj[i];
1053         ierr = PetscSortRemoveDupsInt(&k,&new_adjncy[new_xadj[i]]);CHKERRQ(ierr);
1054         new_xadj[i+1] = new_xadj[i]+k;
1055       }
1056       /* set new CSR into graph */
1057       ierr = PetscFree(graph->xadj);CHKERRQ(ierr);
1058       ierr = PetscFree(graph->adjncy);CHKERRQ(ierr);
1059       graph->xadj = new_xadj;
1060       graph->adjncy = new_adjncy;
1061     }
1062   }
1063 
1064   /* mark special nodes (if any) -> each will become a single node equivalence class */
1065   if (custom_primal_vertices) {
1066     ierr = ISGetLocalSize(custom_primal_vertices,&is_size);CHKERRQ(ierr);
1067     ierr = ISGetIndices(custom_primal_vertices,(const PetscInt**)&is_indices);CHKERRQ(ierr);
1068     for (i=0,j=0;i<is_size;i++){
1069       if (is_indices[i] > -1 && is_indices[i] < graph->nvtxs  && graph->special_dof[is_indices[i]] != PCBDDCGRAPH_DIRICHLET_MARK) { /* out of bounds indices (if any) are skipped */
1070         graph->special_dof[is_indices[i]] = PCBDDCGRAPH_SPECIAL_MARK-j;
1071         j++;
1072       }
1073     }
1074     ierr = ISRestoreIndices(custom_primal_vertices,(const PetscInt**)&is_indices);CHKERRQ(ierr);
1075   }
1076 
1077   /* mark interior nodes (if commsize > graph->commsizelimit) as touched and belonging to partition number 0 */
1078   if (commsize > graph->commsizelimit) {
1079     for (i=0;i<graph->nvtxs;i++) {
1080       if (!graph->count[i]) {
1081         ierr = PetscBTSet(graph->touched,i);CHKERRQ(ierr);
1082         graph->subset[i] = 0;
1083       }
1084     }
1085   }
1086 
1087   /* init graph structure and compute default subsets */
1088   nodes_touched = 0;
1089   for (i=0;i<graph->nvtxs;i++) {
1090     if (PetscBTLookup(graph->touched,i)) {
1091       nodes_touched++;
1092     }
1093   }
1094   i = 0;
1095   graph->ncc = 0;
1096   total_counts = 0;
1097 
1098   /* allocated space for queues */
1099   if (commsize == graph->commsizelimit) {
1100     ierr = PetscMalloc2(graph->nvtxs+1,&graph->cptr,graph->nvtxs,&graph->queue);CHKERRQ(ierr);
1101   } else {
1102     PetscInt nused = graph->nvtxs - nodes_touched;
1103     ierr = PetscMalloc2(nused+1,&graph->cptr,nused,&graph->queue);CHKERRQ(ierr);
1104   }
1105 
1106   while (nodes_touched<graph->nvtxs) {
1107     /*  find first untouched node in local ordering */
1108     while (PetscBTLookup(graph->touched,i)) i++;
1109     ierr = PetscBTSet(graph->touched,i);CHKERRQ(ierr);
1110     graph->subset[i] = graph->ncc+1;
1111     graph->cptr[graph->ncc] = total_counts;
1112     graph->queue[total_counts] = i;
1113     total_counts++;
1114     nodes_touched++;
1115     /* now find all other nodes having the same set of sharing subdomains */
1116     for (j=i+1;j<graph->nvtxs;j++) {
1117       /* check for same number of sharing subdomains, dof number and same special mark */
1118       if (!PetscBTLookup(graph->touched,j) && graph->count[i] == graph->count[j] && graph->which_dof[i] == graph->which_dof[j] && graph->special_dof[i] == graph->special_dof[j]) {
1119         /* check for same set of sharing subdomains */
1120         same_set = PETSC_TRUE;
1121         for (k=0;k<graph->count[j];k++){
1122           if (graph->neighbours_set[i][k] != graph->neighbours_set[j][k]) {
1123             same_set = PETSC_FALSE;
1124           }
1125         }
1126         /* I have found a friend of mine */
1127         if (same_set) {
1128           ierr = PetscBTSet(graph->touched,j);CHKERRQ(ierr);
1129           graph->subset[j] = graph->ncc+1;
1130           nodes_touched++;
1131           graph->queue[total_counts] = j;
1132           total_counts++;
1133         }
1134       }
1135     }
1136     graph->ncc++;
1137   }
1138   /* set default number of subsets (at this point no info on csr and/or local_subs has been taken into account, so n_subsets = ncc */
1139   graph->n_subsets = graph->ncc;
1140   ierr = PetscMalloc1(graph->n_subsets,&graph->subset_ncc);CHKERRQ(ierr);
1141   for (i=0;i<graph->n_subsets;i++) {
1142     graph->subset_ncc[i] = 1;
1143   }
1144   /* final pointer */
1145   graph->cptr[graph->ncc] = total_counts;
1146 
1147   /* For consistency reasons (among neighbours), I need to sort (by global ordering) each connected component */
1148   /* Get a reference node (min index in global ordering) for each subset for tagging messages */
1149   ierr = PetscMalloc1(graph->ncc,&graph->subset_ref_node);CHKERRQ(ierr);
1150   ierr = PetscMalloc1(graph->cptr[graph->ncc],&queue_global);CHKERRQ(ierr);
1151   ierr = ISLocalToGlobalMappingApply(graph->l2gmap,graph->cptr[graph->ncc],graph->queue,queue_global);CHKERRQ(ierr);
1152   for (j=0;j<graph->ncc;j++) {
1153     ierr = PetscSortIntWithArray(graph->cptr[j+1]-graph->cptr[j],&queue_global[graph->cptr[j]],&graph->queue[graph->cptr[j]]);CHKERRQ(ierr);
1154     graph->subset_ref_node[j] = graph->queue[graph->cptr[j]];
1155   }
1156   ierr = PetscFree(queue_global);CHKERRQ(ierr);
1157   graph->queue_sorted = PETSC_TRUE;
1158 
1159   /* save information on subsets (needed when analyzing the connected components) */
1160   if (graph->ncc) {
1161     ierr = PetscMalloc2(graph->ncc,&graph->subset_size,graph->ncc,&graph->subset_idxs);CHKERRQ(ierr);
1162     ierr = PetscMalloc1(graph->cptr[graph->ncc],&graph->subset_idxs[0]);CHKERRQ(ierr);
1163     ierr = PetscArrayzero(graph->subset_idxs[0],graph->cptr[graph->ncc]);CHKERRQ(ierr);
1164     for (j=1;j<graph->ncc;j++) {
1165       graph->subset_size[j-1] = graph->cptr[j] - graph->cptr[j-1];
1166       graph->subset_idxs[j] = graph->subset_idxs[j-1] + graph->subset_size[j-1];
1167     }
1168     graph->subset_size[graph->ncc-1] = graph->cptr[graph->ncc] - graph->cptr[graph->ncc-1];
1169     ierr = PetscArraycpy(graph->subset_idxs[0],graph->queue,graph->cptr[graph->ncc]);CHKERRQ(ierr);
1170   }
1171 
1172   /* renumber reference nodes */
1173   ierr = ISCreateGeneral(PetscObjectComm((PetscObject)(graph->l2gmap)),graph->ncc,graph->subset_ref_node,PETSC_COPY_VALUES,&subset_n);CHKERRQ(ierr);
1174   ierr = ISLocalToGlobalMappingApplyIS(graph->l2gmap,subset_n,&subset);CHKERRQ(ierr);
1175   ierr = ISDestroy(&subset_n);CHKERRQ(ierr);
1176   ierr = ISRenumber(subset,NULL,NULL,&subset_n);CHKERRQ(ierr);
1177   ierr = ISDestroy(&subset);CHKERRQ(ierr);
1178   ierr = ISGetLocalSize(subset_n,&k);CHKERRQ(ierr);
1179   if (k != graph->ncc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Invalid size of new subset! %D != %D",k,graph->ncc);
1180   ierr = ISGetIndices(subset_n,&is_indices);CHKERRQ(ierr);
1181   ierr = PetscArraycpy(graph->subset_ref_node,is_indices,graph->ncc);CHKERRQ(ierr);
1182   ierr = ISRestoreIndices(subset_n,&is_indices);CHKERRQ(ierr);
1183   ierr = ISDestroy(&subset_n);CHKERRQ(ierr);
1184 
1185   /* free workspace */
1186   graph->setupcalled = PETSC_TRUE;
1187   PetscFunctionReturn(0);
1188 }
1189 
PCBDDCGraphResetCoords(PCBDDCGraph graph)1190 PetscErrorCode PCBDDCGraphResetCoords(PCBDDCGraph graph)
1191 {
1192   PetscErrorCode ierr;
1193 
1194   PetscFunctionBegin;
1195   if (!graph) PetscFunctionReturn(0);
1196   ierr = PetscFree(graph->coords);CHKERRQ(ierr);
1197   graph->cdim  = 0;
1198   graph->cnloc = 0;
1199   graph->cloc  = PETSC_FALSE;
1200   PetscFunctionReturn(0);
1201 }
1202 
PCBDDCGraphResetCSR(PCBDDCGraph graph)1203 PetscErrorCode PCBDDCGraphResetCSR(PCBDDCGraph graph)
1204 {
1205   PetscErrorCode ierr;
1206 
1207   PetscFunctionBegin;
1208   if (!graph) PetscFunctionReturn(0);
1209   if (graph->freecsr) {
1210     ierr = PetscFree(graph->xadj);CHKERRQ(ierr);
1211     ierr = PetscFree(graph->adjncy);CHKERRQ(ierr);
1212   } else {
1213     graph->xadj = NULL;
1214     graph->adjncy = NULL;
1215   }
1216   graph->freecsr = PETSC_FALSE;
1217   graph->nvtxs_csr = 0;
1218   PetscFunctionReturn(0);
1219 }
1220 
PCBDDCGraphReset(PCBDDCGraph graph)1221 PetscErrorCode PCBDDCGraphReset(PCBDDCGraph graph)
1222 {
1223   PetscErrorCode ierr;
1224 
1225   PetscFunctionBegin;
1226   if (!graph) PetscFunctionReturn(0);
1227   ierr = ISLocalToGlobalMappingDestroy(&graph->l2gmap);CHKERRQ(ierr);
1228   ierr = PetscFree(graph->subset_ncc);CHKERRQ(ierr);
1229   ierr = PetscFree(graph->subset_ref_node);CHKERRQ(ierr);
1230   if (graph->nvtxs) {
1231     ierr = PetscFree(graph->neighbours_set[0]);CHKERRQ(ierr);
1232   }
1233   ierr = PetscBTDestroy(&graph->touched);CHKERRQ(ierr);
1234   ierr = PetscFree5(graph->count,
1235                     graph->neighbours_set,
1236                     graph->subset,
1237                     graph->which_dof,
1238                     graph->special_dof);CHKERRQ(ierr);
1239   ierr = PetscFree2(graph->cptr,graph->queue);CHKERRQ(ierr);
1240   if (graph->mirrors) {
1241     ierr = PetscFree(graph->mirrors_set[0]);CHKERRQ(ierr);
1242   }
1243   ierr = PetscFree2(graph->mirrors,graph->mirrors_set);CHKERRQ(ierr);
1244   if (graph->subset_idxs) {
1245     ierr = PetscFree(graph->subset_idxs[0]);CHKERRQ(ierr);
1246   }
1247   ierr = PetscFree2(graph->subset_size,graph->subset_idxs);CHKERRQ(ierr);
1248   ierr = ISDestroy(&graph->dirdofs);CHKERRQ(ierr);
1249   ierr = ISDestroy(&graph->dirdofsB);CHKERRQ(ierr);
1250   if (graph->n_local_subs) {
1251     ierr = PetscFree(graph->local_subs);CHKERRQ(ierr);
1252   }
1253   graph->has_dirichlet       = PETSC_FALSE;
1254   graph->twodimset           = PETSC_FALSE;
1255   graph->twodim              = PETSC_FALSE;
1256   graph->nvtxs               = 0;
1257   graph->nvtxs_global        = 0;
1258   graph->n_subsets           = 0;
1259   graph->custom_minimal_size = 1;
1260   graph->n_local_subs        = 0;
1261   graph->maxcount            = PETSC_MAX_INT;
1262   graph->setupcalled         = PETSC_FALSE;
1263   PetscFunctionReturn(0);
1264 }
1265 
PCBDDCGraphInit(PCBDDCGraph graph,ISLocalToGlobalMapping l2gmap,PetscInt N,PetscInt maxcount)1266 PetscErrorCode PCBDDCGraphInit(PCBDDCGraph graph, ISLocalToGlobalMapping l2gmap, PetscInt N, PetscInt maxcount)
1267 {
1268   PetscInt       n;
1269   PetscErrorCode ierr;
1270 
1271   PetscFunctionBegin;
1272   PetscValidPointer(graph,1);
1273   PetscValidHeaderSpecific(l2gmap,IS_LTOGM_CLASSID,2);
1274   PetscValidLogicalCollectiveInt(l2gmap,N,3);
1275   PetscValidLogicalCollectiveInt(l2gmap,maxcount,4);
1276   /* raise an error if already allocated */
1277   if (graph->nvtxs_global) SETERRQ(PetscObjectComm((PetscObject)l2gmap),PETSC_ERR_PLIB,"BDDCGraph already initialized");
1278   /* set number of vertices */
1279   ierr = PetscObjectReference((PetscObject)l2gmap);CHKERRQ(ierr);
1280   graph->l2gmap = l2gmap;
1281   ierr = ISLocalToGlobalMappingGetSize(l2gmap,&n);CHKERRQ(ierr);
1282   graph->nvtxs = n;
1283   graph->nvtxs_global = N;
1284   /* allocate used space */
1285   ierr = PetscBTCreate(graph->nvtxs,&graph->touched);CHKERRQ(ierr);
1286   ierr = PetscMalloc5(graph->nvtxs,&graph->count,graph->nvtxs,&graph->neighbours_set,graph->nvtxs,&graph->subset,graph->nvtxs,&graph->which_dof,graph->nvtxs,&graph->special_dof);CHKERRQ(ierr);
1287   /* zeroes memory */
1288   ierr = PetscArrayzero(graph->count,graph->nvtxs);CHKERRQ(ierr);
1289   ierr = PetscArrayzero(graph->subset,graph->nvtxs);CHKERRQ(ierr);
1290   /* use -1 as a default value for which_dof array */
1291   for (n=0;n<graph->nvtxs;n++) graph->which_dof[n] = -1;
1292   ierr = PetscArrayzero(graph->special_dof,graph->nvtxs);CHKERRQ(ierr);
1293   /* zeroes first pointer to neighbour set */
1294   if (graph->nvtxs) {
1295     graph->neighbours_set[0] = NULL;
1296   }
1297   /* zeroes workspace for values of ncc */
1298   graph->subset_ncc = NULL;
1299   graph->subset_ref_node = NULL;
1300   /* maxcount for cc */
1301   graph->maxcount = maxcount;
1302   PetscFunctionReturn(0);
1303 }
1304 
PCBDDCGraphDestroy(PCBDDCGraph * graph)1305 PetscErrorCode PCBDDCGraphDestroy(PCBDDCGraph* graph)
1306 {
1307   PetscErrorCode ierr;
1308 
1309   PetscFunctionBegin;
1310   ierr = PCBDDCGraphResetCSR(*graph);CHKERRQ(ierr);
1311   ierr = PCBDDCGraphResetCoords(*graph);CHKERRQ(ierr);
1312   ierr = PCBDDCGraphReset(*graph);CHKERRQ(ierr);
1313   ierr = PetscFree(*graph);CHKERRQ(ierr);
1314   PetscFunctionReturn(0);
1315 }
1316 
PCBDDCGraphCreate(PCBDDCGraph * graph)1317 PetscErrorCode PCBDDCGraphCreate(PCBDDCGraph *graph)
1318 {
1319   PCBDDCGraph    new_graph;
1320   PetscErrorCode ierr;
1321 
1322   PetscFunctionBegin;
1323   ierr = PetscNew(&new_graph);CHKERRQ(ierr);
1324   new_graph->custom_minimal_size = 1;
1325   new_graph->commsizelimit = 1;
1326   *graph = new_graph;
1327   PetscFunctionReturn(0);
1328 }
1329