Actual source code: bddcgraph.c

  1: #include <petsc/private/petscimpl.h>
  2: #include <petsc/private/pcbddcprivateimpl.h>
  3: #include <petsc/private/pcbddcstructsimpl.h>

  5: PetscErrorCode PCBDDCDestroyGraphCandidatesIS(void *ctx)
  6: {
  7:   PCBDDCGraphCandidates cand = (PCBDDCGraphCandidates)ctx;

  9:   for (PetscInt i = 0; i < cand->nfc; i++) ISDestroy(&cand->Faces[i]);
 10:   for (PetscInt i = 0; i < cand->nec; i++) ISDestroy(&cand->Edges[i]);
 11:   PetscFree(cand->Faces);
 12:   PetscFree(cand->Edges);
 13:   ISDestroy(&cand->Vertices);
 14:   PetscFree(cand);
 15:   return 0;
 16: }

 18: PetscErrorCode PCBDDCGraphGetDirichletDofsB(PCBDDCGraph graph, IS *dirdofs)
 19: {
 20:   if (graph->dirdofsB) {
 21:     PetscObjectReference((PetscObject)graph->dirdofsB);
 22:   } else if (graph->has_dirichlet) {
 23:     PetscInt  i, size;
 24:     PetscInt *dirdofs_idxs;

 26:     size = 0;
 27:     for (i = 0; i < graph->nvtxs; i++) {
 28:       if (graph->count[i] && graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK) size++;
 29:     }

 31:     PetscMalloc1(size, &dirdofs_idxs);
 32:     size = 0;
 33:     for (i = 0; i < graph->nvtxs; i++) {
 34:       if (graph->count[i] && graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK) dirdofs_idxs[size++] = i;
 35:     }
 36:     ISCreateGeneral(PETSC_COMM_SELF, size, dirdofs_idxs, PETSC_OWN_POINTER, &graph->dirdofsB);
 37:     PetscObjectReference((PetscObject)graph->dirdofsB);
 38:   }
 39:   *dirdofs = graph->dirdofsB;
 40:   return 0;
 41: }

 43: PetscErrorCode PCBDDCGraphGetDirichletDofs(PCBDDCGraph graph, IS *dirdofs)
 44: {
 45:   if (graph->dirdofs) {
 46:     PetscObjectReference((PetscObject)graph->dirdofs);
 47:   } else if (graph->has_dirichlet) {
 48:     PetscInt  i, size;
 49:     PetscInt *dirdofs_idxs;

 51:     size = 0;
 52:     for (i = 0; i < graph->nvtxs; i++) {
 53:       if (graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK) size++;
 54:     }

 56:     PetscMalloc1(size, &dirdofs_idxs);
 57:     size = 0;
 58:     for (i = 0; i < graph->nvtxs; i++) {
 59:       if (graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK) dirdofs_idxs[size++] = i;
 60:     }
 61:     ISCreateGeneral(PetscObjectComm((PetscObject)graph->l2gmap), size, dirdofs_idxs, PETSC_OWN_POINTER, &graph->dirdofs);
 62:     PetscObjectReference((PetscObject)graph->dirdofs);
 63:   }
 64:   *dirdofs = graph->dirdofs;
 65:   return 0;
 66: }

 68: PetscErrorCode PCBDDCGraphASCIIView(PCBDDCGraph graph, PetscInt verbosity_level, PetscViewer viewer)
 69: {
 70:   PetscInt  i, j, tabs;
 71:   PetscInt *queue_in_global_numbering;

 73:   PetscViewerASCIIPushSynchronized(viewer);
 74:   PetscViewerASCIIGetTab(viewer, &tabs);
 75:   PetscViewerASCIIPrintf(viewer, "--------------------------------------------------\n");
 76:   PetscViewerFlush(viewer);
 77:   PetscViewerASCIISynchronizedPrintf(viewer, "Local BDDC graph for subdomain %04d\n", PetscGlobalRank);
 78:   PetscViewerASCIISynchronizedPrintf(viewer, "Number of vertices %" PetscInt_FMT "\n", graph->nvtxs);
 79:   PetscViewerASCIISynchronizedPrintf(viewer, "Number of local subdomains %" PetscInt_FMT "\n", graph->n_local_subs ? graph->n_local_subs : 1);
 80:   PetscViewerASCIISynchronizedPrintf(viewer, "Custom minimal size %" PetscInt_FMT "\n", graph->custom_minimal_size);
 81:   if (graph->maxcount != PETSC_MAX_INT) PetscViewerASCIISynchronizedPrintf(viewer, "Max count %" PetscInt_FMT "\n", graph->maxcount);
 82:   PetscViewerASCIISynchronizedPrintf(viewer, "Topological two dim? %s (set %s)\n", PetscBools[graph->twodim], PetscBools[graph->twodimset]);
 83:   if (verbosity_level > 2) {
 84:     for (i = 0; i < graph->nvtxs; i++) {
 85:       PetscViewerASCIISynchronizedPrintf(viewer, "%" PetscInt_FMT ":\n", i);
 86:       PetscViewerASCIISynchronizedPrintf(viewer, "   which_dof: %" PetscInt_FMT "\n", graph->which_dof[i]);
 87:       PetscViewerASCIISynchronizedPrintf(viewer, "   special_dof: %" PetscInt_FMT "\n", graph->special_dof[i]);
 88:       PetscViewerASCIISynchronizedPrintf(viewer, "   neighbours: %" PetscInt_FMT "\n", graph->count[i]);
 89:       PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
 90:       if (graph->count[i]) {
 91:         PetscViewerASCIISynchronizedPrintf(viewer, "     set of neighbours:");
 92:         for (j = 0; j < graph->count[i]; j++) PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, graph->neighbours_set[i][j]);
 93:         PetscViewerASCIISynchronizedPrintf(viewer, "\n");
 94:       }
 95:       PetscViewerASCIISetTab(viewer, tabs);
 96:       PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
 97:       if (graph->mirrors) {
 98:         PetscViewerASCIISynchronizedPrintf(viewer, "   mirrors: %" PetscInt_FMT "\n", graph->mirrors[i]);
 99:         if (graph->mirrors[i]) {
100:           PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
101:           PetscViewerASCIISynchronizedPrintf(viewer, "     set of mirrors:");
102:           for (j = 0; j < graph->mirrors[i]; j++) PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, graph->mirrors_set[i][j]);
103:           PetscViewerASCIISynchronizedPrintf(viewer, "\n");
104:           PetscViewerASCIISetTab(viewer, tabs);
105:           PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
106:         }
107:       }
108:       if (verbosity_level > 3) {
109:         if (graph->xadj) {
110:           PetscViewerASCIISynchronizedPrintf(viewer, "   local adj list:");
111:           PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
112:           for (j = graph->xadj[i]; j < graph->xadj[i + 1]; j++) PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, graph->adjncy[j]);
113:           PetscViewerASCIISynchronizedPrintf(viewer, "\n");
114:           PetscViewerASCIISetTab(viewer, tabs);
115:           PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
116:         } else {
117:           PetscViewerASCIISynchronizedPrintf(viewer, "   no adj info\n");
118:         }
119:       }
120:       if (graph->n_local_subs) PetscViewerASCIISynchronizedPrintf(viewer, "   local sub id: %" PetscInt_FMT "\n", graph->local_subs[i]);
121:       PetscViewerASCIISynchronizedPrintf(viewer, "   interface subset id: %" PetscInt_FMT "\n", graph->subset[i]);
122:       if (graph->subset[i] && graph->subset_ncc) PetscViewerASCIISynchronizedPrintf(viewer, "   ncc for subset: %" PetscInt_FMT "\n", graph->subset_ncc[graph->subset[i] - 1]);
123:     }
124:   }
125:   PetscViewerASCIISynchronizedPrintf(viewer, "Total number of connected components %" PetscInt_FMT "\n", graph->ncc);
126:   PetscMalloc1(graph->cptr[graph->ncc], &queue_in_global_numbering);
127:   ISLocalToGlobalMappingApply(graph->l2gmap, graph->cptr[graph->ncc], graph->queue, queue_in_global_numbering);
128:   for (i = 0; i < graph->ncc; i++) {
129:     PetscInt  node_num = graph->queue[graph->cptr[i]];
130:     PetscBool printcc  = PETSC_FALSE;
131:     PetscViewerASCIISynchronizedPrintf(viewer, "  cc %" PetscInt_FMT " (size %" PetscInt_FMT ", fid %" PetscInt_FMT ", neighs:", i, graph->cptr[i + 1] - graph->cptr[i], graph->which_dof[node_num]);
132:     PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
133:     for (j = 0; j < graph->count[node_num]; j++) PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, graph->neighbours_set[node_num][j]);
134:     if (verbosity_level > 1) {
135:       PetscViewerASCIISynchronizedPrintf(viewer, "):");
136:       if (verbosity_level > 2 || graph->twodim || graph->count[node_num] > 1 || (graph->count[node_num] == 1 && graph->special_dof[node_num] == PCBDDCGRAPH_NEUMANN_MARK)) printcc = PETSC_TRUE;
137:       if (printcc) {
138:         for (j = graph->cptr[i]; j < graph->cptr[i + 1]; j++) PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT " (%" PetscInt_FMT ")", graph->queue[j], queue_in_global_numbering[j]);
139:       }
140:     } else {
141:       PetscViewerASCIISynchronizedPrintf(viewer, ")");
142:     }
143:     PetscViewerASCIISynchronizedPrintf(viewer, "\n");
144:     PetscViewerASCIISetTab(viewer, tabs);
145:     PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
146:   }
147:   PetscFree(queue_in_global_numbering);
148:   PetscViewerFlush(viewer);
149:   return 0;
150: }

