Actual source code: bddcgraph.c
1: #include <petsc/private/petscimpl.h>
2: #include <petsc/private/pcbddcprivateimpl.h>
3: #include <petsc/private/pcbddcstructsimpl.h>
4: #include <petsc/private/hashmapi.h>
5: #include <petsc/private/hashmapijk.h>
6: #include <petsc/private/hashsetij.h>
7: #include <petsc/private/pcbddcgraphhashmap.h>
8: #include <petscpartitioner.h>
9: #include <petscsf.h>
11: typedef enum {
12: PCBDDCGRAPH_COMPONENT_VERTEX,
13: PCBDDCGRAPH_COMPONENT_EDGE,
14: PCBDDCGRAPH_COMPONENT_FACE
15: } PCBDDCGraphComponentType;
17: typedef struct {
18: PetscInt repdof;
19: PetscInt ndofs;
20: PetscInt subset;
21: } PCBDDCGraphLocalEntity;
23: /* Unlike PCBDDCGraphNodeHash(), Dirichlet reconstruction must retain the local incidence of MPI-shared dofs. */
24: static inline PetscHash_t PCBDDCGraphDirichletNodeHash_Private(const PCBDDCGraphNode *node)
25: {
26: PetscHash_t hash;
28: hash = PetscHashCombine(PetscHashInt(node->count), PetscHashInt(node->which_dof));
29: for (PetscInt i = 0; i < node->count; i++) hash = PetscHashCombine(hash, PetscHashInt(node->neighbours_set[i]));
30: hash = PetscHashCombine(hash, PetscHashInt(node->local_groups_count));
31: for (PetscInt i = 0; i < node->local_groups_count; i++) hash = PetscHashCombine(hash, PetscHashInt(node->local_groups[i]));
32: return hash;
33: }
35: static inline int PCBDDCGraphDirichletNodeEqual_Private(const PCBDDCGraphNode *a, const PCBDDCGraphNode *b)
36: {
37: if (a->count != b->count || a->which_dof != b->which_dof || a->local_groups_count != b->local_groups_count) return 0;
38: for (PetscInt i = 0; i < a->count; i++)
39: if (a->neighbours_set[i] != b->neighbours_set[i]) return 0;
40: for (PetscInt i = 0; i < a->local_groups_count; i++)
41: if (a->local_groups[i] != b->local_groups[i]) return 0;
42: return 1;
43: }
45: PETSC_HASH_MAP(HMapPCBDDCDirichletNode, PCBDDCGraphNode *, PetscInt, PCBDDCGraphDirichletNodeHash_Private, PCBDDCGraphDirichletNodeEqual_Private, -1)
47: static PCBDDCGraphComponentType PCBDDCGraphClassifyEntity_Private(PCBDDCGraph graph, const PCBDDCGraphNode *node, PetscInt ndofs)
48: {
49: if (ndofs <= graph->custom_minimal_size || node->count > graph->maxcount) return PCBDDCGRAPH_COMPONENT_VERTEX;
50: if (!graph->twodim && node->count == 2 && node->special_dof != PCBDDCGRAPH_NEUMANN_MARK) return PCBDDCGRAPH_COMPONENT_FACE;
51: return PCBDDCGRAPH_COMPONENT_EDGE;
52: }
54: static PCBDDCGraphComponentType PCBDDCGraphGetComponentType_Private(PCBDDCGraph graph, PetscInt cc)
55: {
56: const PetscInt repdof = graph->queue[graph->cptr[cc]];
57: const PetscInt ccsize = graph->cptr[cc + 1] - graph->cptr[cc];
59: return PCBDDCGraphClassifyEntity_Private(graph, &graph->nodes[repdof], ccsize);
60: }
62: static PetscBool PCBDDCGraphNodeIsCanonical_Private(const PCBDDCGraphNode *node)
63: {
64: return (PetscBool)(node->local_groups_count > 1 && node->local_sub == node->local_groups[0]);
65: }
67: static PetscBool PCBDDCGraphIncidenceStrictSubset_Private(PetscInt na, const PetscInt a[], PetscInt nb, const PetscInt b[])
68: {
69: PetscInt ia = 0, ib = 0;
71: if (na >= nb) return PETSC_FALSE;
72: while (ia < na && ib < nb) {
73: if (a[ia] == b[ib]) {
74: ia++;
75: ib++;
76: } else if (a[ia] > b[ib]) ib++;
77: else return PETSC_FALSE;
78: }
79: return (PetscBool)(ia == na);
80: }
82: static PetscErrorCode PCBDDCGraphAddAdjacencyEvidence_Private(PetscHMapIJK evidence, PetscInt threshold, PetscHSetIJ adjacency)
83: {
84: PetscHashIter iter;
85: PetscHashIJKKey key;
86: PetscInt count;
88: PetscFunctionBegin;
89: PetscHashIterBegin(evidence, iter);
90: while (!PetscHashIterAtEnd(evidence, iter)) {
91: PetscHashIterGetKey(evidence, iter, key);
92: PetscHashIterGetVal(evidence, iter, count);
93: if (count >= threshold) {
94: PetscHashIJKey pair = {key.i, key.j};
96: PetscCall(PetscHSetIJAdd(adjacency, pair));
97: }
98: PetscHashIterNext(evidence, iter);
99: }
100: PetscFunctionReturn(PETSC_SUCCESS);
101: }
103: PetscErrorCode PCBDDCDestroyGraphCandidatesIS(PetscCtxRt ctx)
104: {
105: PCBDDCGraphCandidates cand = *(PCBDDCGraphCandidates *)ctx;
107: PetscFunctionBegin;
108: for (PetscInt i = 0; i < cand->nfc; i++) PetscCall(ISDestroy(&cand->Faces[i]));
109: for (PetscInt i = 0; i < cand->nec; i++) PetscCall(ISDestroy(&cand->Edges[i]));
110: PetscCall(PetscFree(cand->Faces));
111: PetscCall(PetscFree(cand->Edges));
112: PetscCall(ISDestroy(&cand->Vertices));
113: PetscCall(PetscFree(cand));
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: PetscErrorCode PCBDDCGraphGetDirichletDofsB(PCBDDCGraph graph, IS *dirdofs)
118: {
119: PetscFunctionBegin;
120: if (graph->dirdofsB) {
121: PetscCall(PetscObjectReference((PetscObject)graph->dirdofsB));
122: } else if (graph->has_dirichlet) {
123: PetscInt i, size;
124: PetscInt *dirdofs_idxs;
126: size = 0;
127: for (i = 0; i < graph->nvtxs; i++) {
128: if (graph->nodes[i].count > 1 && graph->nodes[i].special_dof == PCBDDCGRAPH_DIRICHLET_MARK) size++;
129: }
131: PetscCall(PetscMalloc1(size, &dirdofs_idxs));
132: size = 0;
133: for (i = 0; i < graph->nvtxs; i++) {
134: if (graph->nodes[i].count > 1 && graph->nodes[i].special_dof == PCBDDCGRAPH_DIRICHLET_MARK) dirdofs_idxs[size++] = i;
135: }
136: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size, dirdofs_idxs, PETSC_OWN_POINTER, &graph->dirdofsB));
137: PetscCall(PetscObjectReference((PetscObject)graph->dirdofsB));
138: }
139: *dirdofs = graph->dirdofsB;
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: PetscErrorCode PCBDDCGraphGetDirichletDofs(PCBDDCGraph graph, IS *dirdofs)
144: {
145: PetscFunctionBegin;
146: if (graph->dirdofs) {
147: PetscCall(PetscObjectReference((PetscObject)graph->dirdofs));
148: } else if (graph->has_dirichlet) {
149: PetscInt i, size;
150: PetscInt *dirdofs_idxs;
152: size = 0;
153: for (i = 0; i < graph->nvtxs; i++) {
154: if (graph->nodes[i].special_dof == PCBDDCGRAPH_DIRICHLET_MARK) size++;
155: }
157: PetscCall(PetscMalloc1(size, &dirdofs_idxs));
158: size = 0;
159: for (i = 0; i < graph->nvtxs; i++) {
160: if (graph->nodes[i].special_dof == PCBDDCGRAPH_DIRICHLET_MARK) dirdofs_idxs[size++] = i;
161: }
162: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)graph->l2gmap), size, dirdofs_idxs, PETSC_OWN_POINTER, &graph->dirdofs));
163: PetscCall(PetscObjectReference((PetscObject)graph->dirdofs));
164: }
165: *dirdofs = graph->dirdofs;
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: PetscErrorCode PCBDDCGraphASCIIView(PCBDDCGraph graph, PetscInt verbosity_level, PetscViewer viewer)
170: {
171: PetscInt i, j, tabs;
172: PetscInt *queue_in_global_numbering;
174: PetscFunctionBegin;
175: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(graph->seq_graph ? PETSC_COMM_SELF : PetscObjectComm((PetscObject)graph->l2gmap), &viewer));
176: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
177: PetscCall(PetscViewerASCIIGetTab(viewer, &tabs));
178: PetscCall(PetscViewerASCIIPrintf(viewer, "--------------------------------------------------\n"));
179: PetscCall(PetscViewerFlush(viewer));
180: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Local BDDC graph for subdomain %04d (seq %d)\n", PetscGlobalRank, graph->seq_graph));
181: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Number of vertices %" PetscInt_FMT "\n", graph->nvtxs));
182: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Number of local subdomains %" PetscInt_FMT "\n", graph->n_local_subs ? graph->n_local_subs : 1));
183: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Custom minimal size %" PetscInt_FMT "\n", graph->custom_minimal_size));
184: if (graph->maxcount != PETSC_INT_MAX) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Max count %" PetscInt_FMT "\n", graph->maxcount));
185: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Topological two dim? %s (set %s)\n", PetscBools[graph->twodim], PetscBools[graph->twodimset]));
186: if (verbosity_level > 2) {
187: for (i = 0; i < graph->nvtxs; i++) {
188: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%" PetscInt_FMT ":\n", i));
189: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " which_dof: %" PetscInt_FMT "\n", graph->nodes[i].which_dof));
190: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " special_dof: %" PetscInt_FMT "\n", graph->nodes[i].special_dof));
191: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " shared by: %" PetscInt_FMT "\n", graph->nodes[i].count));
192: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
193: if (graph->nodes[i].count) {
194: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " set of neighbours:"));
195: for (j = 0; j < graph->nodes[i].count; j++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, graph->nodes[i].neighbours_set[j]));
196: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
197: }
198: PetscCall(PetscViewerASCIISetTab(viewer, tabs));
199: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
200: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " number of local groups: %" PetscInt_FMT "\n", graph->nodes[i].local_groups_count));
201: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
202: if (graph->nodes[i].local_groups_count) {
203: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " groups:"));
204: for (j = 0; j < graph->nodes[i].local_groups_count; j++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, graph->nodes[i].local_groups[j]));
205: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
206: }
207: PetscCall(PetscViewerASCIISetTab(viewer, tabs));
208: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
210: if (verbosity_level > 3) {
211: if (graph->xadj) {
212: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " local adj list:"));
213: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
214: for (j = graph->xadj[i]; j < graph->xadj[i + 1]; j++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, graph->adjncy[j]));
215: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
216: PetscCall(PetscViewerASCIISetTab(viewer, tabs));
217: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
218: } else {
219: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " no adj info\n"));
220: }
221: }
222: if (graph->n_local_subs) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " local sub id: %" PetscInt_FMT "\n", graph->local_subs[i]));
223: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " interface subset id: %" PetscInt_FMT "\n", graph->nodes[i].subset));
224: if (graph->nodes[i].subset && graph->subset_ncc) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " ncc for subset: %" PetscInt_FMT "\n", graph->subset_ncc[graph->nodes[i].subset - 1]));
225: }
226: }
227: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Total number of connected components %" PetscInt_FMT "\n", graph->ncc));
228: PetscCall(PetscMalloc1(graph->cptr[graph->ncc], &queue_in_global_numbering));
229: PetscCall(ISLocalToGlobalMappingApply(graph->l2gmap, graph->cptr[graph->ncc], graph->queue, queue_in_global_numbering));
230: for (i = 0; i < graph->ncc; i++) {
231: PetscInt node_num = graph->queue[graph->cptr[i]];
232: PetscBool printcc = PETSC_FALSE;
233: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " cc %" PetscInt_FMT " (size %" PetscInt_FMT ", fid %" PetscInt_FMT ", neighs:", i, graph->cptr[i + 1] - graph->cptr[i], graph->nodes[node_num].which_dof));
234: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
235: for (j = 0; j < graph->nodes[node_num].count; j++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT, graph->nodes[node_num].neighbours_set[j]));
236: if (verbosity_level > 1) {
237: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "):"));
238: if (verbosity_level > 2 || graph->twodim || graph->nodes[node_num].count > 2 || (graph->nodes[node_num].count == 2 && graph->nodes[node_num].special_dof == PCBDDCGRAPH_NEUMANN_MARK)) printcc = PETSC_TRUE;
239: if (printcc) {
240: for (j = graph->cptr[i]; j < graph->cptr[i + 1]; j++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT " (%" PetscInt_FMT ")", graph->queue[j], queue_in_global_numbering[j]));
241: }
242: } else {
243: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, ")"));
244: }
245: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
246: PetscCall(PetscViewerASCIISetTab(viewer, tabs));
247: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
248: }
249: PetscCall(PetscFree(queue_in_global_numbering));
250: PetscCall(PetscViewerFlush(viewer));
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }
254: PetscErrorCode PCBDDCGraphRestoreCandidatesIS(PCBDDCGraph graph, PetscInt *n_faces, IS *FacesIS[], PetscInt *n_edges, IS *EdgesIS[], IS *VerticesIS)
255: {
256: PetscInt i;
257: PetscContainer gcand;
259: PetscFunctionBegin;
260: PetscCall(PetscObjectQuery((PetscObject)graph->l2gmap, "_PCBDDCGraphCandidatesIS", (PetscObject *)&gcand));
261: if (gcand) {
262: if (n_faces) *n_faces = 0;
263: if (n_edges) *n_edges = 0;
264: if (FacesIS) *FacesIS = NULL;
265: if (EdgesIS) *EdgesIS = NULL;
266: if (VerticesIS) *VerticesIS = NULL;
267: }
268: if (n_faces) {
269: if (FacesIS) {
270: for (i = 0; i < *n_faces; i++) PetscCall(ISDestroy(&(*FacesIS)[i]));
271: PetscCall(PetscFree(*FacesIS));
272: }
273: *n_faces = 0;
274: }
275: if (n_edges) {
276: if (EdgesIS) {
277: for (i = 0; i < *n_edges; i++) PetscCall(ISDestroy(&(*EdgesIS)[i]));
278: PetscCall(PetscFree(*EdgesIS));
279: }
280: *n_edges = 0;
281: }
282: if (VerticesIS) PetscCall(ISDestroy(VerticesIS));
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: PetscErrorCode PCBDDCGraphGetCandidatesIS(PCBDDCGraph graph, PetscInt *n_faces, IS *FacesIS[], PetscInt *n_edges, IS *EdgesIS[], IS *VerticesIS)
287: {
288: IS *ISForFaces, *ISForEdges, ISForVertices;
289: PetscInt i, nfc, nec, nvc, *idx, *mark;
290: PetscContainer gcand;
292: PetscFunctionBegin;
293: PetscCall(PetscObjectQuery((PetscObject)graph->l2gmap, "_PCBDDCGraphCandidatesIS", (PetscObject *)&gcand));
294: if (gcand) {
295: PCBDDCGraphCandidates cand;
297: PetscCall(PetscContainerGetPointer(gcand, &cand));
298: if (n_faces) *n_faces = cand->nfc;
299: if (FacesIS) *FacesIS = cand->Faces;
300: if (n_edges) *n_edges = cand->nec;
301: if (EdgesIS) *EdgesIS = cand->Edges;
302: if (VerticesIS) *VerticesIS = cand->Vertices;
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
305: PetscCall(PetscCalloc1(graph->ncc, &mark));
306: /* loop on ccs to evaluate number of faces, edges and vertices */
307: nfc = 0;
308: nec = 0;
309: nvc = 0;
310: for (i = 0; i < graph->ncc; i++) {
311: const PCBDDCGraphComponentType type = PCBDDCGraphGetComponentType_Private(graph, i);
313: if (type == PCBDDCGRAPH_COMPONENT_FACE) {
314: nfc++;
315: mark[i] = 2;
316: } else if (type == PCBDDCGRAPH_COMPONENT_EDGE) {
317: nec++;
318: mark[i] = 1;
319: } else {
320: nvc += graph->cptr[i + 1] - graph->cptr[i];
321: }
322: }
324: /* allocate IS arrays for faces, edges. Vertices need a single index set. */
325: if (FacesIS) PetscCall(PetscMalloc1(nfc, &ISForFaces));
326: if (EdgesIS) PetscCall(PetscMalloc1(nec, &ISForEdges));
327: if (VerticesIS) PetscCall(PetscMalloc1(nvc, &idx));
329: /* loop on ccs to compute index sets for faces and edges */
330: if (!graph->queue_sorted) {
331: PetscInt *queue_global;
333: PetscCall(PetscMalloc1(graph->cptr[graph->ncc], &queue_global));
334: PetscCall(ISLocalToGlobalMappingApply(graph->l2gmap, graph->cptr[graph->ncc], graph->queue, queue_global));
335: for (i = 0; i < graph->ncc; i++) PetscCall(PetscSortIntWithArray(graph->cptr[i + 1] - graph->cptr[i], &queue_global[graph->cptr[i]], &graph->queue[graph->cptr[i]]));
336: PetscCall(PetscFree(queue_global));
337: graph->queue_sorted = PETSC_TRUE;
338: }
339: nfc = 0;
340: nec = 0;
341: for (i = 0; i < graph->ncc; i++) {
342: if (mark[i] == 2) {
343: if (FacesIS) PetscCall(ISCreateGeneral(PETSC_COMM_SELF, graph->cptr[i + 1] - graph->cptr[i], &graph->queue[graph->cptr[i]], PETSC_USE_POINTER, &ISForFaces[nfc]));
344: nfc++;
345: } else if (mark[i] == 1) {
346: if (EdgesIS) PetscCall(ISCreateGeneral(PETSC_COMM_SELF, graph->cptr[i + 1] - graph->cptr[i], &graph->queue[graph->cptr[i]], PETSC_USE_POINTER, &ISForEdges[nec]));
347: nec++;
348: }
349: }
351: /* index set for vertices */
352: if (VerticesIS) {
353: nvc = 0;
354: for (i = 0; i < graph->ncc; i++) {
355: if (!mark[i]) {
356: for (PetscInt j = graph->cptr[i]; j < graph->cptr[i + 1]; j++) {
357: idx[nvc] = graph->queue[j];
358: nvc++;
359: }
360: }
361: }
362: /* sort vertex set (by local ordering) */
363: PetscCall(PetscSortInt(nvc, idx));
364: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nvc, idx, PETSC_OWN_POINTER, &ISForVertices));
365: }
366: PetscCall(PetscFree(mark));
368: /* get back info */
369: if (n_faces) *n_faces = nfc;
370: if (FacesIS) *FacesIS = ISForFaces;
371: if (n_edges) *n_edges = nec;
372: if (EdgesIS) *EdgesIS = ISForEdges;
373: if (VerticesIS) *VerticesIS = ISForVertices;
374: PetscFunctionReturn(PETSC_SUCCESS);
375: }
377: PetscErrorCode PCBDDCGraphCreateLocalSubdomainAdjacency(PCBDDCGraph graph, PetscInt *n_subs, PetscInt **xadj, PetscInt **adjncy)
378: {
379: PetscHMapPCBDDCDirichletNode dirichlet_nodes;
380: PetscHMapIJK edge_evidence, vertex_evidence, entity_pair_seen;
381: PetscHSetIJ adjacency;
382: PCBDDCGraphLocalEntity *entities;
383: PetscInt *cursor, *sub_entity_ptr, *sub_entities, *sub_entity_cursor;
384: PetscInt n_local_subs = graph->n_local_subs, n_entities, n_dirichlet_entities = 0;
385: PetscBool two_dimensional = graph->twodim;
387: PetscFunctionBegin;
388: PetscAssertPointer(n_subs, 2);
389: PetscAssertPointer(xadj, 3);
390: PetscAssertPointer(adjncy, 4);
391: PetscCheck(graph->setupcalled, PetscObjectComm((PetscObject)graph->l2gmap), PETSC_ERR_ORDER, "PCBDDCGraphSetUp() should be called first");
392: PetscCheck(graph->multi_element, PetscObjectComm((PetscObject)graph->l2gmap), PETSC_ERR_ARG_WRONGSTATE, "Local subdomain adjacency is only available for multi-element graphs");
393: PetscCheck(!graph->nvtxs || graph->local_subs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing local subdomain information");
394: *n_subs = n_local_subs;
396: /*
397: The dof connected components provide topological evidence for adjacency of local subdomains. Each component carries a sorted incidence set in
398: node->local_groups. Incidence to exactly two subdomains proves adjacency directly. For larger incidence sets, one field must provide enough
399: distinct lower-dimensional entities: one edge or two vertices in 2D, and three edges or vertices in 3D. The original subset number is used as
400: the topological-entity identity, so multiple dofs or components for the same entity are not counted twice. Dirichlet dofs are absent from the
401: component queues; their equivalence classes are reconstructed first and then treated by the same rules.
402: */
403: PetscCall(PetscHMapIJKCreate(&edge_evidence));
404: PetscCall(PetscHMapIJKCreate(&vertex_evidence));
405: PetscCall(PetscHMapIJKCreate(&entity_pair_seen));
406: PetscCall(PetscHSetIJCreate(&adjacency));
407: PetscCall(PetscMalloc1(graph->ncc + graph->nvtxs, &entities));
408: for (PetscInt cc = 0; cc < graph->ncc; cc++) {
409: const PetscInt repdof = graph->queue[graph->cptr[cc]];
411: entities[cc].repdof = repdof;
412: entities[cc].ndofs = graph->cptr[cc + 1] - graph->cptr[cc];
413: entities[cc].subset = graph->nodes[repdof].subset;
414: }
415: PetscCall(PetscHMapPCBDDCDirichletNodeCreate(&dirichlet_nodes));
416: for (PetscInt i = 0; i < graph->nvtxs; i++) {
417: PCBDDCGraphNode *node = &graph->nodes[i];
418: PetscHashIter iter;
419: PetscInt entity;
420: PetscBool missing;
422: if (node->special_dof != PCBDDCGRAPH_DIRICHLET_MARK || !PCBDDCGraphNodeIsCanonical_Private(node)) continue;
423: PetscCall(PetscHMapPCBDDCDirichletNodePut(dirichlet_nodes, node, &iter, &missing));
424: if (missing) {
425: entity = graph->ncc + n_dirichlet_entities++;
426: PetscCall(PetscHMapPCBDDCDirichletNodeIterSet(dirichlet_nodes, iter, entity));
427: entities[entity].repdof = i;
428: entities[entity].ndofs = 0;
429: entities[entity].subset = graph->n_subsets + n_dirichlet_entities;
430: } else PetscCall(PetscHMapPCBDDCDirichletNodeIterGet(dirichlet_nodes, iter, &entity));
431: entities[entity].ndofs++;
432: }
433: PetscCall(PetscHMapPCBDDCDirichletNodeDestroy(&dirichlet_nodes));
434: n_entities = graph->ncc + n_dirichlet_entities;
436: /* Invert entity incidence for the containment and Dirichlet corrections below. */
437: PetscCall(PetscCalloc1(n_local_subs + 1, &sub_entity_ptr));
438: for (PetscInt entity = 0; entity < n_entities; entity++) {
439: const PetscInt repdof = entities[entity].repdof;
440: const PCBDDCGraphNode *node = &graph->nodes[repdof];
442: if (!PCBDDCGraphNodeIsCanonical_Private(node)) continue;
443: for (PetscInt i = 0; i < node->local_groups_count; i++) {
444: const PetscInt group = node->local_groups[i];
446: PetscCheck(group >= 0 && group < n_local_subs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid local subdomain index %" PetscInt_FMT " not in [0,%" PetscInt_FMT ")", group, n_local_subs);
447: sub_entity_ptr[group + 1]++;
448: }
449: }
450: for (PetscInt i = 0; i < n_local_subs; i++) sub_entity_ptr[i + 1] += sub_entity_ptr[i];
451: PetscCall(PetscMalloc2(sub_entity_ptr[n_local_subs], &sub_entities, n_local_subs, &sub_entity_cursor));
452: PetscCall(PetscArraycpy(sub_entity_cursor, sub_entity_ptr, n_local_subs));
453: for (PetscInt entity = 0; entity < n_entities; entity++) {
454: const PetscInt repdof = entities[entity].repdof;
455: const PCBDDCGraphNode *node = &graph->nodes[repdof];
457: if (!PCBDDCGraphNodeIsCanonical_Private(node)) continue;
458: for (PetscInt i = 0; i < node->local_groups_count; i++) sub_entities[sub_entity_cursor[node->local_groups[i]]++] = entity;
459: }
460: /* Fold a Dirichlet fragment back into an unconstrained component with the same field and incidence. */
461: for (PetscInt entity = graph->ncc; entity < n_entities; entity++) {
462: const PCBDDCGraphNode *node = &graph->nodes[entities[entity].repdof];
463: const PetscInt *groups = node->local_groups;
464: const PetscInt ng = node->local_groups_count;
466: for (PetscInt p = sub_entity_ptr[groups[0]]; p < sub_entity_ptr[groups[0] + 1]; p++) {
467: const PetscInt other = sub_entities[p];
468: const PCBDDCGraphNode *onode;
469: PetscBool same;
471: if (other >= graph->ncc) continue;
472: onode = &graph->nodes[entities[other].repdof];
473: if (node->which_dof != onode->which_dof || ng != onode->local_groups_count) continue;
474: PetscCall(PetscArraycmp(groups, onode->local_groups, ng, &same));
475: if (same) {
476: entities[other].ndofs += entities[entity].ndofs;
477: entities[entity].subset = entities[other].subset;
478: entities[entity].ndofs = 0;
479: break;
480: }
481: }
482: }
483: for (PetscInt entity = 0; entity < n_entities; entity++) {
484: const PetscInt repdof = entities[entity].repdof;
485: const PCBDDCGraphNode *node = &graph->nodes[repdof];
486: PCBDDCGraphComponentType type = PCBDDCGraphClassifyEntity_Private(graph, node, entities[entity].ndofs);
487: const PetscInt *groups = node->local_groups;
488: const PetscInt ng = node->local_groups_count;
490: if (!entities[entity].ndofs || !PCBDDCGraphNodeIsCanonical_Private(node)) continue;
491: if (type == PCBDDCGRAPH_COMPONENT_VERTEX && ng > 2) {
492: PetscInt container_subset = -1;
494: /* A one-dof high-order edge looks like a vertex by size; its incidence lies in two distinct endpoint subsets. */
495: for (PetscInt p = sub_entity_ptr[groups[0]]; p < sub_entity_ptr[groups[0] + 1]; p++) {
496: const PetscInt other = sub_entities[p];
497: const PetscInt orepdof = entities[other].repdof;
498: const PCBDDCGraphNode *onode = &graph->nodes[orepdof];
500: if (!entities[other].ndofs) continue;
501: if (node->which_dof == onode->which_dof && PCBDDCGraphIncidenceStrictSubset_Private(ng, groups, onode->local_groups_count, onode->local_groups)) {
502: if (container_subset < 0) container_subset = entities[other].subset;
503: else if (entities[other].subset != container_subset) {
504: type = PCBDDCGRAPH_COMPONENT_EDGE;
505: break;
506: }
507: }
508: }
509: }
510: for (PetscInt i = 0; i < ng; i++) {
511: for (PetscInt j = i + 1; j < ng; j++) {
512: PetscHashIJKey pair = {groups[i], groups[j]};
514: if (ng == 2) PetscCall(PetscHSetIJAdd(adjacency, pair));
515: else {
516: PetscHMapIJK evidence;
517: PetscHashIJKKey entity_pair = {groups[i], groups[j], entities[entity].subset};
518: PetscHashIJKKey field_pair = {groups[i], groups[j], node->which_dof};
519: PetscHashIter iter;
520: PetscInt count;
521: PetscBool missing;
523: PetscCall(PetscHMapIJKPut(entity_pair_seen, entity_pair, &iter, &missing));
524: if (!missing) continue;
525: evidence = type == PCBDDCGRAPH_COMPONENT_EDGE ? edge_evidence : vertex_evidence;
526: PetscCall(PetscHMapIJKGetWithDefault(evidence, field_pair, 0, &count));
527: PetscCall(PetscHMapIJKSet(evidence, field_pair, count + 1));
528: }
529: }
530: }
531: }
532: PetscCall(PetscFree2(sub_entities, sub_entity_cursor));
533: PetscCall(PetscFree(sub_entity_ptr));
534: PetscCall(PetscFree(entities));
536: /* With only vertex dofs, the usual component-size test cannot distinguish 2D from 3D. */
537: if (two_dimensional) {
538: PetscHashIter iter;
540: PetscHashIterBegin(vertex_evidence, iter);
541: while (!PetscHashIterAtEnd(vertex_evidence, iter)) {
542: PetscInt count;
544: PetscHashIterGetVal(vertex_evidence, iter, count);
545: if (count >= 3) {
546: two_dimensional = PETSC_FALSE;
547: break;
548: }
549: PetscHashIterNext(vertex_evidence, iter);
550: }
551: }
552: PetscCall(PCBDDCGraphAddAdjacencyEvidence_Private(edge_evidence, two_dimensional ? 1 : 3, adjacency));
553: PetscCall(PCBDDCGraphAddAdjacencyEvidence_Private(vertex_evidence, two_dimensional ? 2 : 3, adjacency));
554: PetscCall(PetscHMapIJKDestroy(&edge_evidence));
555: PetscCall(PetscHMapIJKDestroy(&vertex_evidence));
556: PetscCall(PetscHMapIJKDestroy(&entity_pair_seen));
558: /* Convert the undirected pair set to symmetric CSR. */
559: PetscCall(PetscCalloc1(n_local_subs + 1, xadj));
560: PetscCall(PetscMalloc1(n_local_subs, &cursor));
561: {
562: PetscHashIter iter;
563: PetscHashIJKey key;
565: PetscHashIterBegin(adjacency, iter);
566: while (!PetscHashIterAtEnd(adjacency, iter)) {
567: PetscHashIterGetKey(adjacency, iter, key);
568: (*xadj)[key.i + 1]++;
569: (*xadj)[key.j + 1]++;
570: PetscHashIterNext(adjacency, iter);
571: }
572: }
573: for (PetscInt i = 0; i < n_local_subs; i++) (*xadj)[i + 1] += (*xadj)[i];
574: PetscCall(PetscMalloc1((*xadj)[n_local_subs], adjncy));
575: PetscCall(PetscArraycpy(cursor, *xadj, n_local_subs));
576: {
577: PetscHashIter iter;
578: PetscHashIJKey key;
580: PetscHashIterBegin(adjacency, iter);
581: while (!PetscHashIterAtEnd(adjacency, iter)) {
582: PetscHashIterGetKey(adjacency, iter, key);
583: (*adjncy)[cursor[key.i]++] = key.j;
584: (*adjncy)[cursor[key.j]++] = key.i;
585: PetscHashIterNext(adjacency, iter);
586: }
587: }
588: PetscCall(PetscHSetIJDestroy(&adjacency));
589: PetscCall(PetscFree(cursor));
590: for (PetscInt i = 0; i < n_local_subs; i++) PetscCall(PetscSortInt((*xadj)[i + 1] - (*xadj)[i], PetscSafePointerPlusOffset(*adjncy, (*xadj)[i])));
591: PetscFunctionReturn(PETSC_SUCCESS);
592: }
594: PetscErrorCode PCBDDCGraphPartitionLocalSubdomains(PCBDDCGraph graph, const char *prefix, PetscInt *nparts, PetscInt n_subs, PetscInt *xadj, PetscInt *adjncy, IS *partition)
595: {
596: PetscPartitioner part;
597: PetscSection part_section;
598: IS point_partition;
599: const PetscInt *points;
600: PetscInt *labels;
602: PetscFunctionBegin;
603: PetscAssertPointer(prefix, 2);
604: PetscAssertPointer(nparts, 3);
605: PetscAssertPointer(partition, 7);
606: PetscCheck(n_subs >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of local subdomains must be non-negative, got %" PetscInt_FMT, n_subs);
607: PetscCheck(n_subs == graph->n_local_subs, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of local subdomains %" PetscInt_FMT " does not match graph value %" PetscInt_FMT, n_subs, graph->n_local_subs);
608: if (!n_subs) {
609: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_COPY_VALUES, partition));
610: *nparts = 0;
611: PetscFunctionReturn(PETSC_SUCCESS);
612: }
613: PetscAssertPointer(xadj, 5);
614: if (xadj[n_subs]) PetscAssertPointer(adjncy, 6);
615: *nparts = PetscMin(PetscMax(*nparts, 1), n_subs);
616: if (*nparts == n_subs) {
617: PetscCall(ISCreateStride(PETSC_COMM_SELF, n_subs, 0, 1, partition));
618: PetscFunctionReturn(PETSC_SUCCESS);
619: }
621: PetscCall(PetscPartitionerCreate(PETSC_COMM_SELF, &part));
622: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)part, prefix));
623: PetscCall(PetscPartitionerSetFromOptions(part));
624: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &part_section));
625: PetscCall(PetscPartitionerPartition(part, *nparts, n_subs, xadj, adjncy, NULL, NULL, NULL, part_section, &point_partition));
626: PetscCall(PetscMalloc1(n_subs, &labels));
627: for (PetscInt i = 0; i < n_subs; i++) labels[i] = -1;
628: PetscCall(ISGetIndices(point_partition, &points));
629: for (PetscInt p = 0; p < *nparts; p++) {
630: PetscInt dof, offset;
632: PetscCall(PetscSectionGetDof(part_section, p, &dof));
633: PetscCall(PetscSectionGetOffset(part_section, p, &offset));
634: for (PetscInt i = 0; i < dof; i++) {
635: const PetscInt point = points[offset + i];
637: PetscCheck(point >= 0 && point < n_subs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid partition point %" PetscInt_FMT " not in [0,%" PetscInt_FMT ")", point, n_subs);
638: labels[point] = p;
639: }
640: }
641: PetscCall(ISRestoreIndices(point_partition, &points));
642: for (PetscInt i = 0; i < n_subs; i++) PetscCheck(labels[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Partition point %" PetscInt_FMT " was not assigned", i);
643: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n_subs, labels, PETSC_OWN_POINTER, partition));
644: PetscCall(ISDestroy(&point_partition));
645: PetscCall(PetscSectionDestroy(&part_section));
646: PetscCall(PetscPartitionerDestroy(&part));
647: PetscFunctionReturn(PETSC_SUCCESS);
648: }
650: PetscErrorCode PCBDDCGraphComputeConnectedComponents(PCBDDCGraph graph)
651: {
652: PetscBool adapt_interface;
653: MPI_Comm interface_comm;
654: PetscBT cornerp = NULL;
656: PetscFunctionBegin;
657: PetscCall(PetscObjectGetComm((PetscObject)graph->l2gmap, &interface_comm));
658: /* compute connected components locally */
659: PetscCall(PCBDDCGraphComputeConnectedComponentsLocal(graph));
660: if (graph->seq_graph) PetscFunctionReturn(PETSC_SUCCESS);
662: if (graph->active_coords && !graph->multi_element) { /* face based corner selection XXX multi_element */
663: PetscBT excluded;
664: PetscReal *wdist;
665: PetscInt n_neigh, *neigh, *n_shared, **shared;
666: PetscInt maxc, ns;
668: PetscCall(PetscBTCreate(graph->nvtxs, &cornerp));
669: PetscCall(ISLocalToGlobalMappingGetInfo(graph->l2gmap, &n_neigh, &neigh, &n_shared, &shared));
670: for (ns = 1, maxc = 0; ns < n_neigh; ns++) maxc = PetscMax(maxc, n_shared[ns]);
671: PetscCall(PetscMalloc1(maxc * graph->cdim, &wdist));
672: PetscCall(PetscBTCreate(maxc, &excluded));
674: for (ns = 1; ns < n_neigh; ns++) { /* first proc is self */
675: PetscReal *anchor, mdist;
676: PetscInt fst, j, k, d, cdim = graph->cdim, n = n_shared[ns];
677: PetscInt point1, point2, point3, point4;
679: /* import coordinates on shared interface */
680: PetscCall(PetscBTMemzero(n, excluded));
681: for (j = 0, fst = -1, k = 0; j < n; j++) {
682: PetscBool skip = PETSC_FALSE;
683: for (d = 0; d < cdim; d++) {
684: PetscReal c = graph->coords[shared[ns][j] * cdim + d];
685: skip = (PetscBool)(skip || c == PETSC_MAX_REAL);
686: wdist[k++] = c;
687: }
688: if (skip) PetscCall(PetscBTSet(excluded, j));
689: else if (fst == -1) fst = j;
690: }
691: if (fst == -1) continue;
693: /* the dofs are sorted by global numbering, so each rank starts from the same id
694: and it will detect the same corners from the given set */
696: /* find the farthest point from the starting one */
697: anchor = wdist + fst * cdim;
698: mdist = -1.0;
699: point1 = fst;
700: for (j = fst; j < n; j++) {
701: PetscReal dist = 0.0;
703: if (PetscUnlikely(PetscBTLookup(excluded, j))) continue;
704: for (d = 0; d < cdim; d++) dist += (wdist[j * cdim + d] - anchor[d]) * (wdist[j * cdim + d] - anchor[d]);
705: if (dist > mdist) {
706: mdist = dist;
707: point1 = j;
708: }
709: }
711: /* find the farthest point from point1 */
712: anchor = wdist + point1 * cdim;
713: mdist = -1.0;
714: point2 = point1;
715: for (j = fst; j < n; j++) {
716: PetscReal dist = 0.0;
718: if (PetscUnlikely(PetscBTLookup(excluded, j))) continue;
719: for (d = 0; d < cdim; d++) dist += (wdist[j * cdim + d] - anchor[d]) * (wdist[j * cdim + d] - anchor[d]);
720: if (dist > mdist) {
721: mdist = dist;
722: point2 = j;
723: }
724: }
726: /* find the third point maximizing the triangle area */
727: point3 = point2;
728: if (cdim > 2) {
729: PetscReal a = 0.0;
731: for (d = 0; d < cdim; d++) a += (wdist[point1 * cdim + d] - wdist[point2 * cdim + d]) * (wdist[point1 * cdim + d] - wdist[point2 * cdim + d]);
732: a = PetscSqrtReal(a);
733: mdist = -1.0;
734: for (j = fst; j < n; j++) {
735: PetscReal area, b = 0.0, c = 0.0, s;
737: if (PetscUnlikely(PetscBTLookup(excluded, j))) continue;
738: for (d = 0; d < cdim; d++) {
739: b += (wdist[point1 * cdim + d] - wdist[j * cdim + d]) * (wdist[point1 * cdim + d] - wdist[j * cdim + d]);
740: c += (wdist[point2 * cdim + d] - wdist[j * cdim + d]) * (wdist[point2 * cdim + d] - wdist[j * cdim + d]);
741: }
742: b = PetscSqrtReal(b);
743: c = PetscSqrtReal(c);
744: s = 0.5 * (a + b + c);
746: /* Heron's formula, area squared */
747: area = s * (s - a) * (s - b) * (s - c);
748: if (area > mdist) {
749: mdist = area;
750: point3 = j;
751: }
752: }
753: }
755: /* find the farthest point from point3 different from point1 and point2 */
756: anchor = wdist + point3 * cdim;
757: mdist = -1.0;
758: point4 = point3;
759: for (j = fst; j < n; j++) {
760: PetscReal dist = 0.0;
762: if (PetscUnlikely(PetscBTLookup(excluded, j)) || j == point1 || j == point2 || j == point3) continue;
763: for (d = 0; d < cdim; d++) dist += (wdist[j * cdim + d] - anchor[d]) * (wdist[j * cdim + d] - anchor[d]);
764: if (dist > mdist) {
765: mdist = dist;
766: point4 = j;
767: }
768: }
770: PetscCall(PetscBTSet(cornerp, shared[ns][point1]));
771: PetscCall(PetscBTSet(cornerp, shared[ns][point2]));
772: PetscCall(PetscBTSet(cornerp, shared[ns][point3]));
773: PetscCall(PetscBTSet(cornerp, shared[ns][point4]));
775: /* all dofs having the same coordinates will be primal */
776: for (j = fst; j < n; j++) {
777: PetscBool same[] = {PETSC_TRUE, PETSC_TRUE, PETSC_TRUE, PETSC_TRUE};
779: if (PetscUnlikely(PetscBTLookup(excluded, j))) continue;
780: for (d = 0; d < cdim; d++) {
781: same[0] = (PetscBool)(same[0] && (PetscAbsReal(wdist[j * cdim + d] - wdist[point1 * cdim + d]) < PETSC_SMALL));
782: same[1] = (PetscBool)(same[1] && (PetscAbsReal(wdist[j * cdim + d] - wdist[point2 * cdim + d]) < PETSC_SMALL));
783: same[2] = (PetscBool)(same[2] && (PetscAbsReal(wdist[j * cdim + d] - wdist[point3 * cdim + d]) < PETSC_SMALL));
784: same[3] = (PetscBool)(same[3] && (PetscAbsReal(wdist[j * cdim + d] - wdist[point4 * cdim + d]) < PETSC_SMALL));
785: }
786: if (same[0] || same[1] || same[2] || same[3]) PetscCall(PetscBTSet(cornerp, shared[ns][j]));
787: }
788: }
789: PetscCall(PetscBTDestroy(&excluded));
790: PetscCall(PetscFree(wdist));
791: PetscCall(ISLocalToGlobalMappingRestoreInfo(graph->l2gmap, &n_neigh, &neigh, &n_shared, &shared));
792: }
794: /* Adapt connected components if needed */
795: adapt_interface = (cornerp || graph->multi_element) ? PETSC_TRUE : PETSC_FALSE;
796: for (PetscInt i = 0; i < graph->n_subsets && !adapt_interface; i++) {
797: if (graph->subset_ncc[i] > 1) adapt_interface = PETSC_TRUE;
798: }
799: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &adapt_interface, 1, MPI_C_BOOL, MPI_LOR, interface_comm));
800: if (adapt_interface) {
801: PetscSF msf;
802: const PetscInt *n_ref_sharing;
803: PetscInt *labels, *rootlabels, *mrlabels;
804: PetscInt nr, nmr, nrs, ncc, cum_queue;
806: PetscCall(PetscMalloc1(graph->nvtxs, &labels));
807: PetscCall(PetscArrayzero(labels, graph->nvtxs));
808: for (PetscInt i = 0, k = 0; i < graph->ncc; i++) {
809: PetscInt s = 1;
810: for (PetscInt j = graph->cptr[i]; j < graph->cptr[i + 1]; j++) {
811: if (cornerp && PetscBTLookup(cornerp, graph->queue[j])) {
812: labels[graph->queue[j]] = -(k + s + 1);
813: s += 1;
814: } else {
815: labels[graph->queue[j]] = -(k + 1);
816: }
817: }
818: k += s;
819: }
820: PetscCall(PetscSFGetGraph(graph->interface_ref_sf, &nr, NULL, NULL, NULL));
821: PetscCall(PetscSFGetGraph(graph->interface_subset_sf, &nrs, NULL, NULL, NULL));
822: PetscCall(PetscSFGetMultiSF(graph->interface_subset_sf, &msf));
823: PetscCall(PetscSFGetGraph(msf, &nmr, NULL, NULL, NULL));
824: PetscCall(PetscCalloc2(nmr, &mrlabels, nrs, &rootlabels));
826: PetscCall(PetscSFComputeDegreeBegin(graph->interface_subset_sf, &n_ref_sharing));
827: PetscCall(PetscSFComputeDegreeEnd(graph->interface_subset_sf, &n_ref_sharing));
828: PetscCall(PetscSFGatherBegin(graph->interface_subset_sf, MPIU_INT, labels, mrlabels));
829: PetscCall(PetscSFGatherEnd(graph->interface_subset_sf, MPIU_INT, labels, mrlabels));
831: /* analyze contributions from processes
832: The structure of mrlabels is suitable to find intersections of ccs.
833: supposing the root subset has dimension 5 and leaves with labels:
834: 0: [4 4 7 4 7], (2 connected components)
835: 1: [3 2 2 3 2], (2 connected components)
836: 2: [1 1 6 5 6], (3 connected components)
837: the multiroot data and the new labels corresponding to intersected connected components will be (column major)
839: 4 4 7 4 7
840: mrlabels 3 2 2 3 2
841: 1 1 6 5 6
842: ---------
843: rootlabels 0 1 2 3 2
844: */
845: for (PetscInt i = 0, rcumlabels = 0, mcumlabels = 0; i < nr; i++) {
846: const PetscInt subset_size = graph->interface_ref_rsize[i];
847: const PetscInt *n_sharing = n_ref_sharing + rcumlabels;
848: const PetscInt *mrbuffer = mrlabels + mcumlabels;
849: PetscInt *rbuffer = rootlabels + rcumlabels;
850: PetscInt subset_counter = 0;
852: for (PetscInt j = 0; j < subset_size; j++) {
853: if (!rbuffer[j]) { /* found a new cc */
854: const PetscInt *jlabels = mrbuffer + j * n_sharing[0];
855: rbuffer[j] = ++subset_counter;
857: for (PetscInt k = j + 1; k < subset_size; k++) { /* check for other nodes in new cc */
858: PetscBool same_set = PETSC_TRUE;
859: const PetscInt *klabels = mrbuffer + k * n_sharing[0];
861: for (PetscInt s = 0; s < n_sharing[0]; s++) {
862: if (jlabels[s] != klabels[s]) {
863: same_set = PETSC_FALSE;
864: break;
865: }
866: }
867: if (same_set) rbuffer[k] = subset_counter;
868: }
869: }
870: }
871: if (subset_size) {
872: rcumlabels += subset_size;
873: mcumlabels += n_sharing[0] * subset_size;
874: }
875: }
877: /* Now communicate the intersected labels */
878: PetscCall(PetscSFBcastBegin(graph->interface_subset_sf, MPIU_INT, rootlabels, labels, MPI_REPLACE));
879: PetscCall(PetscSFBcastEnd(graph->interface_subset_sf, MPIU_INT, rootlabels, labels, MPI_REPLACE));
880: PetscCall(PetscFree2(mrlabels, rootlabels));
882: /* and adapt local connected components */
883: PetscInt *ocptr, *oqueue;
884: PetscBool *touched;
886: PetscCall(PetscMalloc3(graph->ncc + 1, &ocptr, graph->cptr[graph->ncc], &oqueue, graph->cptr[graph->ncc], &touched));
887: PetscCall(PetscArraycpy(ocptr, graph->cptr, graph->ncc + 1));
888: PetscCall(PetscArraycpy(oqueue, graph->queue, graph->cptr[graph->ncc]));
889: PetscCall(PetscArrayzero(touched, graph->cptr[graph->ncc]));
891: ncc = 0;
892: cum_queue = 0;
893: for (PetscInt i = 0; i < graph->ncc; i++) {
894: for (PetscInt j = ocptr[i]; j < ocptr[i + 1]; j++) {
895: const PetscInt jlabel = labels[oqueue[j]];
897: if (jlabel) {
898: graph->cptr[ncc] = cum_queue;
899: ncc++;
900: for (PetscInt k = j; k < ocptr[i + 1]; k++) { /* check for other nodes in new cc */
901: if (labels[oqueue[k]] == jlabel) {
902: graph->queue[cum_queue++] = oqueue[k];
903: labels[oqueue[k]] = 0;
904: }
905: }
906: }
907: }
908: }
909: PetscCall(PetscFree3(ocptr, oqueue, touched));
910: PetscCall(PetscFree(labels));
911: graph->cptr[ncc] = cum_queue;
912: graph->queue_sorted = PETSC_FALSE;
913: graph->ncc = ncc;
914: }
915: PetscCall(PetscBTDestroy(&cornerp));
917: /* Determine if we are in 2D or 3D */
918: if (!graph->twodimset) {
919: graph->twodim = PETSC_TRUE;
920: for (PetscInt i = 0; i < graph->ncc; i++) {
921: PetscInt repdof = graph->queue[graph->cptr[i]];
922: PetscInt ccsize = graph->cptr[i + 1] - graph->cptr[i];
923: if (graph->nodes[repdof].count > 2 && ccsize > graph->custom_minimal_size) {
924: graph->twodim = PETSC_FALSE;
925: break;
926: }
927: }
928: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &graph->twodim, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)graph->l2gmap)));
929: graph->twodimset = PETSC_TRUE;
930: }
931: PetscFunctionReturn(PETSC_SUCCESS);
932: }
934: static inline PetscErrorCode PCBDDCGraphComputeCC_Private(PCBDDCGraph graph, PetscInt pid, PetscInt *PETSC_RESTRICT queue_tip, PetscInt n_prev, PetscInt *n_added)
935: {
936: PetscInt i, j, n = 0;
938: const PetscInt *PETSC_RESTRICT xadj = graph->xadj;
939: const PetscInt *PETSC_RESTRICT adjncy = graph->adjncy;
940: const PetscInt *PETSC_RESTRICT subset_idxs = graph->subset_idxs[pid - 1];
941: const PetscInt *PETSC_RESTRICT local_subs = graph->local_subs;
942: const PetscInt subset_size = graph->subset_size[pid - 1];
944: PCBDDCGraphNode *PETSC_RESTRICT nodes = graph->nodes;
946: const PetscBool havecsr = (PetscBool)(!!xadj);
947: const PetscBool havesubs = (PetscBool)(!!graph->n_local_subs);
949: PetscFunctionBegin;
950: if (havecsr && !havesubs) {
951: for (i = -n_prev; i < 0; i++) {
952: const PetscInt start_dof = queue_tip[i];
954: /* we assume that if a dof has a size 1 adjacency list and the corresponding entry is negative, it is connected to all dofs */
955: if (xadj[start_dof + 1] - xadj[start_dof] == 1 && adjncy[xadj[start_dof]] < 0) {
956: for (j = 0; j < subset_size; j++) { /* pid \in [1,graph->n_subsets] */
957: const PetscInt dof = subset_idxs[j];
959: if (!nodes[dof].touched && nodes[dof].subset == pid) {
960: nodes[dof].touched = PETSC_TRUE;
961: queue_tip[n] = dof;
962: n++;
963: }
964: }
965: } else {
966: for (j = xadj[start_dof]; j < xadj[start_dof + 1]; j++) {
967: const PetscInt dof = adjncy[j];
969: if (!nodes[dof].touched && nodes[dof].subset == pid) {
970: nodes[dof].touched = PETSC_TRUE;
971: queue_tip[n] = dof;
972: n++;
973: }
974: }
975: }
976: }
977: } else if (havecsr && havesubs) {
978: const PetscInt sid = local_subs[queue_tip[-n_prev]];
980: for (i = -n_prev; i < 0; i++) {
981: const PetscInt start_dof = queue_tip[i];
983: /* 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 */
984: if (xadj[start_dof + 1] - xadj[start_dof] == 1 && adjncy[xadj[start_dof]] < 0) {
985: for (j = 0; j < subset_size; j++) { /* pid \in [1,graph->n_subsets] */
986: const PetscInt dof = subset_idxs[j];
988: if (!nodes[dof].touched && nodes[dof].subset == pid && local_subs[dof] == sid) {
989: nodes[dof].touched = PETSC_TRUE;
990: queue_tip[n] = dof;
991: n++;
992: }
993: }
994: } else {
995: for (j = xadj[start_dof]; j < xadj[start_dof + 1]; j++) {
996: const PetscInt dof = adjncy[j];
998: if (!nodes[dof].touched && nodes[dof].subset == pid && local_subs[dof] == sid) {
999: nodes[dof].touched = PETSC_TRUE;
1000: queue_tip[n] = dof;
1001: n++;
1002: }
1003: }
1004: }
1005: }
1006: } else if (havesubs) { /* sub info only */
1007: const PetscInt sid = local_subs[queue_tip[-n_prev]];
1009: for (j = 0; j < subset_size; j++) { /* pid \in [1,graph->n_subsets] */
1010: const PetscInt dof = subset_idxs[j];
1012: if (!nodes[dof].touched && nodes[dof].subset == pid && local_subs[dof] == sid) {
1013: nodes[dof].touched = PETSC_TRUE;
1014: queue_tip[n] = dof;
1015: n++;
1016: }
1017: }
1018: } else {
1019: for (j = 0; j < subset_size; j++) { /* pid \in [1,graph->n_subsets] */
1020: const PetscInt dof = subset_idxs[j];
1022: if (!nodes[dof].touched && nodes[dof].subset == pid) {
1023: nodes[dof].touched = PETSC_TRUE;
1024: queue_tip[n] = dof;
1025: n++;
1026: }
1027: }
1028: }
1029: *n_added = n;
1030: PetscFunctionReturn(PETSC_SUCCESS);
1031: }
1033: PetscErrorCode PCBDDCGraphComputeConnectedComponentsLocal(PCBDDCGraph graph)
1034: {
1035: PetscInt ncc, cum_queue;
1037: PetscFunctionBegin;
1038: PetscCheck(graph->setupcalled, PetscObjectComm((PetscObject)graph->l2gmap), PETSC_ERR_ORDER, "PCBDDCGraphSetUp should be called first");
1039: /* quiet return if there isn't any local info */
1040: if (!graph->xadj && !graph->n_local_subs) PetscFunctionReturn(PETSC_SUCCESS);
1042: /* reset any previous search of connected components */
1043: for (PetscInt i = 0; i < graph->nvtxs; i++) graph->nodes[i].touched = PETSC_FALSE;
1044: if (!graph->seq_graph) {
1045: for (PetscInt i = 0; i < graph->nvtxs; i++) {
1046: if (graph->nodes[i].special_dof == PCBDDCGRAPH_DIRICHLET_MARK || graph->nodes[i].count < 2) graph->nodes[i].touched = PETSC_TRUE;
1047: }
1048: }
1050: /* begin search for connected components */
1051: cum_queue = 0;
1052: ncc = 0;
1053: for (PetscInt n = 0; n < graph->n_subsets; n++) {
1054: const PetscInt *subset_idxs = graph->subset_idxs[n];
1055: const PetscInt pid = n + 1; /* partition labeled by 0 is discarded */
1057: PetscInt found = 0, prev = 0, first = 0, ncc_pid = 0;
1059: while (found != graph->subset_size[n]) {
1060: PetscInt added = 0;
1062: if (!prev) { /* search for new starting dof */
1063: while (graph->nodes[subset_idxs[first]].touched) first++;
1064: graph->nodes[subset_idxs[first]].touched = PETSC_TRUE;
1065: graph->queue[cum_queue] = subset_idxs[first];
1066: graph->cptr[ncc] = cum_queue;
1067: prev = 1;
1068: cum_queue++;
1069: found++;
1070: ncc_pid++;
1071: ncc++;
1072: }
1073: PetscCall(PCBDDCGraphComputeCC_Private(graph, pid, graph->queue + cum_queue, prev, &added));
1074: if (!added) {
1075: graph->subset_ncc[n] = ncc_pid;
1076: graph->cptr[ncc] = cum_queue;
1077: }
1078: prev = added;
1079: found += added;
1080: cum_queue += added;
1081: if (added && found == graph->subset_size[n]) {
1082: graph->subset_ncc[n] = ncc_pid;
1083: graph->cptr[ncc] = cum_queue;
1084: }
1085: }
1086: }
1087: graph->ncc = ncc;
1088: graph->queue_sorted = PETSC_FALSE;
1089: PetscFunctionReturn(PETSC_SUCCESS);
1090: }
1092: PetscErrorCode PCBDDCGraphSetUp(PCBDDCGraph graph, PetscInt custom_minimal_size, IS neumann_is, IS dirichlet_is, PetscInt n_ISForDofs, IS ISForDofs[], IS custom_primal_vertices)
1093: {
1094: IS subset;
1095: MPI_Comm comm;
1096: PetscHMapPCBDDCGraphNode subsetmaps;
1097: const PetscInt *is_indices;
1098: PetscInt *queue_global, *nodecount, **nodeneighs, *subset_sizes;
1099: PetscInt i, j, k, nodes_touched, is_size, nvtxs = graph->nvtxs;
1100: PetscMPIInt size, rank;
1102: PetscFunctionBegin;
1104: if (neumann_is) {
1106: PetscCheckSameComm(graph->l2gmap, 1, neumann_is, 3);
1107: }
1108: graph->has_dirichlet = PETSC_FALSE;
1109: if (dirichlet_is) {
1111: PetscCheckSameComm(graph->l2gmap, 1, dirichlet_is, 4);
1112: graph->has_dirichlet = PETSC_TRUE;
1113: }
1115: for (i = 0; i < n_ISForDofs; i++) {
1117: PetscCheckSameComm(graph->l2gmap, 1, ISForDofs[i], 6);
1118: }
1119: if (custom_primal_vertices) {
1121: PetscCheckSameComm(graph->l2gmap, 1, custom_primal_vertices, 7);
1122: }
1123: for (i = 0; i < nvtxs; i++) graph->nodes[i].touched = PETSC_FALSE;
1125: PetscCall(PetscObjectGetComm((PetscObject)graph->l2gmap, &comm));
1126: PetscCallMPI(MPI_Comm_size(comm, &size));
1127: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1129: /* custom_minimal_size */
1130: graph->custom_minimal_size = custom_minimal_size;
1132: /* get node info from l2gmap */
1133: PetscCall(ISLocalToGlobalMappingGetNodeInfo(graph->l2gmap, NULL, &nodecount, &nodeneighs));
1135: /* Allocate space for storing the set of neighbours for each node */
1136: graph->multi_element = PETSC_FALSE;
1137: for (i = 0; i < nvtxs; i++) {
1138: graph->nodes[i].count = nodecount[i];
1139: if (!graph->seq_graph) {
1140: PetscCall(PetscMalloc1(nodecount[i], &graph->nodes[i].neighbours_set));
1141: PetscCall(PetscArraycpy(graph->nodes[i].neighbours_set, nodeneighs[i], nodecount[i]));
1143: if (!graph->multi_element) {
1144: PetscInt nself;
1145: for (j = 0, nself = 0; j < graph->nodes[i].count; j++)
1146: if (graph->nodes[i].neighbours_set[j] == rank) nself++;
1147: if (nself > 1) graph->multi_element = PETSC_TRUE;
1148: }
1149: } else {
1150: PetscCall(PetscCalloc1(nodecount[i], &graph->nodes[i].neighbours_set));
1151: }
1152: }
1153: PetscCall(ISLocalToGlobalMappingRestoreNodeInfo(graph->l2gmap, NULL, &nodecount, &nodeneighs));
1154: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &graph->multi_element, 1, MPI_C_BOOL, MPI_LOR, comm));
1156: /* compute local groups */
1157: if (graph->multi_element) {
1158: const PetscInt *idxs, *indegree;
1159: IS is, lis;
1160: PetscLayout layout;
1161: PetscSF sf, multisf;
1162: PetscInt n, nmulti, c, *multi_root_subs, *start;
1164: PetscCheck(!nvtxs || graph->local_subs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing local subdomain information");
1166: PetscCall(ISLocalToGlobalMappingGetIndices(graph->l2gmap, &idxs));
1167: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nvtxs, idxs, PETSC_USE_POINTER, &is));
1168: PetscCall(ISRenumber(is, NULL, &n, &lis));
1169: PetscCall(ISDestroy(&is));
1171: PetscCall(ISLocalToGlobalMappingRestoreIndices(graph->l2gmap, &idxs));
1172: PetscCall(ISGetIndices(lis, &idxs));
1173: PetscCall(PetscLayoutCreate(PETSC_COMM_SELF, &layout));
1174: PetscCall(PetscLayoutSetSize(layout, n));
1175: PetscCall(PetscSFCreate(PETSC_COMM_SELF, &sf));
1176: PetscCall(PetscSFSetGraphLayout(sf, layout, nvtxs, NULL, PETSC_OWN_POINTER, idxs));
1177: PetscCall(PetscLayoutDestroy(&layout));
1178: PetscCall(PetscSFGetMultiSF(sf, &multisf));
1179: PetscCall(PetscSFComputeDegreeBegin(sf, &indegree));
1180: PetscCall(PetscSFComputeDegreeEnd(sf, &indegree));
1181: PetscCall(PetscSFGetGraph(multisf, &nmulti, NULL, NULL, NULL));
1182: PetscCall(PetscMalloc2(nmulti, &multi_root_subs, n + 1, &start));
1183: start[0] = 0;
1184: for (i = 0; i < n; i++) start[i + 1] = start[i] + indegree[i];
1185: PetscCall(PetscSFGatherBegin(sf, MPIU_INT, graph->local_subs, multi_root_subs));
1186: PetscCall(PetscSFGatherEnd(sf, MPIU_INT, graph->local_subs, multi_root_subs));
1187: for (i = 0; i < nvtxs; i++) {
1188: PetscInt gid = idxs[i];
1190: graph->nodes[i].local_sub = graph->local_subs[i];
1191: for (j = 0, c = 0; j < graph->nodes[i].count; j++) {
1192: if (graph->nodes[i].neighbours_set[j] == rank) c++;
1193: }
1194: PetscCheck(c == indegree[idxs[i]], PETSC_COMM_SELF, PETSC_ERR_PLIB, "%" PetscInt_FMT " != %" PetscInt_FMT, c, indegree[idxs[i]]);
1195: PetscCall(PetscMalloc1(c, &graph->nodes[i].local_groups));
1196: for (j = 0; j < c; j++) graph->nodes[i].local_groups[j] = multi_root_subs[start[gid] + j];
1197: PetscCall(PetscSortInt(c, graph->nodes[i].local_groups));
1198: graph->nodes[i].local_groups_count = c;
1199: }
1200: PetscCall(PetscFree2(multi_root_subs, start));
1201: PetscCall(ISRestoreIndices(lis, &idxs));
1202: PetscCall(ISDestroy(&lis));
1203: PetscCall(PetscSFDestroy(&sf));
1204: }
1206: /*
1207: Get info for dofs splitting
1208: User can specify just a subset; an additional field is considered as a complementary field
1209: */
1210: for (i = 0, k = 0; i < n_ISForDofs; i++) {
1211: PetscInt bs;
1213: PetscCall(ISGetBlockSize(ISForDofs[i], &bs));
1214: k += bs;
1215: }
1216: for (i = 0; i < nvtxs; i++) graph->nodes[i].which_dof = k; /* by default a dof belongs to the complement set */
1217: for (i = 0, k = 0; i < n_ISForDofs; i++) {
1218: PetscInt bs;
1220: PetscCall(ISGetLocalSize(ISForDofs[i], &is_size));
1221: PetscCall(ISGetBlockSize(ISForDofs[i], &bs));
1222: PetscCall(ISGetIndices(ISForDofs[i], &is_indices));
1223: for (j = 0; j < is_size / bs; j++) {
1224: for (PetscInt b = 0; b < bs; b++) {
1225: PetscInt jj = bs * j + b;
1227: if (is_indices[jj] > -1 && is_indices[jj] < nvtxs) { /* out of bounds indices (if any) are skipped */
1228: graph->nodes[is_indices[jj]].which_dof = k + b;
1229: }
1230: }
1231: }
1232: PetscCall(ISRestoreIndices(ISForDofs[i], &is_indices));
1233: k += bs;
1234: }
1236: /* Take into account Neumann nodes */
1237: if (neumann_is) {
1238: PetscCall(ISGetLocalSize(neumann_is, &is_size));
1239: PetscCall(ISGetIndices(neumann_is, &is_indices));
1240: for (i = 0; i < is_size; i++) {
1241: if (is_indices[i] > -1 && is_indices[i] < nvtxs) { /* out of bounds indices (if any) are skipped */
1242: graph->nodes[is_indices[i]].special_dof = PCBDDCGRAPH_NEUMANN_MARK;
1243: }
1244: }
1245: PetscCall(ISRestoreIndices(neumann_is, &is_indices));
1246: }
1248: /* Take into account Dirichlet nodes (they overwrite any mark previously set) */
1249: if (dirichlet_is) {
1250: PetscCall(ISGetLocalSize(dirichlet_is, &is_size));
1251: PetscCall(ISGetIndices(dirichlet_is, &is_indices));
1252: for (i = 0; i < is_size; i++) {
1253: if (is_indices[i] > -1 && is_indices[i] < nvtxs) { /* out of bounds indices (if any) are skipped */
1254: if (!graph->seq_graph) { /* dirichlet nodes treated as internal */
1255: graph->nodes[is_indices[i]].touched = PETSC_TRUE;
1256: graph->nodes[is_indices[i]].subset = 0;
1257: }
1258: graph->nodes[is_indices[i]].special_dof = PCBDDCGRAPH_DIRICHLET_MARK;
1259: }
1260: }
1261: PetscCall(ISRestoreIndices(dirichlet_is, &is_indices));
1262: }
1264: /* mark special nodes (if any) -> each will become a single dof equivalence class (i.e. point constraint for BDDC) */
1265: if (custom_primal_vertices) {
1266: PetscCall(ISGetLocalSize(custom_primal_vertices, &is_size));
1267: PetscCall(ISGetIndices(custom_primal_vertices, &is_indices));
1268: for (i = 0, j = 0; i < is_size; i++) {
1269: if (is_indices[i] > -1 && is_indices[i] < nvtxs && graph->nodes[is_indices[i]].special_dof != PCBDDCGRAPH_DIRICHLET_MARK) { /* out of bounds indices (if any) are skipped */
1270: graph->nodes[is_indices[i]].special_dof = PCBDDCGRAPH_SPECIAL_MARK - j;
1271: j++;
1272: }
1273: }
1274: PetscCall(ISRestoreIndices(custom_primal_vertices, &is_indices));
1275: }
1277: /* mark interior nodes as touched and belonging to partition number 0 */
1278: if (!graph->seq_graph) {
1279: for (i = 0; i < nvtxs; i++) {
1280: const PetscInt icount = graph->nodes[i].count;
1281: if (graph->nodes[i].count < 2) {
1282: graph->nodes[i].touched = PETSC_TRUE;
1283: graph->nodes[i].subset = 0;
1284: } else {
1285: if (graph->multi_element) {
1286: graph->nodes[i].shared = PETSC_FALSE;
1287: for (k = 0; k < icount; k++)
1288: if (graph->nodes[i].neighbours_set[k] != rank) {
1289: graph->nodes[i].shared = PETSC_TRUE;
1290: break;
1291: }
1292: } else {
1293: graph->nodes[i].shared = PETSC_TRUE;
1294: }
1295: }
1296: }
1297: } else {
1298: for (i = 0; i < nvtxs; i++) graph->nodes[i].shared = PETSC_TRUE;
1299: }
1301: /* init graph structure and compute default subsets */
1302: nodes_touched = 0;
1303: for (i = 0; i < nvtxs; i++)
1304: if (graph->nodes[i].touched) nodes_touched++;
1306: /* allocated space for queues */
1307: if (graph->seq_graph) {
1308: PetscCall(PetscMalloc2(nvtxs + 1, &graph->cptr, nvtxs, &graph->queue));
1309: } else {
1310: PetscInt nused = nvtxs - nodes_touched;
1311: PetscCall(PetscMalloc2(nused + 1, &graph->cptr, nused, &graph->queue));
1312: }
1314: graph->ncc = 0;
1315: PetscCall(PetscHMapPCBDDCGraphNodeCreate(&subsetmaps));
1316: PetscCall(PetscCalloc1(nvtxs, &subset_sizes));
1317: for (i = 0; i < nvtxs; i++) {
1318: PetscHashIter iter;
1319: PetscBool missing;
1320: PetscInt subset;
1322: if (graph->nodes[i].touched) continue;
1323: graph->nodes[i].touched = PETSC_TRUE;
1324: PetscCall(PetscHMapPCBDDCGraphNodePut(subsetmaps, &graph->nodes[i], &iter, &missing));
1325: if (missing) {
1326: graph->ncc++;
1327: PetscCall(PetscHMapPCBDDCGraphNodeIterSet(subsetmaps, iter, graph->ncc));
1328: subset = graph->ncc;
1329: } else PetscCall(PetscHMapPCBDDCGraphNodeIterGet(subsetmaps, iter, &subset));
1330: subset_sizes[subset - 1] += 1;
1331: graph->nodes[i].subset = subset;
1332: }
1334: graph->cptr[0] = 0;
1335: for (i = 0; i < graph->ncc; i++) graph->cptr[i + 1] = graph->cptr[i] + subset_sizes[i];
1336: for (i = 0; i < graph->ncc; i++) subset_sizes[i] = 0;
1338: for (i = 0; i < nvtxs; i++) {
1339: const PetscInt subset = graph->nodes[i].subset - 1;
1340: if (subset < 0) continue;
1341: PetscCheck(subset_sizes[subset] + graph->cptr[subset] < graph->cptr[subset + 1], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error for subset %" PetscInt_FMT, subset);
1342: graph->queue[subset_sizes[subset] + graph->cptr[subset]] = i;
1343: subset_sizes[subset] += 1;
1344: }
1345: for (i = 0; i < graph->ncc; i++) PetscCheck(subset_sizes[i] + graph->cptr[i] == graph->cptr[i + 1], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error for subset %" PetscInt_FMT, i);
1347: PetscCall(PetscHMapPCBDDCGraphNodeDestroy(&subsetmaps));
1348: PetscCall(PetscFree(subset_sizes));
1350: /* set default number of subsets */
1351: graph->n_subsets = graph->ncc;
1352: PetscCall(PetscMalloc1(graph->n_subsets, &graph->subset_ncc));
1353: for (i = 0; i < graph->n_subsets; i++) graph->subset_ncc[i] = 1;
1355: PetscCall(PetscMalloc1(graph->ncc, &graph->subset_ref_node));
1356: PetscCall(PetscMalloc1(graph->cptr[graph->ncc], &queue_global));
1357: PetscCall(PetscMalloc2(graph->ncc, &graph->subset_size, graph->ncc, &graph->subset_idxs));
1358: if (graph->multi_element) PetscCall(PetscMalloc1(graph->ncc, &graph->gsubset_size));
1359: else graph->gsubset_size = graph->subset_size;
1360: PetscCall(ISLocalToGlobalMappingApply(graph->l2gmap, graph->cptr[graph->ncc], graph->queue, queue_global));
1362: PetscHMapI cnt_unique;
1364: PetscCall(PetscHMapICreate(&cnt_unique));
1365: for (j = 0; j < graph->ncc; j++) {
1366: PetscInt c = 0, ref_node = PETSC_INT_MAX;
1368: for (k = graph->cptr[j]; k < graph->cptr[j + 1]; k++) {
1369: ref_node = PetscMin(ref_node, queue_global[k]);
1370: if (graph->multi_element) {
1371: PetscBool missing;
1372: PetscHashIter iter;
1374: PetscCall(PetscHMapIPut(cnt_unique, queue_global[k], &iter, &missing));
1375: if (missing) c++;
1376: }
1377: }
1378: graph->gsubset_size[j] = c;
1379: graph->subset_size[j] = graph->cptr[j + 1] - graph->cptr[j];
1380: graph->subset_ref_node[j] = ref_node;
1381: if (graph->multi_element) PetscCall(PetscHMapIClear(cnt_unique));
1382: }
1383: PetscCall(PetscHMapIDestroy(&cnt_unique));
1385: /* save information on subsets (needed when analyzing the connected components) */
1386: if (graph->ncc) {
1387: PetscCall(PetscMalloc1(graph->cptr[graph->ncc], &graph->subset_idxs[0]));
1388: PetscCall(PetscArrayzero(graph->subset_idxs[0], graph->cptr[graph->ncc]));
1389: for (j = 1; j < graph->ncc; j++) graph->subset_idxs[j] = graph->subset_idxs[j - 1] + graph->subset_size[j - 1];
1390: PetscCall(PetscArraycpy(graph->subset_idxs[0], graph->queue, graph->cptr[graph->ncc]));
1391: }
1393: /* check consistency and create SF to analyze components on the interface between subdomains */
1394: if (!graph->seq_graph) {
1395: PetscSF msf;
1396: PetscLayout map;
1397: const PetscInt *degree;
1398: PetscInt nr, nmr, *rdata;
1399: PetscBool valid = PETSC_TRUE;
1400: PetscInt subset_N;
1401: IS subset_n;
1402: const PetscInt *idxs;
1404: PetscCall(ISCreateGeneral(comm, graph->n_subsets, graph->subset_ref_node, PETSC_USE_POINTER, &subset));
1405: PetscCall(ISRenumber(subset, NULL, &subset_N, &subset_n));
1406: PetscCall(ISDestroy(&subset));
1408: PetscCall(PetscSFCreate(comm, &graph->interface_ref_sf));
1409: PetscCall(PetscLayoutCreateFromSizes(comm, PETSC_DECIDE, subset_N, 1, &map));
1410: PetscCall(ISGetIndices(subset_n, &idxs));
1411: PetscCall(PetscSFSetGraphLayout(graph->interface_ref_sf, map, graph->n_subsets, NULL, PETSC_OWN_POINTER, idxs));
1412: PetscCall(ISRestoreIndices(subset_n, &idxs));
1413: PetscCall(ISDestroy(&subset_n));
1414: PetscCall(PetscLayoutDestroy(&map));
1416: PetscCall(PetscSFComputeDegreeBegin(graph->interface_ref_sf, °ree));
1417: PetscCall(PetscSFComputeDegreeEnd(graph->interface_ref_sf, °ree));
1418: PetscCall(PetscSFGetMultiSF(graph->interface_ref_sf, &msf));
1419: PetscCall(PetscSFGetGraph(graph->interface_ref_sf, &nr, NULL, NULL, NULL));
1420: PetscCall(PetscSFGetGraph(msf, &nmr, NULL, NULL, NULL));
1421: PetscCall(PetscCalloc1(nmr, &rdata));
1422: PetscCall(PetscSFGatherBegin(graph->interface_ref_sf, MPIU_INT, graph->gsubset_size, rdata));
1423: PetscCall(PetscSFGatherEnd(graph->interface_ref_sf, MPIU_INT, graph->gsubset_size, rdata));
1424: for (PetscInt i = 0, c = 0; i < nr && valid; i++) {
1425: for (PetscInt j = 0; j < degree[i]; j++) {
1426: if (rdata[j + c] != rdata[c]) valid = PETSC_FALSE;
1427: }
1428: c += degree[i];
1429: }
1430: PetscCall(PetscFree(rdata));
1431: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &valid, 1, MPI_C_BOOL, MPI_LAND, comm));
1432: PetscCheck(valid, comm, PETSC_ERR_PLIB, "Initial local subsets are not consistent");
1434: /* Now create SF with each root extended to gsubset_size roots */
1435: PetscInt mss = 0;
1436: const PetscSFNode *subs_remote;
1438: PetscCall(PetscSFGetGraph(graph->interface_ref_sf, NULL, NULL, NULL, &subs_remote));
1439: for (PetscInt i = 0; i < graph->n_subsets; i++) mss = PetscMax(graph->subset_size[i], mss);
1441: PetscInt nri, nli, *start_rsize, *cum_rsize;
1442: PetscCall(PetscCalloc1(graph->n_subsets + 1, &start_rsize));
1443: PetscCall(PetscCalloc1(nr, &graph->interface_ref_rsize));
1444: PetscCall(PetscMalloc1(nr + 1, &cum_rsize));
1445: PetscCall(PetscSFReduceBegin(graph->interface_ref_sf, MPIU_INT, graph->gsubset_size, graph->interface_ref_rsize, MPI_REPLACE));
1446: PetscCall(PetscSFReduceEnd(graph->interface_ref_sf, MPIU_INT, graph->gsubset_size, graph->interface_ref_rsize, MPI_REPLACE));
1448: nri = 0;
1449: cum_rsize[0] = 0;
1450: for (PetscInt i = 0; i < nr; i++) {
1451: nri += graph->interface_ref_rsize[i];
1452: cum_rsize[i + 1] = cum_rsize[i] + graph->interface_ref_rsize[i];
1453: }
1454: nli = graph->cptr[graph->ncc];
1455: PetscCall(PetscSFBcastBegin(graph->interface_ref_sf, MPIU_INT, cum_rsize, start_rsize, MPI_REPLACE));
1456: PetscCall(PetscSFBcastEnd(graph->interface_ref_sf, MPIU_INT, cum_rsize, start_rsize, MPI_REPLACE));
1457: PetscCall(PetscFree(cum_rsize));
1459: PetscInt *ilocal, *queue_global_uniq;
1460: PetscSFNode *iremote;
1461: PetscBool *touched;
1463: PetscCall(PetscSFCreate(comm, &graph->interface_subset_sf));
1464: PetscCall(PetscMalloc1(nli, &ilocal));
1465: PetscCall(PetscMalloc1(nli, &iremote));
1466: PetscCall(PetscMalloc2(mss, &queue_global_uniq, mss, &touched));
1467: for (PetscInt i = 0, nli = 0; i < graph->n_subsets; i++) {
1468: const PetscMPIInt rr = (PetscMPIInt)subs_remote[i].rank;
1469: const PetscInt start = start_rsize[i];
1470: const PetscInt subset_size = graph->subset_size[i];
1471: const PetscInt gsubset_size = graph->gsubset_size[i];
1472: const PetscInt *subset_idxs = graph->subset_idxs[i];
1473: const PetscInt *lsub_queue_global = queue_global + graph->cptr[i];
1475: k = subset_size;
1476: PetscCall(PetscArrayzero(touched, subset_size));
1477: PetscCall(PetscArraycpy(queue_global_uniq, lsub_queue_global, subset_size));
1478: PetscCall(PetscSortRemoveDupsInt(&k, queue_global_uniq));
1479: PetscCheck(k == gsubset_size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid local subset %" PetscInt_FMT " size %" PetscInt_FMT " != %" PetscInt_FMT, i, k, gsubset_size);
1481: PetscInt t = 0, j = 0;
1482: while (t < subset_size) {
1483: while (j < subset_size && touched[j]) j++;
1484: PetscCheck(j < subset_size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected %" PetscInt_FMT " >= %" PetscInt_FMT, j, subset_size);
1485: const PetscInt ls = graph->nodes[subset_idxs[j]].local_sub;
1487: for (k = j; k < subset_size; k++) {
1488: if (graph->nodes[subset_idxs[k]].local_sub == ls) {
1489: PetscInt ig;
1491: PetscCall(PetscFindInt(lsub_queue_global[k], gsubset_size, queue_global_uniq, &ig));
1492: ilocal[nli] = subset_idxs[k];
1493: iremote[nli].rank = rr;
1494: iremote[nli].index = start + ig;
1495: touched[k] = PETSC_TRUE;
1496: nli++;
1497: t++;
1498: }
1499: }
1500: }
1501: }
1502: PetscCheck(nli == graph->cptr[graph->ncc], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid ilocal size %" PetscInt_FMT " != %" PetscInt_FMT, nli, graph->cptr[graph->ncc]);
1503: PetscCall(PetscSFSetGraph(graph->interface_subset_sf, nri, nli, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
1504: PetscCall(PetscFree(start_rsize));
1505: PetscCall(PetscFree2(queue_global_uniq, touched));
1506: }
1507: PetscCall(PetscFree(queue_global));
1509: /* free workspace */
1510: graph->setupcalled = PETSC_TRUE;
1511: PetscFunctionReturn(PETSC_SUCCESS);
1512: }
1514: PetscErrorCode PCBDDCGraphResetCoords(PCBDDCGraph graph)
1515: {
1516: PetscFunctionBegin;
1517: if (!graph) PetscFunctionReturn(PETSC_SUCCESS);
1518: PetscCall(PetscFree(graph->coords));
1519: graph->cdim = 0;
1520: graph->cnloc = 0;
1521: graph->cloc = PETSC_FALSE;
1522: PetscFunctionReturn(PETSC_SUCCESS);
1523: }
1525: PetscErrorCode PCBDDCGraphResetCSR(PCBDDCGraph graph)
1526: {
1527: PetscFunctionBegin;
1528: if (!graph) PetscFunctionReturn(PETSC_SUCCESS);
1529: if (graph->freecsr) {
1530: PetscCall(PetscFree(graph->xadj));
1531: PetscCall(PetscFree(graph->adjncy));
1532: } else {
1533: graph->xadj = NULL;
1534: graph->adjncy = NULL;
1535: }
1536: graph->freecsr = PETSC_FALSE;
1537: graph->nvtxs_csr = 0;
1538: PetscFunctionReturn(PETSC_SUCCESS);
1539: }
1541: PetscErrorCode PCBDDCGraphReset(PCBDDCGraph graph)
1542: {
1543: PetscFunctionBegin;
1544: if (!graph) PetscFunctionReturn(PETSC_SUCCESS);
1545: PetscCall(ISLocalToGlobalMappingDestroy(&graph->l2gmap));
1546: PetscCall(PetscFree(graph->subset_ncc));
1547: PetscCall(PetscFree(graph->subset_ref_node));
1548: for (PetscInt i = 0; i < graph->nvtxs; i++) {
1549: PetscCall(PetscFree(graph->nodes[i].neighbours_set));
1550: PetscCall(PetscFree(graph->nodes[i].local_groups));
1551: }
1552: PetscCall(PetscFree(graph->nodes));
1553: PetscCall(PetscFree2(graph->cptr, graph->queue));
1554: if (graph->subset_idxs) PetscCall(PetscFree(graph->subset_idxs[0]));
1555: PetscCall(PetscFree2(graph->subset_size, graph->subset_idxs));
1556: if (graph->multi_element) PetscCall(PetscFree(graph->gsubset_size));
1557: PetscCall(PetscFree(graph->interface_ref_rsize));
1558: PetscCall(PetscSFDestroy(&graph->interface_subset_sf));
1559: PetscCall(PetscSFDestroy(&graph->interface_ref_sf));
1560: PetscCall(ISDestroy(&graph->dirdofs));
1561: PetscCall(ISDestroy(&graph->dirdofsB));
1562: if (graph->n_local_subs) PetscCall(PetscFree(graph->local_subs));
1563: graph->multi_element = PETSC_FALSE;
1564: graph->has_dirichlet = PETSC_FALSE;
1565: graph->twodimset = PETSC_FALSE;
1566: graph->twodim = PETSC_FALSE;
1567: graph->nvtxs = 0;
1568: graph->nvtxs_global = 0;
1569: graph->n_subsets = 0;
1570: graph->custom_minimal_size = 1;
1571: graph->n_local_subs = 0;
1572: graph->maxcount = PETSC_INT_MAX;
1573: graph->seq_graph = PETSC_FALSE;
1574: graph->setupcalled = PETSC_FALSE;
1575: PetscFunctionReturn(PETSC_SUCCESS);
1576: }
1578: PetscErrorCode PCBDDCGraphInit(PCBDDCGraph graph, ISLocalToGlobalMapping l2gmap, PetscInt N, PetscInt maxcount)
1579: {
1580: PetscInt n;
1582: PetscFunctionBegin;
1583: PetscAssertPointer(graph, 1);
1587: /* raise an error if already allocated */
1588: PetscCheck(!graph->nvtxs_global, PetscObjectComm((PetscObject)l2gmap), PETSC_ERR_PLIB, "BDDCGraph already initialized");
1589: /* set number of vertices */
1590: PetscCall(PetscObjectReference((PetscObject)l2gmap));
1591: graph->l2gmap = l2gmap;
1592: PetscCall(ISLocalToGlobalMappingGetSize(l2gmap, &n));
1593: graph->nvtxs = n;
1594: graph->nvtxs_global = N;
1595: /* allocate used space */
1596: PetscCall(PetscCalloc1(graph->nvtxs, &graph->nodes));
1597: /* use -1 as a default value for which_dof array */
1598: for (n = 0; n < graph->nvtxs; n++) graph->nodes[n].which_dof = -1;
1600: /* zeroes workspace for values of ncc */
1601: graph->subset_ncc = NULL;
1602: graph->subset_ref_node = NULL;
1603: /* maxcount for cc */
1604: graph->maxcount = maxcount;
1605: PetscFunctionReturn(PETSC_SUCCESS);
1606: }
1608: PetscErrorCode PCBDDCGraphDestroy(PCBDDCGraph *graph)
1609: {
1610: PetscFunctionBegin;
1611: PetscCall(PCBDDCGraphResetCSR(*graph));
1612: PetscCall(PCBDDCGraphResetCoords(*graph));
1613: PetscCall(PCBDDCGraphReset(*graph));
1614: PetscCall(PetscFree(*graph));
1615: PetscFunctionReturn(PETSC_SUCCESS);
1616: }
1618: PetscErrorCode PCBDDCGraphCreate(PCBDDCGraph *graph)
1619: {
1620: PCBDDCGraph new_graph;
1622: PetscFunctionBegin;
1623: PetscCall(PetscNew(&new_graph));
1624: new_graph->custom_minimal_size = 1;
1625: *graph = new_graph;
1626: PetscFunctionReturn(PETSC_SUCCESS);
1627: }