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