152: PetscErrorCode PCBDDCGraphRestoreCandidatesIS(PCBDDCGraph graph, PetscInt *n_faces, IS *FacesIS[], PetscInt *n_edges, IS *EdgesIS[], IS *VerticesIS)
153: {
154:   PetscInt       i;
155:   PetscContainer gcand;

157:   PetscObjectQuery((PetscObject)graph->l2gmap, "_PCBDDCGraphCandidatesIS", (PetscObject *)&gcand);
158:   if (gcand) {
159:     if (n_faces) *n_faces = 0;
160:     if (n_edges) *n_edges = 0;
161:     if (FacesIS) *FacesIS = NULL;
162:     if (EdgesIS) *EdgesIS = NULL;
163:     if (VerticesIS) *VerticesIS = NULL;
164:   }
165:   if (n_faces) {
166:     if (FacesIS) {
167:       for (i = 0; i < *n_faces; i++) ISDestroy(&((*FacesIS)[i]));
168:       PetscFree(*FacesIS);
169:     }
170:     *n_faces = 0;
171:   }
172:   if (n_edges) {
173:     if (EdgesIS) {
174:       for (i = 0; i < *n_edges; i++) ISDestroy(&((*EdgesIS)[i]));
175:       PetscFree(*EdgesIS);
176:     }
177:     *n_edges = 0;
178:   }
179:   if (VerticesIS) ISDestroy(VerticesIS);
180:   return 0;
181: }

183: PetscErrorCode PCBDDCGraphGetCandidatesIS(PCBDDCGraph graph, PetscInt *n_faces, IS *FacesIS[], PetscInt *n_edges, IS *EdgesIS[], IS *VerticesIS)
184: {
185:   IS            *ISForFaces, *ISForEdges, ISForVertices;
186:   PetscInt       i, nfc, nec, nvc, *idx, *mark;
187:   PetscContainer gcand;

189:   PetscObjectQuery((PetscObject)graph->l2gmap, "_PCBDDCGraphCandidatesIS", (PetscObject *)&gcand);
190:   if (gcand) {
191:     PCBDDCGraphCandidates cand;

193:     PetscContainerGetPointer(gcand, (void **)&cand);
194:     if (n_faces) *n_faces = cand->nfc;
195:     if (FacesIS) *FacesIS = cand->Faces;
196:     if (n_edges) *n_edges = cand->nec;
197:     if (EdgesIS) *EdgesIS = cand->Edges;
198:     if (VerticesIS) *VerticesIS = cand->Vertices;
199:     return 0;
200:   }
201:   PetscCalloc1(graph->ncc, &mark);
202:   /* loop on ccs to evaluate 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:   }

221:   /* allocate IS arrays for faces, edges. Vertices need a single index set. */
222:   if (FacesIS) PetscMalloc1(nfc, &ISForFaces);
223:   if (EdgesIS) PetscMalloc1(nec, &ISForEdges);
224:   if (VerticesIS) PetscMalloc1(nvc, &idx);

226:   /* loop on ccs to compute index sets for faces and edges */
227:   if (!graph->queue_sorted) {
228:     PetscInt *queue_global;

230:     PetscMalloc1(graph->cptr[graph->ncc], &queue_global);
231:     ISLocalToGlobalMappingApply(graph->l2gmap, graph->cptr[graph->ncc], graph->queue, queue_global);
232:     for (i = 0; i < graph->ncc; i++) PetscSortIntWithArray(graph->cptr[i + 1] - graph->cptr[i], &queue_global[graph->cptr[i]], &graph->queue[graph->cptr[i]]);
233:     PetscFree(queue_global);
234:     graph->queue_sorted = PETSC_TRUE;
235:   }
236:   nfc = 0;
237:   nec = 0;
238:   for (i = 0; i < graph->ncc; i++) {
239:     if (mark[i] == 2) {
240:       if (FacesIS) ISCreateGeneral(PETSC_COMM_SELF, graph->cptr[i + 1] - graph->cptr[i], &graph->queue[graph->cptr[i]], PETSC_USE_POINTER, &ISForFaces[nfc]);
241:       nfc++;
242:     } else if (mark[i] == 1) {
243:       if (EdgesIS) ISCreateGeneral(PETSC_COMM_SELF, graph->cptr[i + 1] - graph->cptr[i], &graph->queue[graph->cptr[i]], PETSC_USE_POINTER, &ISForEdges[nec]);
244:       nec++;
245:     }
246:   }

248:   /* index set for vertices */
249:   if (VerticesIS) {
250:     nvc = 0;
251:     for (i = 0; i < graph->ncc; i++) {
252:       if (!mark[i]) {
253:         PetscInt j;

255:         for (j = graph->cptr[i]; j < graph->cptr[i + 1]; j++) {
256:           idx[nvc] = graph->queue[j];
257:           nvc++;
258:         }
259:       }
260:     }
261:     /* sort vertex set (by local ordering) */
262:     PetscSortInt(nvc, idx);
263:     ISCreateGeneral(PETSC_COMM_SELF, nvc, idx, PETSC_OWN_POINTER, &ISForVertices);
264:   }
265:   PetscFree(mark);

267:   /* get back info */
268:   if (n_faces) *n_faces = nfc;
269:   if (FacesIS) *FacesIS = ISForFaces;
270:   if (n_edges) *n_edges = nec;
271:   if (EdgesIS) *EdgesIS = ISForEdges;
272:   if (VerticesIS) *VerticesIS = ISForVertices;
273:   return 0;
274: }

