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: }