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