276: PetscErrorCode PCBDDCGraphComputeConnectedComponents(PCBDDCGraph graph)
277: {
278:   PetscBool   adapt_interface_reduced;
279:   MPI_Comm    interface_comm;
280:   PetscMPIInt size;
281:   PetscInt    i;
282:   PetscBT     cornerp;

284:   /* compute connected components locally */
285:   PetscObjectGetComm((PetscObject)(graph->l2gmap), &interface_comm);
286:   PCBDDCGraphComputeConnectedComponentsLocal(graph);

288:   cornerp = NULL;
289:   if (graph->active_coords) { /* face based corner selection */
290:     PetscBT    excluded;
291:     PetscReal *wdist;
292:     PetscInt   n_neigh, *neigh, *n_shared, **shared;
293:     PetscInt   maxc, ns;

295:     PetscBTCreate(graph->nvtxs, &cornerp);
296:     ISLocalToGlobalMappingGetInfo(graph->l2gmap, &n_neigh, &neigh, &n_shared, &shared);
297:     for (ns = 1, maxc = 0; ns < n_neigh; ns++) maxc = PetscMax(maxc, n_shared[ns]);
298:     PetscMalloc1(maxc * graph->cdim, &wdist);
299:     PetscBTCreate(maxc, &excluded);

301:     for (ns = 1; ns < n_neigh; ns++) { /* first proc is self */
302:       PetscReal *anchor, mdist;
303:       PetscInt   fst, j, k, d, cdim = graph->cdim, n = n_shared[ns];
304:       PetscInt   point1, point2, point3, point4;

306:       /* import coordinates on shared interface */
307:       PetscBTMemzero(n, excluded);
308:       for (j = 0, fst = -1, k = 0; j < n; j++) {
309:         PetscBool skip = PETSC_FALSE;
310:         for (d = 0; d < cdim; d++) {
311:           PetscReal c = graph->coords[shared[ns][j] * cdim + d];
312:           skip        = (PetscBool)(skip || c == PETSC_MAX_REAL);
313:           wdist[k++]  = c;
314:         }
315:         if (skip) {
316:           PetscBTSet(excluded, j);
317:         } else if (fst == -1) fst = j;
318:       }
319:       if (fst == -1) continue;

321:       /* the dofs are sorted by global numbering, so each rank starts from the same id
322:          and it will detect the same corners from the given set */

324:       /* find the farthest point from the starting one */
325:       anchor = wdist + fst * cdim;
326:       mdist  = -1.0;
327:       point1 = fst;
328:       for (j = fst; j < n; j++) {
329:         PetscReal dist = 0.0;

331:         if (PetscUnlikely(PetscBTLookup(excluded, j))) continue;
332:         for (d = 0; d < cdim; d++) dist += (wdist[j * cdim + d] - anchor[d]) * (wdist[j * cdim + d] - anchor[d]);
333:         if (dist > mdist) {
334:           mdist  = dist;
335:           point1 = j;
336:         }
337:       }

339:       /* find the farthest point from point1 */
340:       anchor = wdist + point1 * cdim;
341:       mdist  = -1.0;
342:       point2 = point1;
343:       for (j = fst; j < n; j++) {
344:         PetscReal dist = 0.0;

346:         if (PetscUnlikely(PetscBTLookup(excluded, j))) continue;
347:         for (d = 0; d < cdim; d++) dist += (wdist[j * cdim + d] - anchor[d]) * (wdist[j * cdim + d] - anchor[d]);
348:         if (dist > mdist) {
349:           mdist  = dist;
350:           point2 = j;
351:         }
352:       }

354:       /* find the third point maximizing the triangle area */
355:       point3 = point2;
356:       if (cdim > 2) {
357:         PetscReal a = 0.0;

359:         for (d = 0; d < cdim; d++) a += (wdist[point1 * cdim + d] - wdist[point2 * cdim + d]) * (wdist[point1 * cdim + d] - wdist[point2 * cdim + d]);
360:         a     = PetscSqrtReal(a);
361:         mdist = -1.0;
362:         for (j = fst; j < n; j++) {
363:           PetscReal area, b = 0.0, c = 0.0, s;

365:           if (PetscUnlikely(PetscBTLookup(excluded, j))) continue;
366:           for (d = 0; d < cdim; d++) {
367:             b += (wdist[point1 * cdim + d] - wdist[j * cdim + d]) * (wdist[point1 * cdim + d] - wdist[j * cdim + d]);
368:             c += (wdist[point2 * cdim + d] - wdist[j * cdim + d]) * (wdist[point2 * cdim + d] - wdist[j * cdim + d]);
369:           }
370:           b = PetscSqrtReal(b);
371:           c = PetscSqrtReal(c);
372:           s = 0.5 * (a + b + c);

374:           /* Heron's formula, area squared */
375:           area = s * (s - a) * (s - b) * (s - c);
376:           if (area > mdist) {
377:             mdist  = area;
378:             point3 = j;
379:           }
380:         }
381:       }

383:       /* find the farthest point from point3 different from point1 and point2 */
384:       anchor = wdist + point3 * cdim;
385:       mdist  = -1.0;
386:       point4 = point3;
387:       for (j = fst; j < n; j++) {
388:         PetscReal dist = 0.0;

390:         if (PetscUnlikely(PetscBTLookup(excluded, j)) || j == point1 || j == point2 || j == point3) continue;
391:         for (d = 0; d < cdim; d++) dist += (wdist[j * cdim + d] - anchor[d]) * (wdist[j * cdim + d] - anchor[d]);
392:         if (dist > mdist) {
393:           mdist  = dist;
394:           point4 = j;
395:         }
396:       }

398:       PetscBTSet(cornerp, shared[ns][point1]);
399:       PetscBTSet(cornerp, shared[ns][point2]);
400:       PetscBTSet(cornerp, shared[ns][point3]);
401:       PetscBTSet(cornerp, shared[ns][point4]);

403:       /* all dofs having the same coordinates will be primal */
404:       for (j = fst; j < n; j++) {
405:         PetscBool same[] = {PETSC_TRUE, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE};

407:         if (PetscUnlikely(PetscBTLookup(excluded, j))) continue;
408:         for (d = 0; d < cdim; d++) {
409:           same[0] = (PetscBool)(same[0] && (PetscAbsReal(wdist[j * cdim + d] - wdist[point1 * cdim + d]) < PETSC_SMALL));
410:           same[1] = (PetscBool)(same[1] && (PetscAbsReal(wdist[j * cdim + d] - wdist[point2 * cdim + d]) < PETSC_SMALL));
411:           same[2] = (PetscBool)(same[2] && (PetscAbsReal(wdist[j * cdim + d] - wdist[point3 * cdim + d]) < PETSC_SMALL));
412:           same[3] = (PetscBool)(same[3] && (PetscAbsReal(wdist[j * cdim + d] - wdist[point4 * cdim + d]) < PETSC_SMALL));
413:         }
414:         if (same[0] || same[1] || same[2] || same[3]) PetscBTSet(cornerp, shared[ns][j]);
415:       }
416:     }
417:     PetscBTDestroy(&excluded);
418:     PetscFree(wdist);
419:     ISLocalToGlobalMappingRestoreInfo(graph->l2gmap, &n_neigh, &neigh, &n_shared, &shared);
420:   }

422:   /* check consistency of connected components among neighbouring subdomains -> it adapt them in case it is needed */
423:   MPI_Comm_size(interface_comm, &size);
424:   adapt_interface_reduced = PETSC_FALSE;
425:   if (size > 1) {
426:     PetscInt  i;
427:     PetscBool adapt_interface = cornerp ? PETSC_TRUE : PETSC_FALSE;
428:     for (i = 0; i < graph->n_subsets && !adapt_interface; i++) {
429:       /* We are not sure that on a given subset of the local interface,
430:          with two connected components, the latters be the same among sharing subdomains */
431:       if (graph->subset_ncc[i] > 1) adapt_interface = PETSC_TRUE;
432:     }
433:     MPIU_Allreduce(&adapt_interface, &adapt_interface_reduced, 1, MPIU_BOOL, MPI_LOR, interface_comm);
434:   }

436:   if (graph->n_subsets && adapt_interface_reduced) {
437:     PetscBT      subset_cc_adapt;
438:     MPI_Request *send_requests, *recv_requests;
439:     PetscInt    *send_buffer, *recv_buffer;
440:     PetscInt     sum_requests, start_of_recv, start_of_send;
441:     PetscInt    *cum_recv_counts;
442:     PetscInt    *labels;
443:     PetscInt     ncc, cum_queue, mss, mns, j, k, s;
444:     PetscInt   **refine_buffer = NULL, *private_labels = NULL;
445:     PetscBool   *subset_has_corn, *recv_buffer_bool, *send_buffer_bool;

447:     PetscCalloc1(graph->n_subsets, &subset_has_corn);
448:     if (cornerp) {
449:       for (i = 0; i < graph->n_subsets; i++) {
450:         for (j = 0; j < graph->subset_size[i]; j++) {
451:           if (PetscBTLookup(cornerp, graph->subset_idxs[i][j])) {
452:             subset_has_corn[i] = PETSC_TRUE;
453:             break;
454:           }
455:         }
456:       }
457:     }
458:     PetscMalloc1(graph->nvtxs, &labels);
459:     PetscArrayzero(labels, graph->nvtxs);
460:     for (i = 0, k = 0; i < graph->ncc; i++) {
461:       PetscInt s = 1;
462:       for (j = graph->cptr[i]; j < graph->cptr[i + 1]; j++) {
463:         if (cornerp && PetscBTLookup(cornerp, graph->queue[j])) {
464:           labels[graph->queue[j]] = k + s;
465:           s += 1;
466:         } else {
467:           labels[graph->queue[j]] = k;
468:         }
469:       }
470:       k += s;
471:     }

473:     /* allocate some space */
474:     PetscMalloc1(graph->n_subsets + 1, &cum_recv_counts);
475:     PetscArrayzero(cum_recv_counts, graph->n_subsets + 1);

477:     /* first count how many neighbours per connected component I will receive from */
478:     cum_recv_counts[0] = 0;
479:     for (i = 0; i < graph->n_subsets; i++) cum_recv_counts[i + 1] = cum_recv_counts[i] + graph->count[graph->subset_idxs[i][0]];
480:     PetscMalloc1(graph->n_subsets, &send_buffer_bool);
481:     PetscMalloc1(cum_recv_counts[graph->n_subsets], &recv_buffer_bool);
482:     PetscMalloc2(cum_recv_counts[graph->n_subsets], &send_requests, cum_recv_counts[graph->n_subsets], &recv_requests);
483:     for (i = 0; i < cum_recv_counts[graph->n_subsets]; i++) {
484:       send_requests[i] = MPI_REQUEST_NULL;
485:       recv_requests[i] = MPI_REQUEST_NULL;
486:     }

488:     /* exchange with my neighbours the number of my connected components on the subset of interface */
489:     sum_requests = 0;
490:     for (i = 0; i < graph->n_subsets; i++) send_buffer_bool[i] = (PetscBool)(graph->subset_ncc[i] > 1 || subset_has_corn[i]);
491:     for (i = 0; i < graph->n_subsets; i++) {
492:       PetscMPIInt neigh, tag;
493:       PetscInt    count, *neighs;

495:       count  = graph->count[graph->subset_idxs[i][0]];
496:       neighs = graph->neighbours_set[graph->subset_idxs[i][0]];
497:       PetscMPIIntCast(2 * graph->subset_ref_node[i], &tag);
498:       for (k = 0; k < count; k++) {
499:         PetscMPIIntCast(neighs[k], &neigh);
500:         MPI_Isend(send_buffer_bool + i, 1, MPIU_BOOL, neigh, tag, interface_comm, &send_requests[sum_requests]);
501:         MPI_Irecv(recv_buffer_bool + sum_requests, 1, MPIU_BOOL, neigh, tag, interface_comm, &recv_requests[sum_requests]);
502:         sum_requests++;
503:       }
504:     }
505:     MPI_Waitall(sum_requests, recv_requests, MPI_STATUSES_IGNORE);
506:     MPI_Waitall(sum_requests, send_requests, MPI_STATUSES_IGNORE);

508:     /* determine the subsets I have to adapt (those having more than 1 cc) */
509:     PetscBTCreate(graph->n_subsets, &subset_cc_adapt);
510:     PetscBTMemzero(graph->n_subsets, subset_cc_adapt);
511:     for (i = 0; i < graph->n_subsets; i++) {
512:       if (graph->subset_ncc[i] > 1 || subset_has_corn[i]) {
513:         PetscBTSet(subset_cc_adapt, i);
514:         continue;
515:       }
516:       for (j = cum_recv_counts[i]; j < cum_recv_counts[i + 1]; j++) {
517:         if (recv_buffer_bool[j]) {
518:           PetscBTSet(subset_cc_adapt, i);
519:           break;
520:         }
521:       }
522:     }
523:     PetscFree(send_buffer_bool);
524:     PetscFree(recv_buffer_bool);
525:     PetscFree(subset_has_corn);

527:     /* determine send/recv buffers sizes */
528:     j   = 0;
529:     mss = 0;
530:     for (i = 0; i < graph->n_subsets; i++) {
531:       if (PetscBTLookup(subset_cc_adapt, i)) {
532:         j += graph->subset_size[i];
533:         mss = PetscMax(graph->subset_size[i], mss);
534:       }
535:     }
536:     k   = 0;
537:     mns = 0;
538:     for (i = 0; i < graph->n_subsets; i++) {
539:       if (PetscBTLookup(subset_cc_adapt, i)) {
540:         k += (cum_recv_counts[i + 1] - cum_recv_counts[i]) * graph->subset_size[i];
541:         mns = PetscMax(cum_recv_counts[i + 1] - cum_recv_counts[i], mns);
542:       }
543:     }
544:     PetscMalloc2(j, &send_buffer, k, &recv_buffer);

546:     /* fill send buffer (order matters: subset_idxs ordered by global ordering) */
547:     j = 0;
548:     for (i = 0; i < graph->n_subsets; i++)
549:       if (PetscBTLookup(subset_cc_adapt, i))
550:         for (k = 0; k < graph->subset_size[i]; k++) send_buffer[j++] = labels[graph->subset_idxs[i][k]];

552:     /* now exchange the data */
553:     start_of_recv = 0;
554:     start_of_send = 0;
555:     sum_requests  = 0;
556:     for (i = 0; i < graph->n_subsets; i++) {
557:       if (PetscBTLookup(subset_cc_adapt, i)) {
558:         PetscMPIInt neigh, tag;
559:         PetscInt    size_of_send = graph->subset_size[i];

561:         j = graph->subset_idxs[i][0];
562:         PetscMPIIntCast(2 * graph->subset_ref_node[i] + 1, &tag);
563:         for (k = 0; k < graph->count[j]; k++) {
564:           PetscMPIIntCast(graph->neighbours_set[j][k], &neigh);
565:           MPI_Isend(&send_buffer[start_of_send], size_of_send, MPIU_INT, neigh, tag, interface_comm, &send_requests[sum_requests]);
566:           MPI_Irecv(&recv_buffer[start_of_recv], size_of_send, MPIU_INT, neigh, tag, interface_comm, &recv_requests[sum_requests]);
567:           start_of_recv += size_of_send;
568:           sum_requests++;
569:         }
570:         start_of_send += size_of_send;
571:       }
572:     }
573:     MPI_Waitall(sum_requests, recv_requests, MPI_STATUSES_IGNORE);

575:     /* refine connected components */
576:     start_of_recv = 0;
577:     /* allocate some temporary space */
578:     if (mss) {
579:       PetscMalloc1(mss, &refine_buffer);
580:       PetscMalloc2(mss * (mns + 1), &refine_buffer[0], mss, &private_labels);
581:     }
582:     ncc            = 0;
583:     cum_queue      = 0;
584:     graph->cptr[0] = 0;
585:     for (i = 0; i < graph->n_subsets; i++) {
586:       if (PetscBTLookup(subset_cc_adapt, i)) {
587:         PetscInt subset_counter = 0;
588:         PetscInt sharingprocs   = cum_recv_counts[i + 1] - cum_recv_counts[i] + 1; /* count myself */
589:         PetscInt buffer_size    = graph->subset_size[i];

591:         /* compute pointers */
592:         for (j = 1; j < buffer_size; j++) refine_buffer[j] = refine_buffer[j - 1] + sharingprocs;
593:         /* analyze contributions from subdomains that share the i-th subset
594:            The structure of refine_buffer is suitable to find intersections of ccs among sharingprocs.
595:            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)
596:            sharing procs connected components:
597:              neigh 0: [0 1 4], [2 3], labels [4,7]  (2 connected components)
598:              neigh 1: [0 1], [2 3 4], labels [3 2]  (2 connected components)
599:              neigh 2: [0 4], [1], [2 3], labels [1 5 6] (3 connected components)
600:            refine_buffer will be filled as:
601:              [ 4, 3, 1;
602:                4, 2, 1;
603:                7, 2, 6;
604:                4, 3, 5;
605:                7, 2, 6; ];
606:            The connected components in local ordering are [0], [1], [2 3], [4] */
607:         /* fill temp_buffer */
608:         for (k = 0; k < buffer_size; k++) refine_buffer[k][0] = labels[graph->subset_idxs[i][k]];
609:         for (j = 0; j < sharingprocs - 1; j++) {
610:           for (k = 0; k < buffer_size; k++) refine_buffer[k][j + 1] = recv_buffer[start_of_recv + k];
611:           start_of_recv += buffer_size;
612:         }
613:         PetscArrayzero(private_labels, buffer_size);
614:         for (j = 0; j < buffer_size; j++) {
615:           if (!private_labels[j]) { /* found a new cc  */
616:             PetscBool same_set;

618:             graph->cptr[ncc] = cum_queue;
619:             ncc++;
620:             subset_counter++;
621:             private_labels[j]         = subset_counter;
622:             graph->queue[cum_queue++] = graph->subset_idxs[i][j];
623:             for (k = j + 1; k < buffer_size; k++) { /* check for other nodes in new cc */
624:               same_set = PETSC_TRUE;
625:               for (s = 0; s < sharingprocs; s++) {
626:                 if (refine_buffer[j][s] != refine_buffer[k][s]) {
627:                   same_set = PETSC_FALSE;
628:                   break;
629:                 }
630:               }
631:               if (same_set) {
632:                 private_labels[k]         = subset_counter;
633:                 graph->queue[cum_queue++] = graph->subset_idxs[i][k];
634:               }
635:             }
636:           }
637:         }
638:         graph->cptr[ncc]     = cum_queue;
639:         graph->subset_ncc[i] = subset_counter;
640:         graph->queue_sorted  = PETSC_FALSE;
641:       } else { /* this subset does not need to be adapted */
642:         PetscArraycpy(graph->queue + cum_queue, graph->subset_idxs[i], graph->subset_size[i]);
643:         ncc++;
644:         cum_queue += graph->subset_size[i];
645:         graph->cptr[ncc] = cum_queue;
646:       }
647:     }
648:     graph->cptr[ncc] = cum_queue;
649:     graph->ncc       = ncc;
650:     if (mss) {
651:       PetscFree2(refine_buffer[0], private_labels);
652:       PetscFree(refine_buffer);
653:     }
654:     PetscFree(labels);
655:     MPI_Waitall(sum_requests, send_requests, MPI_STATUSES_IGNORE);
656:     PetscFree2(send_requests, recv_requests);
657:     PetscFree2(send_buffer, recv_buffer);
658:     PetscFree(cum_recv_counts);
659:     PetscBTDestroy(&subset_cc_adapt);
660:   }
661:   PetscBTDestroy(&cornerp);

663:   /* Determine if we are in 2D or 3D */
664:   if (!graph->twodimset) {
665:     PetscBool twodim = PETSC_TRUE;
666:     for (i = 0; i < graph->ncc; i++) {
667:       PetscInt repdof = graph->queue[graph->cptr[i]];
668:       PetscInt ccsize = graph->cptr[i + 1] - graph->cptr[i];
669:       if (graph->count[repdof] > 1 && ccsize > graph->custom_minimal_size) {
670:         twodim = PETSC_FALSE;
671:         break;
672:       }
673:     }
674:     MPIU_Allreduce(&twodim, &graph->twodim, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)graph->l2gmap));
675:     graph->twodimset = PETSC_TRUE;
676:   }
677:   return 0;
678: }

680: static inline PetscErrorCode PCBDDCGraphComputeCC_Private(PCBDDCGraph graph, PetscInt pid, PetscInt *queue_tip, PetscInt n_prev, PetscInt *n_added)
681: {
682:   PetscInt  i, j, n;
683:   PetscInt *xadj = graph->xadj, *adjncy = graph->adjncy;
684:   PetscBT   touched  = graph->touched;
685:   PetscBool havecsr  = (PetscBool)(!!xadj);
686:   PetscBool havesubs = (PetscBool)(!!graph->n_local_subs);

688:   n = 0;
689:   if (havecsr && !havesubs) {
690:     for (i = -n_prev; i < 0; i++) {
691:       PetscInt start_dof = queue_tip[i];
692:       /* we assume that if a dof has a size 1 adjacency list and the corresponding entry is negative, it is connected to all dofs */
693:       if (xadj[start_dof + 1] - xadj[start_dof] == 1 && adjncy[xadj[start_dof]] < 0) {
694:         for (j = 0; j < graph->subset_size[pid - 1]; j++) { /* pid \in [1,graph->n_subsets] */
695:           PetscInt dof = graph->subset_idxs[pid - 1][j];
696:           if (!PetscBTLookup(touched, dof) && graph->subset[dof] == pid) {
697:             PetscBTSet(touched, dof);
698:             queue_tip[n] = dof;
699:             n++;
700:           }
701:         }
702:       } else {
703:         for (j = xadj[start_dof]; j < xadj[start_dof + 1]; j++) {
704:           PetscInt dof = adjncy[j];
705:           if (!PetscBTLookup(touched, dof) && graph->subset[dof] == pid) {
706:             PetscBTSet(touched, dof);
707:             queue_tip[n] = dof;
708:             n++;
709:           }
710:         }
711:       }
712:     }
713:   } else if (havecsr && havesubs) {
714:     PetscInt sid = graph->local_subs[queue_tip[-n_prev]];
715:     for (i = -n_prev; i < 0; i++) {
716:       PetscInt start_dof = queue_tip[i];
717:       /* 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 */
718:       if (xadj[start_dof + 1] - xadj[start_dof] == 1 && adjncy[xadj[start_dof]] < 0) {
719:         for (j = 0; j < graph->subset_size[pid - 1]; j++) { /* pid \in [1,graph->n_subsets] */
720:           PetscInt dof = graph->subset_idxs[pid - 1][j];
721:           if (!PetscBTLookup(touched, dof) && graph->subset[dof] == pid && graph->local_subs[dof] == sid) {
722:             PetscBTSet(touched, dof);
723:             queue_tip[n] = dof;
724:             n++;
725:           }
726:         }
727:       } else {
728:         for (j = xadj[start_dof]; j < xadj[start_dof + 1]; j++) {
729:           PetscInt dof = adjncy[j];
730:           if (!PetscBTLookup(touched, dof) && graph->subset[dof] == pid && graph->local_subs[dof] == sid) {
731:             PetscBTSet(touched, dof);
732:             queue_tip[n] = dof;
733:             n++;
734:           }
735:         }
736:       }
737:     }
738:   } else if (havesubs) { /* sub info only */
739:     PetscInt sid = graph->local_subs[queue_tip[-n_prev]];
740:     for (j = 0; j < graph->subset_size[pid - 1]; j++) { /* pid \in [1,graph->n_subsets] */
741:       PetscInt dof = graph->subset_idxs[pid - 1][j];
742:       if (!PetscBTLookup(touched, dof) && graph->subset[dof] == pid && graph->local_subs[dof] == sid) {
743:         PetscBTSet(touched, dof);
744:         queue_tip[n] = dof;
745:         n++;
746:       }
747:     }
748:   } else {
749:     for (j = 0; j < graph->subset_size[pid - 1]; j++) { /* pid \in [1,graph->n_subsets] */
750:       PetscInt dof = graph->subset_idxs[pid - 1][j];
751:       if (!PetscBTLookup(touched, dof) && graph->subset[dof] == pid) {
752:         PetscBTSet(touched, dof);
753:         queue_tip[n] = dof;
754:         n++;
755:       }
756:     }
757:   }
758:   *n_added = n;
759:   return 0;
760: }

762: PetscErrorCode PCBDDCGraphComputeConnectedComponentsLocal(PCBDDCGraph graph)
763: {
764:   PetscInt    ncc, cum_queue, n;
765:   PetscMPIInt commsize;

768:   /* quiet return if there isn't any local info */
769:   if (!graph->xadj && !graph->n_local_subs) return 0;

771:   /* reset any previous search of connected components */
772:   PetscBTMemzero(graph->nvtxs, graph->touched);
773:   MPI_Comm_size(PetscObjectComm((PetscObject)graph->l2gmap), &commsize);
774:   if (commsize > graph->commsizelimit) {
775:     PetscInt i;
776:     for (i = 0; i < graph->nvtxs; i++) {
777:       if (graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK || !graph->count[i]) PetscBTSet(graph->touched, i);
778:     }
779:   }

781:   /* begin search for connected components */
782:   cum_queue = 0;
783:   ncc       = 0;
784:   for (n = 0; n < graph->n_subsets; n++) {
785:     PetscInt pid   = n + 1; /* partition labeled by 0 is discarded */
786:     PetscInt found = 0, prev = 0, first = 0, ncc_pid = 0;
787:     while (found != graph->subset_size[n]) {
788:       PetscInt added = 0;
789:       if (!prev) { /* search for new starting dof */
790:         while (PetscBTLookup(graph->touched, graph->subset_idxs[n][first])) first++;
791:         PetscBTSet(graph->touched, graph->subset_idxs[n][first]);
792:         graph->queue[cum_queue] = graph->subset_idxs[n][first];
793:         graph->cptr[ncc]        = cum_queue;
794:         prev                    = 1;
795:         cum_queue++;
796:         found++;
797:         ncc_pid++;
798:         ncc++;
799:       }
800:       PCBDDCGraphComputeCC_Private(graph, pid, graph->queue + cum_queue, prev, &added);
801:       if (!added) {
802:         graph->subset_ncc[n] = ncc_pid;
803:         graph->cptr[ncc]     = cum_queue;
804:       }
805:       prev = added;
806:       found += added;
807:       cum_queue += added;
808:       if (added && found == graph->subset_size[n]) {
809:         graph->subset_ncc[n] = ncc_pid;
810:         graph->cptr[ncc]     = cum_queue;
811:       }
812:     }
813:   }
814:   graph->ncc          = ncc;
815:   graph->queue_sorted = PETSC_FALSE;
816:   return 0;
817: }

819: PetscErrorCode PCBDDCGraphSetUp(PCBDDCGraph graph, PetscInt custom_minimal_size, IS neumann_is, IS dirichlet_is, PetscInt n_ISForDofs, IS ISForDofs[], IS custom_primal_vertices)
820: {
821:   IS              subset, subset_n;
822:   MPI_Comm        comm;
823:   const PetscInt *is_indices;
824:   PetscInt        n_neigh, *neigh, *n_shared, **shared, *queue_global;
825:   PetscInt        i, j, k, s, total_counts, nodes_touched, is_size;
826:   PetscMPIInt     commsize;
827:   PetscBool       same_set, mirrors_found;

830:   if (neumann_is) {
833:   }
834:   graph->has_dirichlet = PETSC_FALSE;
835:   if (dirichlet_is) {
838:     graph->has_dirichlet = PETSC_TRUE;
839:   }
841:   for (i = 0; i < n_ISForDofs; i++) {
844:   }
845:   if (custom_primal_vertices) {
848:   }
849:   PetscObjectGetComm((PetscObject)(graph->l2gmap), &comm);
850:   MPI_Comm_size(comm, &commsize);

852:   /* custom_minimal_size */
853:   graph->custom_minimal_size = custom_minimal_size;
854:   /* get info l2gmap and allocate work vectors  */
855:   ISLocalToGlobalMappingGetInfo(graph->l2gmap, &n_neigh, &neigh, &n_shared, &shared);
856:   /* check if we have any local periodic nodes (periodic BCs) */
857:   mirrors_found = PETSC_FALSE;
858:   if (graph->nvtxs && n_neigh) {
859:     for (i = 0; i < n_shared[0]; i++) graph->count[shared[0][i]] += 1;
860:     for (i = 0; i < n_shared[0]; i++) {
861:       if (graph->count[shared[0][i]] > 1) {
862:         mirrors_found = PETSC_TRUE;
863:         break;
864:       }
865:     }
866:   }
867:   /* compute local mirrors (if any) */
868:   if (mirrors_found) {
869:     IS        to, from;
870:     PetscInt *local_indices, *global_indices;

872:     ISCreateStride(PETSC_COMM_SELF, graph->nvtxs, 0, 1, &to);
873:     ISLocalToGlobalMappingApplyIS(graph->l2gmap, to, &from);
874:     /* get arrays of local and global indices */
875:     PetscMalloc1(graph->nvtxs, &local_indices);
876:     ISGetIndices(to, (const PetscInt **)&is_indices);
877:     PetscArraycpy(local_indices, is_indices, graph->nvtxs);
878:     ISRestoreIndices(to, (const PetscInt **)&is_indices);
879:     PetscMalloc1(graph->nvtxs, &global_indices);
880:     ISGetIndices(from, (const PetscInt **)&is_indices);
881:     PetscArraycpy(global_indices, is_indices, graph->nvtxs);
882:     ISRestoreIndices(from, (const PetscInt **)&is_indices);
883:     /* allocate space for mirrors */
884:     PetscMalloc2(graph->nvtxs, &graph->mirrors, graph->nvtxs, &graph->mirrors_set);
885:     PetscArrayzero(graph->mirrors, graph->nvtxs);
886:     graph->mirrors_set[0] = NULL;

888:     k = 0;
889:     for (i = 0; i < n_shared[0]; i++) {
890:       j = shared[0][i];
891:       if (graph->count[j] > 1) {
892:         graph->mirrors[j]++;
893:         k++;
894:       }
895:     }
896:     /* allocate space for set of mirrors */
897:     PetscMalloc1(k, &graph->mirrors_set[0]);
898:     for (i = 1; i < graph->nvtxs; i++) graph->mirrors_set[i] = graph->mirrors_set[i - 1] + graph->mirrors[i - 1];

900:     /* fill arrays */
901:     PetscArrayzero(graph->mirrors, graph->nvtxs);
902:     for (j = 0; j < n_shared[0]; j++) {
903:       i = shared[0][j];
904:       if (graph->count[i] > 1) graph->mirrors_set[i][graph->mirrors[i]++] = global_indices[i];
905:     }
906:     PetscSortIntWithArray(graph->nvtxs, global_indices, local_indices);
907:     for (i = 0; i < graph->nvtxs; i++) {
908:       if (graph->mirrors[i] > 0) {
909:         PetscFindInt(graph->mirrors_set[i][0], graph->nvtxs, global_indices, &k);
910:         j = global_indices[k];
911:         while (k > 0 && global_indices[k - 1] == j) k--;
912:         for (j = 0; j < graph->mirrors[i]; j++) graph->mirrors_set[i][j] = local_indices[k + j];
913:         PetscSortInt(graph->mirrors[i], graph->mirrors_set[i]);
914:       }
915:     }
916:     PetscFree(local_indices);
917:     PetscFree(global_indices);
918:     ISDestroy(&to);
919:     ISDestroy(&from);
920:   }
921:   PetscArrayzero(graph->count, graph->nvtxs);

923:   /* Count total number of neigh per node */
924:   k = 0;
925:   for (i = 1; i < n_neigh; i++) {
926:     k += n_shared[i];
927:     for (j = 0; j < n_shared[i]; j++) graph->count[shared[i][j]] += 1;
928:   }
929:   /* Allocate space for storing the set of neighbours for each node */
930:   if (graph->nvtxs) PetscMalloc1(k, &graph->neighbours_set[0]);
931:   for (i = 1; i < graph->nvtxs; i++) { /* dont count myself */
932:     graph->neighbours_set[i] = graph->neighbours_set[i - 1] + graph->count[i - 1];
933:   }
934:   /* Get information for sharing subdomains */
935:   PetscArrayzero(graph->count, graph->nvtxs);
936:   for (i = 1; i < n_neigh; i++) { /* dont count myself */
937:     s = n_shared[i];
938:     for (j = 0; j < s; j++) {
939:       k                                         = shared[i][j];
940:       graph->neighbours_set[k][graph->count[k]] = neigh[i];
941:       graph->count[k] += 1;
942:     }
943:   }
944:   /* sort set of sharing subdomains */
945:   for (i = 0; i < graph->nvtxs; i++) PetscSortRemoveDupsInt(&graph->count[i], graph->neighbours_set[i]);
946:   /* free memory allocated by ISLocalToGlobalMappingGetInfo */
947:   ISLocalToGlobalMappingRestoreInfo(graph->l2gmap, &n_neigh, &neigh, &n_shared, &shared);

949:   /*
950:      Get info for dofs splitting
951:      User can specify just a subset; an additional field is considered as a complementary field
952:   */
953:   for (i = 0, k = 0; i < n_ISForDofs; i++) {
954:     PetscInt bs;

956:     ISGetBlockSize(ISForDofs[i], &bs);
957:     k += bs;
958:   }
959:   for (i = 0; i < graph->nvtxs; i++) graph->which_dof[i] = k; /* by default a dof belongs to the complement set */
960:   for (i = 0, k = 0; i < n_ISForDofs; i++) {
961:     PetscInt bs;

963:     ISGetLocalSize(ISForDofs[i], &is_size);
964:     ISGetBlockSize(ISForDofs[i], &bs);
965:     ISGetIndices(ISForDofs[i], (const PetscInt **)&is_indices);
966:     for (j = 0; j < is_size / bs; j++) {
967:       PetscInt b;

969:       for (b = 0; b < bs; b++) {
970:         PetscInt jj = bs * j + b;

972:         if (is_indices[jj] > -1 && is_indices[jj] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */
973:           graph->which_dof[is_indices[jj]] = k + b;
974:         }
975:       }
976:     }
977:     ISRestoreIndices(ISForDofs[i], (const PetscInt **)&is_indices);
978:     k += bs;
979:   }

981:   /* Take into account Neumann nodes */
982:   if (neumann_is) {
983:     ISGetLocalSize(neumann_is, &is_size);
984:     ISGetIndices(neumann_is, (const PetscInt **)&is_indices);
985:     for (i = 0; i < is_size; i++) {
986:       if (is_indices[i] > -1 && is_indices[i] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */
987:         graph->special_dof[is_indices[i]] = PCBDDCGRAPH_NEUMANN_MARK;
988:       }
989:     }
990:     ISRestoreIndices(neumann_is, (const PetscInt **)&is_indices);
991:   }
992:   /* Take into account Dirichlet nodes (they overwrite any neumann boundary mark previously set) */
993:   if (dirichlet_is) {
994:     ISGetLocalSize(dirichlet_is, &is_size);
995:     ISGetIndices(dirichlet_is, (const PetscInt **)&is_indices);
996:     for (i = 0; i < is_size; i++) {
997:       if (is_indices[i] > -1 && is_indices[i] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */
998:         if (commsize > graph->commsizelimit) {                  /* dirichlet nodes treated as internal */
999:           PetscBTSet(graph->touched, is_indices[i]);
1000:           graph->subset[is_indices[i]] = 0;
1001:         }
1002:         graph->special_dof[is_indices[i]] = PCBDDCGRAPH_DIRICHLET_MARK;
1003:       }
1004:     }
1005:     ISRestoreIndices(dirichlet_is, (const PetscInt **)&is_indices);
1006:   }
1007:   /* mark local periodic nodes (if any) and adapt CSR graph (if any) */
1008:   if (graph->mirrors) {
1009:     for (i = 0; i < graph->nvtxs; i++)
1010:       if (graph->mirrors[i]) graph->special_dof[i] = PCBDDCGRAPH_LOCAL_PERIODIC_MARK;

1012:     if (graph->xadj) {
1013:       PetscInt *new_xadj, *new_adjncy;
1014:       /* sort CSR graph */
1015:       for (i = 0; i < graph->nvtxs; i++) PetscSortInt(graph->xadj[i + 1] - graph->xadj[i], &graph->adjncy[graph->xadj[i]]);
1016:       /* adapt local CSR graph in case of local periodicity */
1017:       k = 0;
1018:       for (i = 0; i < graph->nvtxs; i++)
1019:         for (j = graph->xadj[i]; j < graph->xadj[i + 1]; j++) k += graph->mirrors[graph->adjncy[j]];

1021:       PetscMalloc1(graph->nvtxs + 1, &new_xadj);
1022:       PetscMalloc1(k + graph->xadj[graph->nvtxs], &new_adjncy);
1023:       new_xadj[0] = 0;
1024:       for (i = 0; i < graph->nvtxs; i++) {
1025:         k = graph->xadj[i + 1] - graph->xadj[i];
1026:         PetscArraycpy(&new_adjncy[new_xadj[i]], &graph->adjncy[graph->xadj[i]], k);
1027:         new_xadj[i + 1] = new_xadj[i] + k;
1028:         for (j = graph->xadj[i]; j < graph->xadj[i + 1]; j++) {
1029:           k = graph->mirrors[graph->adjncy[j]];
1030:           PetscArraycpy(&new_adjncy[new_xadj[i + 1]], graph->mirrors_set[graph->adjncy[j]], k);
1031:           new_xadj[i + 1] += k;
1032:         }
1033:         k = new_xadj[i + 1] - new_xadj[i];
1034:         PetscSortRemoveDupsInt(&k, &new_adjncy[new_xadj[i]]);
1035:         new_xadj[i + 1] = new_xadj[i] + k;
1036:       }
1037:       /* set new CSR into graph */
1038:       PetscFree(graph->xadj);
1039:       PetscFree(graph->adjncy);
1040:       graph->xadj   = new_xadj;
1041:       graph->adjncy = new_adjncy;
1042:     }
1043:   }

1045:   /* mark special nodes (if any) -> each will become a single node equivalence class */
1046:   if (custom_primal_vertices) {
1047:     ISGetLocalSize(custom_primal_vertices, &is_size);
1048:     ISGetIndices(custom_primal_vertices, (const PetscInt **)&is_indices);
1049:     for (i = 0, j = 0; i < is_size; i++) {
1050:       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 */
1051:         graph->special_dof[is_indices[i]] = PCBDDCGRAPH_SPECIAL_MARK - j;
1052:         j++;
1053:       }
1054:     }
1055:     ISRestoreIndices(custom_primal_vertices, (const PetscInt **)&is_indices);
1056:   }

1058:   /* mark interior nodes (if commsize > graph->commsizelimit) as touched and belonging to partition number 0 */
1059:   if (commsize > graph->commsizelimit) {
1060:     for (i = 0; i < graph->nvtxs; i++) {
1061:       if (!graph->count[i]) {
1062:         PetscBTSet(graph->touched, i);
1063:         graph->subset[i] = 0;
1064:       }
1065:     }
1066:   }

1068:   /* init graph structure and compute default subsets */
1069:   nodes_touched = 0;
1070:   for (i = 0; i < graph->nvtxs; i++) {
1071:     if (PetscBTLookup(graph->touched, i)) nodes_touched++;
1072:   }
1073:   i            = 0;
1074:   graph->ncc   = 0;
1075:   total_counts = 0;

1077:   /* allocated space for queues */
1078:   if (commsize == graph->commsizelimit) {
1079:     PetscMalloc2(graph->nvtxs + 1, &graph->cptr, graph->nvtxs, &graph->queue);
1080:   } else {
1081:     PetscInt nused = graph->nvtxs - nodes_touched;
1082:     PetscMalloc2(nused + 1, &graph->cptr, nused, &graph->queue);
1083:   }

1085:   while (nodes_touched < graph->nvtxs) {
1086:     /*  find first untouched node in local ordering */
1087:     while (PetscBTLookup(graph->touched, i)) i++;
1088:     PetscBTSet(graph->touched, i);
1089:     graph->subset[i]           = graph->ncc + 1;
1090:     graph->cptr[graph->ncc]    = total_counts;
1091:     graph->queue[total_counts] = i;
1092:     total_counts++;
1093:     nodes_touched++;
1094:     /* now find all other nodes having the same set of sharing subdomains */
1095:     for (j = i + 1; j < graph->nvtxs; j++) {
1096:       /* check for same number of sharing subdomains, dof number and same special mark */
1097:       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]) {
1098:         /* check for same set of sharing subdomains */
1099:         same_set = PETSC_TRUE;
1100:         for (k = 0; k < graph->count[j]; k++) {
1101:           if (graph->neighbours_set[i][k] != graph->neighbours_set[j][k]) same_set = PETSC_FALSE;
1102:         }
1103:         /* I have found a friend of mine */
1104:         if (same_set) {
1105:           PetscBTSet(graph->touched, j);
1106:           graph->subset[j] = graph->ncc + 1;
1107:           nodes_touched++;
1108:           graph->queue[total_counts] = j;
1109:           total_counts++;
1110:         }
1111:       }
1112:     }
1113:     graph->ncc++;
1114:   }
1115:   /* 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 */
1116:   graph->n_subsets = graph->ncc;
1117:   PetscMalloc1(graph->n_subsets, &graph->subset_ncc);
1118:   for (i = 0; i < graph->n_subsets; i++) graph->subset_ncc[i] = 1;
1119:   /* final pointer */
1120:   graph->cptr[graph->ncc] = total_counts;

1122:   /* For consistency reasons (among neighbours), I need to sort (by global ordering) each connected component */
1123:   /* Get a reference node (min index in global ordering) for each subset for tagging messages */
1124:   PetscMalloc1(graph->ncc, &graph->subset_ref_node);
1125:   PetscMalloc1(graph->cptr[graph->ncc], &queue_global);
1126:   ISLocalToGlobalMappingApply(graph->l2gmap, graph->cptr[graph->ncc], graph->queue, queue_global);
1127:   for (j = 0; j < graph->ncc; j++) {
1128:     PetscSortIntWithArray(graph->cptr[j + 1] - graph->cptr[j], &queue_global[graph->cptr[j]], &graph->queue[graph->cptr[j]]);
1129:     graph->subset_ref_node[j] = graph->queue[graph->cptr[j]];
1130:   }
1131:   PetscFree(queue_global);
1132:   graph->queue_sorted = PETSC_TRUE;

1134:   /* save information on subsets (needed when analyzing the connected components) */
1135:   if (graph->ncc) {
1136:     PetscMalloc2(graph->ncc, &graph->subset_size, graph->ncc, &graph->subset_idxs);
1137:     PetscMalloc1(graph->cptr[graph->ncc], &graph->subset_idxs[0]);
1138:     PetscArrayzero(graph->subset_idxs[0], graph->cptr[graph->ncc]);
1139:     for (j = 1; j < graph->ncc; j++) {
1140:       graph->subset_size[j - 1] = graph->cptr[j] - graph->cptr[j - 1];
1141:       graph->subset_idxs[j]     = graph->subset_idxs[j - 1] + graph->subset_size[j - 1];
1142:     }
1143:     graph->subset_size[graph->ncc - 1] = graph->cptr[graph->ncc] - graph->cptr[graph->ncc - 1];
1144:     PetscArraycpy(graph->subset_idxs[0], graph->queue, graph->cptr[graph->ncc]);
1145:   }

1147:   /* renumber reference nodes */
1148:   ISCreateGeneral(PetscObjectComm((PetscObject)(graph->l2gmap)), graph->ncc, graph->subset_ref_node, PETSC_COPY_VALUES, &subset_n);
1149:   ISLocalToGlobalMappingApplyIS(graph->l2gmap, subset_n, &subset);
1150:   ISDestroy(&subset_n);
1151:   ISRenumber(subset, NULL, NULL, &subset_n);
1152:   ISDestroy(&subset);
1153:   ISGetLocalSize(subset_n, &k);
1155:   ISGetIndices(subset_n, &is_indices);
1156:   PetscArraycpy(graph->subset_ref_node, is_indices, graph->ncc);
1157:   ISRestoreIndices(subset_n, &is_indices);
1158:   ISDestroy(&subset_n);

1160:   /* free workspace */
1161:   graph->setupcalled = PETSC_TRUE;
1162:   return 0;
1163: }

1165: PetscErrorCode PCBDDCGraphResetCoords(PCBDDCGraph graph)
1166: {
1167:   if (!graph) return 0;
1168:   PetscFree(graph->coords);
1169:   graph->cdim  = 0;
1170:   graph->cnloc = 0;
1171:   graph->cloc  = PETSC_FALSE;
1172:   return 0;
1173: }

1175: PetscErrorCode PCBDDCGraphResetCSR(PCBDDCGraph graph)
1176: {
1177:   if (!graph) return 0;
1178:   if (graph->freecsr) {
1179:     PetscFree(graph->xadj);
1180:     PetscFree(graph->adjncy);
1181:   } else {
1182:     graph->xadj   = NULL;
1183:     graph->adjncy = NULL;
1184:   }
1185:   graph->freecsr   = PETSC_FALSE;
1186:   graph->nvtxs_csr = 0;
1187:   return 0;
1188: }

1190: PetscErrorCode PCBDDCGraphReset(PCBDDCGraph graph)
1191: {
1192:   if (!graph) return 0;
1193:   ISLocalToGlobalMappingDestroy(&graph->l2gmap);
1194:   PetscFree(graph->subset_ncc);
1195:   PetscFree(graph->subset_ref_node);
1196:   if (graph->nvtxs) PetscFree(graph->neighbours_set[0]);
1197:   PetscBTDestroy(&graph->touched);
1198:   PetscFree5(graph->count, graph->neighbours_set, graph->subset, graph->which_dof, graph->special_dof);
1199:   PetscFree2(graph->cptr, graph->queue);
1200:   if (graph->mirrors) PetscFree(graph->mirrors_set[0]);
1201:   PetscFree2(graph->mirrors, graph->mirrors_set);
1202:   if (graph->subset_idxs) PetscFree(graph->subset_idxs[0]);
1203:   PetscFree2(graph->subset_size, graph->subset_idxs);
1204:   ISDestroy(&graph->dirdofs);
1205:   ISDestroy(&graph->dirdofsB);
1206:   if (graph->n_local_subs) PetscFree(graph->local_subs);
1207:   graph->has_dirichlet       = PETSC_FALSE;
1208:   graph->twodimset           = PETSC_FALSE;
1209:   graph->twodim              = PETSC_FALSE;
1210:   graph->nvtxs               = 0;
1211:   graph->nvtxs_global        = 0;
1212:   graph->n_subsets           = 0;
1213:   graph->custom_minimal_size = 1;
1214:   graph->n_local_subs        = 0;
1215:   graph->maxcount            = PETSC_MAX_INT;
1216:   graph->setupcalled         = PETSC_FALSE;
1217:   return 0;
1218: }

1220: PetscErrorCode PCBDDCGraphInit(PCBDDCGraph graph, ISLocalToGlobalMapping l2gmap, PetscInt N, PetscInt maxcount)
1221: {
1222:   PetscInt n;

1228:   /* raise an error if already allocated */
1230:   /* set number of vertices */
1231:   PetscObjectReference((PetscObject)l2gmap);
1232:   graph->l2gmap = l2gmap;
1233:   ISLocalToGlobalMappingGetSize(l2gmap, &n);
1234:   graph->nvtxs        = n;
1235:   graph->nvtxs_global = N;
1236:   /* allocate used space */
1237:   PetscBTCreate(graph->nvtxs, &graph->touched);
1238:   PetscMalloc5(graph->nvtxs, &graph->count, graph->nvtxs, &graph->neighbours_set, graph->nvtxs, &graph->subset, graph->nvtxs, &graph->which_dof, graph->nvtxs, &graph->special_dof);
1239:   /* zeroes memory */
1240:   PetscArrayzero(graph->count, graph->nvtxs);
1241:   PetscArrayzero(graph->subset, graph->nvtxs);
1242:   /* use -1 as a default value for which_dof array */
1243:   for (n = 0; n < graph->nvtxs; n++) graph->which_dof[n] = -1;
1244:   PetscArrayzero(graph->special_dof, graph->nvtxs);
1245:   /* zeroes first pointer to neighbour set */
1246:   if (graph->nvtxs) graph->neighbours_set[0] = NULL;
1247:   /* zeroes workspace for values of ncc */
1248:   graph->subset_ncc      = NULL;
1249:   graph->subset_ref_node = NULL;
1250:   /* maxcount for cc */
1251:   graph->maxcount = maxcount;
1252:   return 0;
1253: }

1255: PetscErrorCode PCBDDCGraphDestroy(PCBDDCGraph *graph)
1256: {
1257:   PCBDDCGraphResetCSR(*graph);
1258:   PCBDDCGraphResetCoords(*graph);
1259:   PCBDDCGraphReset(*graph);
1260:   PetscFree(*graph);
1261:   return 0;
1262: }

1264: PetscErrorCode PCBDDCGraphCreate(PCBDDCGraph *graph)
1265: {
1266:   PCBDDCGraph new_graph;

1268:   PetscNew(&new_graph);
1269:   new_graph->custom_minimal_size = 1;
1270:   new_graph->commsizelimit       = 1;
1271:   *graph                         = new_graph;
1272:   return 0;
1273: }