Actual source code: bddcgraph.c
petsc-3.10.5 2019-03-28
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: PetscReal *wdist;
305: PetscInt n_neigh,*neigh,*n_shared,**shared;
306: PetscInt maxc, ns;
308: PetscBTCreate(graph->nvtxs,&cornerp);
309: ISLocalToGlobalMappingGetInfo(graph->l2gmap,&n_neigh,&neigh,&n_shared,&shared);
310: for (ns = 1, maxc = 0; ns < n_neigh; ns++) maxc = PetscMax(maxc,n_shared[ns]);
311: PetscMalloc1(maxc*graph->cdim,&wdist);
313: for (ns = 1; ns < n_neigh; ns++) { /* first proc is self */
314: PetscReal *anchor,mdist;
315: PetscInt j,k,d,cdim = graph->cdim;
316: PetscInt point1,point2,point3;
317: /*
318: PetscBool isface = PETSC_FALSE;
320: for (i=0;i<n_shared[ns];i++) {
321: if (graph->count[shared[ns][i]] == 1) {
322: isface = PETSC_TRUE;
323: break;
324: }
325: }
326: if (!isface) continue;
327: */
328: /* import coordinates on shared interface */
329: for (j=0,k=0;j<n_shared[ns];j++)
330: for (d=0;d<cdim;d++)
331: wdist[k++] = graph->coords[shared[ns][j]*cdim+d];
333: /* the dofs are sorted by global numbering, so each rank start from the same id and will detect the same corners */
334: anchor = wdist;
336: /* find the farthest point from the starting one */
337: mdist = -1.0;
338: for (j=0,point1=0;j<n_shared[ns];j++) {
339: PetscReal dist = 0.0;
341: for (d=0;d<cdim;d++) dist += (wdist[j*cdim+d]-anchor[d])*(wdist[j*cdim+d]-anchor[d]);
342: if (dist > mdist) { mdist = dist; point1 = j; }
343: }
345: /* find the farthest point from point1 */
346: anchor = wdist + point1*cdim;
347: mdist = -1.0;
348: for (j=0,point2=0;j<n_shared[ns];j++) {
349: PetscReal dist = 0.0;
351: for (d=0;d<cdim;d++) dist += (wdist[j*cdim+d]-anchor[d])*(wdist[j*cdim+d]-anchor[d]);
352: if (dist > mdist) { mdist = dist; point2 = j; }
353: }
355: /* find the third point maximizing the triangle area */
356: point3 = point2;
357: if (cdim > 2) {
358: PetscReal a = 0.0;
360: for (d=0;d<cdim;d++) a += (wdist[point1*cdim+d]-wdist[point2*cdim+d])*(wdist[point1*cdim+d]-wdist[point2*cdim+d]);
361: mdist = -1.0;
362: for (j=0,point3=0;j<n_shared[ns];j++) {
363: PetscReal area,b = 0.0, c = 0.0;
365: for (d=0;d<cdim;d++) {
366: b += (wdist[point1*cdim+d]-wdist[j*cdim+d])*(wdist[point1*cdim+d]-wdist[j*cdim+d]);
367: c += (wdist[point2*cdim+d]-wdist[j*cdim+d])*(wdist[point2*cdim+d]-wdist[j*cdim+d]);
368: }
369: area = (a+b+c)*(-a+b+c)*(a-b+c)*(a+b-c); /* Heron's formula without divisions by 2 */
370: if (area > mdist) { mdist = area; point3 = j; }
371: }
372: }
374: /* all dofs having the same coordinates will be primal */
375: for (j=0;j<n_shared[ns];j++) {
376: PetscBool same[3] = {PETSC_TRUE,PETSC_TRUE,PETSC_TRUE};
378: for (d=0;d<cdim;d++) {
379: same[0] = (PetscBool)(same[0] && (PetscAbsReal(wdist[j*cdim + d]-wdist[point1*cdim+d]) < PETSC_SMALL));
380: same[1] = (PetscBool)(same[1] && (PetscAbsReal(wdist[j*cdim + d]-wdist[point2*cdim+d]) < PETSC_SMALL));
381: same[2] = (PetscBool)(same[2] && (PetscAbsReal(wdist[j*cdim + d]-wdist[point3*cdim+d]) < PETSC_SMALL));
382: }
383: if (same[0] || same[1] || same[2]) {
384: PetscBTSet(cornerp,shared[ns][j]);
385: }
386: }
387: }
388: PetscFree(wdist);
389: ISLocalToGlobalMappingRestoreInfo(graph->l2gmap,&n_neigh,&neigh,&n_shared,&shared);
390: }
392: /* check consistency of connected components among neighbouring subdomains -> it adapt them in case it is needed */
393: MPI_Comm_size(interface_comm,&size);
394: adapt_interface_reduced = PETSC_FALSE;
395: if (size > 1) {
396: PetscInt i;
397: PetscBool adapt_interface = cornerp ? PETSC_TRUE : PETSC_FALSE;
398: for (i=0;i<graph->n_subsets && !adapt_interface;i++) {
399: /* We are not sure that on a given subset of the local interface,
400: with two connected components, the latters be the same among sharing subdomains */
401: if (graph->subset_ncc[i] > 1) adapt_interface = PETSC_TRUE;
402: }
403: MPIU_Allreduce(&adapt_interface,&adapt_interface_reduced,1,MPIU_BOOL,MPI_LOR,interface_comm);
404: }
406: if (graph->n_subsets && adapt_interface_reduced) {
407: PetscBT subset_cc_adapt;
408: MPI_Request *send_requests,*recv_requests;
409: PetscInt *send_buffer,*recv_buffer;
410: PetscInt sum_requests,start_of_recv,start_of_send;
411: PetscInt *cum_recv_counts;
412: PetscInt *labels;
413: PetscInt ncc,cum_queue,mss,mns,j,k,s;
414: PetscInt **refine_buffer=NULL,*private_labels = NULL;
415: PetscBool *subset_has_corn,*recv_buffer_bool,*send_buffer_bool;
417: PetscCalloc1(graph->n_subsets,&subset_has_corn);
418: if (cornerp) {
419: for (i=0;i<graph->n_subsets;i++) {
420: for (j=0;j<graph->subset_size[i];j++) {
421: if (PetscBTLookup(cornerp,graph->subset_idxs[i][j])) {
422: subset_has_corn[i] = PETSC_TRUE;
423: break;
424: }
425: }
426: }
427: }
428: PetscMalloc1(graph->nvtxs,&labels);
429: PetscMemzero(labels,graph->nvtxs*sizeof(*labels));
430: for (i=0,k=0;i<graph->ncc;i++) {
431: PetscInt s = 1;
432: for (j=graph->cptr[i];j<graph->cptr[i+1];j++) {
433: if (cornerp && PetscBTLookup(cornerp,graph->queue[j])) {
434: labels[graph->queue[j]] = k+s;
435: s += 1;
436: } else {
437: labels[graph->queue[j]] = k;
438: }
439: }
440: k += s;
441: }
443: /* allocate some space */
444: PetscMalloc1(graph->n_subsets+1,&cum_recv_counts);
445: PetscMemzero(cum_recv_counts,(graph->n_subsets+1)*sizeof(*cum_recv_counts));
447: /* first count how many neighbours per connected component I will receive from */
448: cum_recv_counts[0] = 0;
449: for (i=0;i<graph->n_subsets;i++) cum_recv_counts[i+1] = cum_recv_counts[i]+graph->count[graph->subset_idxs[i][0]];
450: PetscMalloc1(graph->n_subsets,&send_buffer_bool);
451: PetscMalloc1(cum_recv_counts[graph->n_subsets],&recv_buffer_bool);
452: PetscMalloc2(cum_recv_counts[graph->n_subsets],&send_requests,cum_recv_counts[graph->n_subsets],&recv_requests);
453: for (i=0;i<cum_recv_counts[graph->n_subsets];i++) {
454: send_requests[i] = MPI_REQUEST_NULL;
455: recv_requests[i] = MPI_REQUEST_NULL;
456: }
458: /* exchange with my neighbours the number of my connected components on the subset of interface */
459: sum_requests = 0;
460: for (i=0;i<graph->n_subsets;i++) {
461: send_buffer_bool[i] = (PetscBool)(graph->subset_ncc[i] > 1 || subset_has_corn[i]);
462: }
463: for (i=0;i<graph->n_subsets;i++) {
464: PetscMPIInt neigh,tag;
465: PetscInt count,*neighs;
467: count = graph->count[graph->subset_idxs[i][0]];
468: neighs = graph->neighbours_set[graph->subset_idxs[i][0]];
469: PetscMPIIntCast(2*graph->subset_ref_node[i],&tag);
470: for (k=0;k<count;k++) {
472: PetscMPIIntCast(neighs[k],&neigh);
473: MPI_Isend(send_buffer_bool + i, 1,MPIU_BOOL,neigh,tag,interface_comm,&send_requests[sum_requests]);
474: MPI_Irecv(recv_buffer_bool + sum_requests,1,MPIU_BOOL,neigh,tag,interface_comm,&recv_requests[sum_requests]);
475: sum_requests++;
476: }
477: }
478: MPI_Waitall(sum_requests,recv_requests,MPI_STATUSES_IGNORE);
479: MPI_Waitall(sum_requests,send_requests,MPI_STATUSES_IGNORE);
481: /* determine the subsets I have to adapt (those having more than 1 cc) */
482: PetscBTCreate(graph->n_subsets,&subset_cc_adapt);
483: PetscBTMemzero(graph->n_subsets,subset_cc_adapt);
484: for (i=0;i<graph->n_subsets;i++) {
485: if (graph->subset_ncc[i] > 1 || subset_has_corn[i]) {
486: PetscBTSet(subset_cc_adapt,i);
487: continue;
488: }
489: for (j=cum_recv_counts[i];j<cum_recv_counts[i+1];j++){
490: if (recv_buffer_bool[j]) {
491: PetscBTSet(subset_cc_adapt,i);
492: break;
493: }
494: }
495: }
496: PetscFree(send_buffer_bool);
497: PetscFree(recv_buffer_bool);
498: PetscFree(subset_has_corn);
500: /* determine send/recv buffers sizes */
501: j = 0;
502: mss = 0;
503: for (i=0;i<graph->n_subsets;i++) {
504: if (PetscBTLookup(subset_cc_adapt,i)) {
505: j += graph->subset_size[i];
506: mss = PetscMax(graph->subset_size[i],mss);
507: }
508: }
509: k = 0;
510: mns = 0;
511: for (i=0;i<graph->n_subsets;i++) {
512: if (PetscBTLookup(subset_cc_adapt,i)) {
513: k += (cum_recv_counts[i+1]-cum_recv_counts[i])*graph->subset_size[i];
514: mns = PetscMax(cum_recv_counts[i+1]-cum_recv_counts[i],mns);
515: }
516: }
517: PetscMalloc2(j,&send_buffer,k,&recv_buffer);
519: /* fill send buffer (order matters: subset_idxs ordered by global ordering) */
520: j = 0;
521: for (i=0;i<graph->n_subsets;i++)
522: if (PetscBTLookup(subset_cc_adapt,i))
523: for (k=0;k<graph->subset_size[i];k++)
524: send_buffer[j++] = labels[graph->subset_idxs[i][k]];
526: /* now exchange the data */
527: start_of_recv = 0;
528: start_of_send = 0;
529: sum_requests = 0;
530: for (i=0;i<graph->n_subsets;i++) {
531: if (PetscBTLookup(subset_cc_adapt,i)) {
532: PetscMPIInt neigh,tag;
533: PetscInt size_of_send = graph->subset_size[i];
535: j = graph->subset_idxs[i][0];
536: PetscMPIIntCast(2*graph->subset_ref_node[i]+1,&tag);
537: for (k=0;k<graph->count[j];k++) {
538: PetscMPIIntCast(graph->neighbours_set[j][k],&neigh);
539: MPI_Isend(&send_buffer[start_of_send],size_of_send,MPIU_INT,neigh,tag,interface_comm,&send_requests[sum_requests]);
540: MPI_Irecv(&recv_buffer[start_of_recv],size_of_send,MPIU_INT,neigh,tag,interface_comm,&recv_requests[sum_requests]);
541: start_of_recv += size_of_send;
542: sum_requests++;
543: }
544: start_of_send += size_of_send;
545: }
546: }
547: MPI_Waitall(sum_requests,recv_requests,MPI_STATUSES_IGNORE);
549: /* refine connected components */
550: start_of_recv = 0;
551: /* allocate some temporary space */
552: if (mss) {
553: PetscMalloc1(mss,&refine_buffer);
554: PetscMalloc2(mss*(mns+1),&refine_buffer[0],mss,&private_labels);
555: }
556: ncc = 0;
557: cum_queue = 0;
558: graph->cptr[0] = 0;
559: for (i=0;i<graph->n_subsets;i++) {
560: if (PetscBTLookup(subset_cc_adapt,i)) {
561: PetscInt subset_counter = 0;
562: PetscInt sharingprocs = cum_recv_counts[i+1]-cum_recv_counts[i]+1; /* count myself */
563: PetscInt buffer_size = graph->subset_size[i];
565: /* compute pointers */
566: for (j=1;j<buffer_size;j++) refine_buffer[j] = refine_buffer[j-1] + sharingprocs;
567: /* analyze contributions from subdomains that share the i-th subset
568: The structure of refine_buffer is suitable to find intersections of ccs among sharingprocs.
569: 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)
570: sharing procs connected components:
571: neigh 0: [0 1 4], [2 3], labels [4,7] (2 connected components)
572: neigh 1: [0 1], [2 3 4], labels [3 2] (2 connected components)
573: neigh 2: [0 4], [1], [2 3], labels [1 5 6] (3 connected components)
574: refine_buffer will be filled as:
575: [ 4, 3, 1;
576: 4, 2, 1;
577: 7, 2, 6;
578: 4, 3, 5;
579: 7, 2, 6; ];
580: The connected components in local ordering are [0], [1], [2 3], [4] */
581: /* fill temp_buffer */
582: for (k=0;k<buffer_size;k++) refine_buffer[k][0] = labels[graph->subset_idxs[i][k]];
583: for (j=0;j<sharingprocs-1;j++) {
584: for (k=0;k<buffer_size;k++) refine_buffer[k][j+1] = recv_buffer[start_of_recv+k];
585: start_of_recv += buffer_size;
586: }
587: PetscMemzero(private_labels,buffer_size*sizeof(PetscInt));
588: for (j=0;j<buffer_size;j++) {
589: if (!private_labels[j]) { /* found a new cc */
590: PetscBool same_set;
592: graph->cptr[ncc] = cum_queue;
593: ncc++;
594: subset_counter++;
595: private_labels[j] = subset_counter;
596: graph->queue[cum_queue++] = graph->subset_idxs[i][j];
597: for (k=j+1;k<buffer_size;k++) { /* check for other nodes in new cc */
598: same_set = PETSC_TRUE;
599: for (s=0;s<sharingprocs;s++) {
600: if (refine_buffer[j][s] != refine_buffer[k][s]) {
601: same_set = PETSC_FALSE;
602: break;
603: }
604: }
605: if (same_set) {
606: private_labels[k] = subset_counter;
607: graph->queue[cum_queue++] = graph->subset_idxs[i][k];
608: }
609: }
610: }
611: }
612: graph->cptr[ncc] = cum_queue;
613: graph->subset_ncc[i] = subset_counter;
614: graph->queue_sorted = PETSC_FALSE;
615: } else { /* this subset does not need to be adapted */
616: PetscMemcpy(graph->queue+cum_queue,graph->subset_idxs[i],graph->subset_size[i]*sizeof(PetscInt));
617: ncc++;
618: cum_queue += graph->subset_size[i];
619: graph->cptr[ncc] = cum_queue;
620: }
621: }
622: graph->cptr[ncc] = cum_queue;
623: graph->ncc = ncc;
624: if (mss) {
625: PetscFree2(refine_buffer[0],private_labels);
626: PetscFree(refine_buffer);
627: }
628: PetscFree(labels);
629: MPI_Waitall(sum_requests,send_requests,MPI_STATUSES_IGNORE);
630: PetscFree2(send_requests,recv_requests);
631: PetscFree2(send_buffer,recv_buffer);
632: PetscFree(cum_recv_counts);
633: PetscBTDestroy(&subset_cc_adapt);
634: }
635: PetscBTDestroy(&cornerp);
637: /* Determine if we are in 2D or 3D */
638: if (!graph->twodimset) {
639: PetscBool twodim = PETSC_TRUE;
640: for (i=0;i<graph->ncc;i++) {
641: PetscInt repdof = graph->queue[graph->cptr[i]];
642: PetscInt ccsize = graph->cptr[i+1]-graph->cptr[i];
643: if (graph->count[repdof] > 1 && ccsize > graph->custom_minimal_size) {
644: twodim = PETSC_FALSE;
645: break;
646: }
647: }
648: MPIU_Allreduce(&twodim,&graph->twodim,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)graph->l2gmap));
649: graph->twodimset = PETSC_TRUE;
650: }
651: return(0);
652: }
655: PETSC_STATIC_INLINE PetscErrorCode PCBDDCGraphComputeCC_Private(PCBDDCGraph graph,PetscInt pid,PetscInt* queue_tip,PetscInt n_prev,PetscInt* n_added)
656: {
657: PetscInt i,j,n;
658: PetscInt *xadj = graph->xadj,*adjncy = graph->adjncy;
659: PetscBT touched = graph->touched;
660: PetscBool havecsr = (PetscBool)(!!xadj);
661: PetscBool havesubs = (PetscBool)(!!graph->n_local_subs);
665: n = 0;
666: if (havecsr && !havesubs) {
667: for (i=-n_prev;i<0;i++) {
668: PetscInt start_dof = queue_tip[i];
669: /* we assume that if a dof has a size 1 adjacency list and the corresponding entry is negative, it is connected to all dofs */
670: if (xadj[start_dof+1]-xadj[start_dof] == 1 && adjncy[xadj[start_dof]] < 0) {
671: for (j=0;j<graph->subset_size[pid-1];j++) { /* pid \in [1,graph->n_subsets] */
672: PetscInt dof = graph->subset_idxs[pid-1][j];
673: if (!PetscBTLookup(touched,dof) && graph->subset[dof] == pid) {
674: PetscBTSet(touched,dof);
675: queue_tip[n] = dof;
676: n++;
677: }
678: }
679: } else {
680: for (j=xadj[start_dof];j<xadj[start_dof+1];j++) {
681: PetscInt dof = adjncy[j];
682: if (!PetscBTLookup(touched,dof) && graph->subset[dof] == pid) {
683: PetscBTSet(touched,dof);
684: queue_tip[n] = dof;
685: n++;
686: }
687: }
688: }
689: }
690: } else if (havecsr && havesubs) {
691: PetscInt sid = graph->local_subs[queue_tip[-n_prev]];
692: for (i=-n_prev;i<0;i++) {
693: PetscInt start_dof = queue_tip[i];
694: /* 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 */
695: if (xadj[start_dof+1]-xadj[start_dof] == 1 && adjncy[xadj[start_dof]] < 0) {
696: for (j=0;j<graph->subset_size[pid-1];j++) { /* pid \in [1,graph->n_subsets] */
697: PetscInt dof = graph->subset_idxs[pid-1][j];
698: if (!PetscBTLookup(touched,dof) && graph->subset[dof] == pid && graph->local_subs[dof] == sid) {
699: PetscBTSet(touched,dof);
700: queue_tip[n] = dof;
701: n++;
702: }
703: }
704: } else {
705: for (j=xadj[start_dof];j<xadj[start_dof+1];j++) {
706: PetscInt dof = adjncy[j];
707: if (!PetscBTLookup(touched,dof) && graph->subset[dof] == pid && graph->local_subs[dof] == sid) {
708: PetscBTSet(touched,dof);
709: queue_tip[n] = dof;
710: n++;
711: }
712: }
713: }
714: }
715: } else if (havesubs) { /* sub info only */
716: PetscInt sid = graph->local_subs[queue_tip[-n_prev]];
717: for (j=0;j<graph->subset_size[pid-1];j++) { /* pid \in [1,graph->n_subsets] */
718: PetscInt dof = graph->subset_idxs[pid-1][j];
719: if (!PetscBTLookup(touched,dof) && graph->subset[dof] == pid && graph->local_subs[dof] == sid) {
720: PetscBTSet(touched,dof);
721: queue_tip[n] = dof;
722: n++;
723: }
724: }
725: } else {
726: for (j=0;j<graph->subset_size[pid-1];j++) { /* pid \in [1,graph->n_subsets] */
727: PetscInt dof = graph->subset_idxs[pid-1][j];
728: if (!PetscBTLookup(touched,dof) && graph->subset[dof] == pid) {
729: PetscBTSet(touched,dof);
730: queue_tip[n] = dof;
731: n++;
732: }
733: }
734: }
735: *n_added = n;
736: return(0);
737: }
739: PetscErrorCode PCBDDCGraphComputeConnectedComponentsLocal(PCBDDCGraph graph)
740: {
741: PetscInt ncc,cum_queue,n;
742: PetscMPIInt commsize;
746: if (!graph->setupcalled) SETERRQ(PetscObjectComm((PetscObject)graph->l2gmap),PETSC_ERR_ORDER,"PCBDDCGraphSetUp should be called first");
747: /* quiet return if there isn't any local info */
748: if (!graph->xadj && !graph->n_local_subs) {
749: return(0);
750: }
752: /* reset any previous search of connected components */
753: PetscBTMemzero(graph->nvtxs,graph->touched);
754: MPI_Comm_size(PetscObjectComm((PetscObject)graph->l2gmap),&commsize);
755: if (commsize > graph->commsizelimit) {
756: PetscInt i;
757: for (i=0;i<graph->nvtxs;i++) {
758: if (graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK || !graph->count[i]) {
759: PetscBTSet(graph->touched,i);
760: }
761: }
762: }
764: /* begin search for connected components */
765: cum_queue = 0;
766: ncc = 0;
767: for (n=0;n<graph->n_subsets;n++) {
768: PetscInt pid = n+1; /* partition labeled by 0 is discarded */
769: PetscInt found = 0,prev = 0,first = 0,ncc_pid = 0;
770: while (found != graph->subset_size[n]) {
771: PetscInt added = 0;
772: if (!prev) { /* search for new starting dof */
773: while (PetscBTLookup(graph->touched,graph->subset_idxs[n][first])) first++;
774: PetscBTSet(graph->touched,graph->subset_idxs[n][first]);
775: graph->queue[cum_queue] = graph->subset_idxs[n][first];
776: graph->cptr[ncc] = cum_queue;
777: prev = 1;
778: cum_queue++;
779: found++;
780: ncc_pid++;
781: ncc++;
782: }
783: PCBDDCGraphComputeCC_Private(graph,pid,graph->queue + cum_queue,prev,&added);
784: if (!added) {
785: graph->subset_ncc[n] = ncc_pid;
786: graph->cptr[ncc] = cum_queue;
787: }
788: prev = added;
789: found += added;
790: cum_queue += added;
791: if (added && found == graph->subset_size[n]) {
792: graph->subset_ncc[n] = ncc_pid;
793: graph->cptr[ncc] = cum_queue;
794: }
795: }
796: }
797: graph->ncc = ncc;
798: graph->queue_sorted = PETSC_FALSE;
799: return(0);
800: }
802: PetscErrorCode PCBDDCGraphSetUp(PCBDDCGraph graph, PetscInt custom_minimal_size, IS neumann_is, IS dirichlet_is, PetscInt n_ISForDofs, IS ISForDofs[], IS custom_primal_vertices)
803: {
804: IS subset,subset_n;
805: MPI_Comm comm;
806: const PetscInt *is_indices;
807: PetscInt n_neigh,*neigh,*n_shared,**shared,*queue_global;
808: PetscInt i,j,k,s,total_counts,nodes_touched,is_size;
809: PetscMPIInt commsize;
810: PetscBool same_set,mirrors_found;
815: if (neumann_is) {
818: }
819: graph->has_dirichlet = PETSC_FALSE;
820: if (dirichlet_is) {
823: graph->has_dirichlet = PETSC_TRUE;
824: }
826: for (i=0;i<n_ISForDofs;i++) {
829: }
830: if (custom_primal_vertices) {
833: }
834: PetscObjectGetComm((PetscObject)(graph->l2gmap),&comm);
835: MPI_Comm_size(comm,&commsize);
837: /* custom_minimal_size */
838: graph->custom_minimal_size = custom_minimal_size;
839: /* get info l2gmap and allocate work vectors */
840: ISLocalToGlobalMappingGetInfo(graph->l2gmap,&n_neigh,&neigh,&n_shared,&shared);
841: /* check if we have any local periodic nodes (periodic BCs) */
842: mirrors_found = PETSC_FALSE;
843: if (graph->nvtxs && n_neigh) {
844: for (i=0; i<n_shared[0]; i++) graph->count[shared[0][i]] += 1;
845: for (i=0; i<n_shared[0]; i++) {
846: if (graph->count[shared[0][i]] > 1) {
847: mirrors_found = PETSC_TRUE;
848: break;
849: }
850: }
851: }
852: /* compute local mirrors (if any) */
853: if (mirrors_found) {
854: IS to,from;
855: PetscInt *local_indices,*global_indices;
857: ISCreateStride(PETSC_COMM_SELF,graph->nvtxs,0,1,&to);
858: ISLocalToGlobalMappingApplyIS(graph->l2gmap,to,&from);
859: /* get arrays of local and global indices */
860: PetscMalloc1(graph->nvtxs,&local_indices);
861: ISGetIndices(to,(const PetscInt**)&is_indices);
862: PetscMemcpy(local_indices,is_indices,graph->nvtxs*sizeof(PetscInt));
863: ISRestoreIndices(to,(const PetscInt**)&is_indices);
864: PetscMalloc1(graph->nvtxs,&global_indices);
865: ISGetIndices(from,(const PetscInt**)&is_indices);
866: PetscMemcpy(global_indices,is_indices,graph->nvtxs*sizeof(PetscInt));
867: ISRestoreIndices(from,(const PetscInt**)&is_indices);
868: /* allocate space for mirrors */
869: PetscMalloc2(graph->nvtxs,&graph->mirrors,graph->nvtxs,&graph->mirrors_set);
870: PetscMemzero(graph->mirrors,graph->nvtxs*sizeof(PetscInt));
871: graph->mirrors_set[0] = 0;
873: k=0;
874: for (i=0;i<n_shared[0];i++) {
875: j=shared[0][i];
876: if (graph->count[j] > 1) {
877: graph->mirrors[j]++;
878: k++;
879: }
880: }
881: /* allocate space for set of mirrors */
882: PetscMalloc1(k,&graph->mirrors_set[0]);
883: for (i=1;i<graph->nvtxs;i++)
884: graph->mirrors_set[i]=graph->mirrors_set[i-1]+graph->mirrors[i-1];
886: /* fill arrays */
887: PetscMemzero(graph->mirrors,graph->nvtxs*sizeof(PetscInt));
888: for (j=0;j<n_shared[0];j++) {
889: i=shared[0][j];
890: if (graph->count[i] > 1)
891: graph->mirrors_set[i][graph->mirrors[i]++]=global_indices[i];
892: }
893: PetscSortIntWithArray(graph->nvtxs,global_indices,local_indices);
894: for (i=0;i<graph->nvtxs;i++) {
895: if (graph->mirrors[i] > 0) {
896: PetscFindInt(graph->mirrors_set[i][0],graph->nvtxs,global_indices,&k);
897: j = global_indices[k];
898: while ( k > 0 && global_indices[k-1] == j) k--;
899: for (j=0;j<graph->mirrors[i];j++) {
900: graph->mirrors_set[i][j]=local_indices[k+j];
901: }
902: PetscSortInt(graph->mirrors[i],graph->mirrors_set[i]);
903: }
904: }
905: PetscFree(local_indices);
906: PetscFree(global_indices);
907: ISDestroy(&to);
908: ISDestroy(&from);
909: }
910: PetscMemzero(graph->count,graph->nvtxs*sizeof(*graph->count));
912: /* Count total number of neigh per node */
913: k = 0;
914: for (i=1;i<n_neigh;i++) {
915: k += n_shared[i];
916: for (j=0;j<n_shared[i];j++) {
917: graph->count[shared[i][j]] += 1;
918: }
919: }
920: /* Allocate space for storing the set of neighbours for each node */
921: if (graph->nvtxs) {
922: PetscMalloc1(k,&graph->neighbours_set[0]);
923: }
924: for (i=1;i<graph->nvtxs;i++) { /* dont count myself */
925: graph->neighbours_set[i]=graph->neighbours_set[i-1]+graph->count[i-1];
926: }
927: /* Get information for sharing subdomains */
928: PetscMemzero(graph->count,graph->nvtxs*sizeof(*graph->count));
929: for (i=1;i<n_neigh;i++) { /* dont count myself */
930: s = n_shared[i];
931: for (j=0;j<s;j++) {
932: k = shared[i][j];
933: graph->neighbours_set[k][graph->count[k]] = neigh[i];
934: graph->count[k] += 1;
935: }
936: }
937: /* sort set of sharing subdomains */
938: for (i=0;i<graph->nvtxs;i++) {
939: PetscSortRemoveDupsInt(&graph->count[i],graph->neighbours_set[i]);
940: }
941: /* free memory allocated by ISLocalToGlobalMappingGetInfo */
942: ISLocalToGlobalMappingRestoreInfo(graph->l2gmap,&n_neigh,&neigh,&n_shared,&shared);
944: /*
945: Get info for dofs splitting
946: User can specify just a subset; an additional field is considered as a complementary field
947: */
948: for (i=0;i<graph->nvtxs;i++) graph->which_dof[i] = n_ISForDofs; /* by default a dof belongs to the complement set */
949: for (i=0;i<n_ISForDofs;i++) {
950: ISGetLocalSize(ISForDofs[i],&is_size);
951: ISGetIndices(ISForDofs[i],(const PetscInt**)&is_indices);
952: for (j=0;j<is_size;j++) {
953: if (is_indices[j] > -1 && is_indices[j] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */
954: graph->which_dof[is_indices[j]] = i;
955: }
956: }
957: ISRestoreIndices(ISForDofs[i],(const PetscInt**)&is_indices);
958: }
960: /* Take into account Neumann nodes */
961: if (neumann_is) {
962: ISGetLocalSize(neumann_is,&is_size);
963: ISGetIndices(neumann_is,(const PetscInt**)&is_indices);
964: for (i=0;i<is_size;i++) {
965: if (is_indices[i] > -1 && is_indices[i] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */
966: graph->special_dof[is_indices[i]] = PCBDDCGRAPH_NEUMANN_MARK;
967: }
968: }
969: ISRestoreIndices(neumann_is,(const PetscInt**)&is_indices);
970: }
971: /* Take into account Dirichlet nodes (they overwrite any neumann boundary mark previously set) */
972: if (dirichlet_is) {
973: ISGetLocalSize(dirichlet_is,&is_size);
974: ISGetIndices(dirichlet_is,(const PetscInt**)&is_indices);
975: for (i=0;i<is_size;i++){
976: if (is_indices[i] > -1 && is_indices[i] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */
977: if (commsize > graph->commsizelimit) { /* dirichlet nodes treated as internal */
978: PetscBTSet(graph->touched,is_indices[i]);
979: graph->subset[is_indices[i]] = 0;
980: }
981: graph->special_dof[is_indices[i]] = PCBDDCGRAPH_DIRICHLET_MARK;
982: }
983: }
984: ISRestoreIndices(dirichlet_is,(const PetscInt**)&is_indices);
985: }
986: /* mark local periodic nodes (if any) and adapt CSR graph (if any) */
987: if (graph->mirrors) {
988: for (i=0;i<graph->nvtxs;i++)
989: if (graph->mirrors[i])
990: graph->special_dof[i] = PCBDDCGRAPH_LOCAL_PERIODIC_MARK;
992: if (graph->xadj) {
993: PetscInt *new_xadj,*new_adjncy;
994: /* sort CSR graph */
995: for (i=0;i<graph->nvtxs;i++)
996: PetscSortInt(graph->xadj[i+1]-graph->xadj[i],&graph->adjncy[graph->xadj[i]]);
998: /* adapt local CSR graph in case of local periodicity */
999: k = 0;
1000: for (i=0;i<graph->nvtxs;i++)
1001: for (j=graph->xadj[i];j<graph->xadj[i+1];j++)
1002: k += graph->mirrors[graph->adjncy[j]];
1004: PetscMalloc1(graph->nvtxs+1,&new_xadj);
1005: PetscMalloc1(k+graph->xadj[graph->nvtxs],&new_adjncy);
1006: new_xadj[0] = 0;
1007: for (i=0;i<graph->nvtxs;i++) {
1008: k = graph->xadj[i+1]-graph->xadj[i];
1009: PetscMemcpy(&new_adjncy[new_xadj[i]],&graph->adjncy[graph->xadj[i]],k*sizeof(PetscInt));
1010: new_xadj[i+1] = new_xadj[i]+k;
1011: for (j=graph->xadj[i];j<graph->xadj[i+1];j++) {
1012: k = graph->mirrors[graph->adjncy[j]];
1013: PetscMemcpy(&new_adjncy[new_xadj[i+1]],graph->mirrors_set[graph->adjncy[j]],k*sizeof(PetscInt));
1014: new_xadj[i+1] += k;
1015: }
1016: k = new_xadj[i+1]-new_xadj[i];
1017: PetscSortRemoveDupsInt(&k,&new_adjncy[new_xadj[i]]);
1018: new_xadj[i+1] = new_xadj[i]+k;
1019: }
1020: /* set new CSR into graph */
1021: PetscFree(graph->xadj);
1022: PetscFree(graph->adjncy);
1023: graph->xadj = new_xadj;
1024: graph->adjncy = new_adjncy;
1025: }
1026: }
1028: /* mark special nodes (if any) -> each will become a single node equivalence class */
1029: if (custom_primal_vertices) {
1030: ISGetLocalSize(custom_primal_vertices,&is_size);
1031: ISGetIndices(custom_primal_vertices,(const PetscInt**)&is_indices);
1032: for (i=0,j=0;i<is_size;i++){
1033: 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 */
1034: graph->special_dof[is_indices[i]] = PCBDDCGRAPH_SPECIAL_MARK-j;
1035: j++;
1036: }
1037: }
1038: ISRestoreIndices(custom_primal_vertices,(const PetscInt**)&is_indices);
1039: }
1041: /* mark interior nodes (if commsize > graph->commsizelimit) as touched and belonging to partition number 0 */
1042: if (commsize > graph->commsizelimit) {
1043: for (i=0;i<graph->nvtxs;i++) {
1044: if (!graph->count[i]) {
1045: PetscBTSet(graph->touched,i);
1046: graph->subset[i] = 0;
1047: }
1048: }
1049: }
1051: /* init graph structure and compute default subsets */
1052: nodes_touched = 0;
1053: for (i=0;i<graph->nvtxs;i++) {
1054: if (PetscBTLookup(graph->touched,i)) {
1055: nodes_touched++;
1056: }
1057: }
1058: i = 0;
1059: graph->ncc = 0;
1060: total_counts = 0;
1062: /* allocated space for queues */
1063: if (commsize == graph->commsizelimit) {
1064: PetscMalloc2(graph->nvtxs+1,&graph->cptr,graph->nvtxs,&graph->queue);
1065: } else {
1066: PetscInt nused = graph->nvtxs - nodes_touched;
1067: PetscMalloc2(nused+1,&graph->cptr,nused,&graph->queue);
1068: }
1070: while (nodes_touched<graph->nvtxs) {
1071: /* find first untouched node in local ordering */
1072: while (PetscBTLookup(graph->touched,i)) i++;
1073: PetscBTSet(graph->touched,i);
1074: graph->subset[i] = graph->ncc+1;
1075: graph->cptr[graph->ncc] = total_counts;
1076: graph->queue[total_counts] = i;
1077: total_counts++;
1078: nodes_touched++;
1079: /* now find all other nodes having the same set of sharing subdomains */
1080: for (j=i+1;j<graph->nvtxs;j++) {
1081: /* check for same number of sharing subdomains, dof number and same special mark */
1082: 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]) {
1083: /* check for same set of sharing subdomains */
1084: same_set = PETSC_TRUE;
1085: for (k=0;k<graph->count[j];k++){
1086: if (graph->neighbours_set[i][k] != graph->neighbours_set[j][k]) {
1087: same_set = PETSC_FALSE;
1088: }
1089: }
1090: /* I have found a friend of mine */
1091: if (same_set) {
1092: PetscBTSet(graph->touched,j);
1093: graph->subset[j] = graph->ncc+1;
1094: nodes_touched++;
1095: graph->queue[total_counts] = j;
1096: total_counts++;
1097: }
1098: }
1099: }
1100: graph->ncc++;
1101: }
1102: /* 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 */
1103: graph->n_subsets = graph->ncc;
1104: PetscMalloc1(graph->n_subsets,&graph->subset_ncc);
1105: for (i=0;i<graph->n_subsets;i++) {
1106: graph->subset_ncc[i] = 1;
1107: }
1108: /* final pointer */
1109: graph->cptr[graph->ncc] = total_counts;
1111: /* For consistency reasons (among neighbours), I need to sort (by global ordering) each connected component */
1112: /* Get a reference node (min index in global ordering) for each subset for tagging messages */
1113: PetscMalloc1(graph->ncc,&graph->subset_ref_node);
1114: PetscMalloc1(graph->cptr[graph->ncc],&queue_global);
1115: ISLocalToGlobalMappingApply(graph->l2gmap,graph->cptr[graph->ncc],graph->queue,queue_global);
1116: for (j=0;j<graph->ncc;j++) {
1117: PetscSortIntWithArray(graph->cptr[j+1]-graph->cptr[j],&queue_global[graph->cptr[j]],&graph->queue[graph->cptr[j]]);
1118: graph->subset_ref_node[j] = graph->queue[graph->cptr[j]];
1119: }
1120: PetscFree(queue_global);
1121: graph->queue_sorted = PETSC_TRUE;
1123: /* save information on subsets (needed when analyzing the connected components) */
1124: if (graph->ncc) {
1125: PetscMalloc2(graph->ncc,&graph->subset_size,graph->ncc,&graph->subset_idxs);
1126: PetscMalloc1(graph->cptr[graph->ncc],&graph->subset_idxs[0]);
1127: PetscMemzero(graph->subset_idxs[0],graph->cptr[graph->ncc]*sizeof(PetscInt));
1128: for (j=1;j<graph->ncc;j++) {
1129: graph->subset_size[j-1] = graph->cptr[j] - graph->cptr[j-1];
1130: graph->subset_idxs[j] = graph->subset_idxs[j-1] + graph->subset_size[j-1];
1131: }
1132: graph->subset_size[graph->ncc-1] = graph->cptr[graph->ncc] - graph->cptr[graph->ncc-1];
1133: PetscMemcpy(graph->subset_idxs[0],graph->queue,graph->cptr[graph->ncc]*sizeof(PetscInt));
1134: }
1136: /* renumber reference nodes */
1137: ISCreateGeneral(PetscObjectComm((PetscObject)(graph->l2gmap)),graph->ncc,graph->subset_ref_node,PETSC_COPY_VALUES,&subset_n);
1138: ISLocalToGlobalMappingApplyIS(graph->l2gmap,subset_n,&subset);
1139: ISDestroy(&subset_n);
1140: ISRenumber(subset,NULL,NULL,&subset_n);
1141: ISDestroy(&subset);
1142: ISGetLocalSize(subset_n,&k);
1143: if (k != graph->ncc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Invalid size of new subset! %D != %D",k,graph->ncc);
1144: ISGetIndices(subset_n,&is_indices);
1145: PetscMemcpy(graph->subset_ref_node,is_indices,graph->ncc*sizeof(PetscInt));
1146: ISRestoreIndices(subset_n,&is_indices);
1147: ISDestroy(&subset_n);
1149: /* free workspace */
1150: graph->setupcalled = PETSC_TRUE;
1151: return(0);
1152: }
1154: PetscErrorCode PCBDDCGraphResetCoords(PCBDDCGraph graph)
1155: {
1159: if (!graph) return(0);
1160: PetscFree(graph->coords);
1161: graph->cdim = 0;
1162: graph->cnloc = 0;
1163: graph->cloc = PETSC_FALSE;
1164: return(0);
1165: }
1167: PetscErrorCode PCBDDCGraphResetCSR(PCBDDCGraph graph)
1168: {
1172: if (!graph) return(0);
1173: if (graph->freecsr) {
1174: PetscFree(graph->xadj);
1175: PetscFree(graph->adjncy);
1176: } else {
1177: graph->xadj = NULL;
1178: graph->adjncy = NULL;
1179: }
1180: graph->freecsr = PETSC_FALSE;
1181: graph->nvtxs_csr = 0;
1182: return(0);
1183: }
1185: PetscErrorCode PCBDDCGraphReset(PCBDDCGraph graph)
1186: {
1190: if (!graph) return(0);
1191: ISLocalToGlobalMappingDestroy(&graph->l2gmap);
1192: PetscFree(graph->subset_ncc);
1193: PetscFree(graph->subset_ref_node);
1194: if (graph->nvtxs) {
1195: PetscFree(graph->neighbours_set[0]);
1196: }
1197: PetscBTDestroy(&graph->touched);
1198: PetscFree5(graph->count,
1199: graph->neighbours_set,
1200: graph->subset,
1201: graph->which_dof,
1202: graph->special_dof);
1203: PetscFree2(graph->cptr,graph->queue);
1204: if (graph->mirrors) {
1205: PetscFree(graph->mirrors_set[0]);
1206: }
1207: PetscFree2(graph->mirrors,graph->mirrors_set);
1208: if (graph->subset_idxs) {
1209: PetscFree(graph->subset_idxs[0]);
1210: }
1211: PetscFree2(graph->subset_size,graph->subset_idxs);
1212: ISDestroy(&graph->dirdofs);
1213: ISDestroy(&graph->dirdofsB);
1214: if (graph->n_local_subs) {
1215: PetscFree(graph->local_subs);
1216: }
1217: graph->has_dirichlet = PETSC_FALSE;
1218: graph->twodimset = PETSC_FALSE;
1219: graph->twodim = PETSC_FALSE;
1220: graph->nvtxs = 0;
1221: graph->nvtxs_global = 0;
1222: graph->n_subsets = 0;
1223: graph->custom_minimal_size = 1;
1224: graph->n_local_subs = 0;
1225: graph->maxcount = PETSC_MAX_INT;
1226: graph->setupcalled = PETSC_FALSE;
1227: return(0);
1228: }
1230: PetscErrorCode PCBDDCGraphInit(PCBDDCGraph graph, ISLocalToGlobalMapping l2gmap, PetscInt N, PetscInt maxcount)
1231: {
1232: PetscInt n;
1240: /* raise an error if already allocated */
1241: if (graph->nvtxs_global) SETERRQ(PetscObjectComm((PetscObject)l2gmap),PETSC_ERR_PLIB,"BDDCGraph already initialized");
1242: /* set number of vertices */
1243: PetscObjectReference((PetscObject)l2gmap);
1244: graph->l2gmap = l2gmap;
1245: ISLocalToGlobalMappingGetSize(l2gmap,&n);
1246: graph->nvtxs = n;
1247: graph->nvtxs_global = N;
1248: /* allocate used space */
1249: PetscBTCreate(graph->nvtxs,&graph->touched);
1250: PetscMalloc5(graph->nvtxs,&graph->count,
1251: graph->nvtxs,&graph->neighbours_set,
1252: graph->nvtxs,&graph->subset,
1253: graph->nvtxs,&graph->which_dof,
1254: graph->nvtxs,&graph->special_dof);
1255: /* zeroes memory */
1256: PetscMemzero(graph->count,graph->nvtxs*sizeof(PetscInt));
1257: PetscMemzero(graph->subset,graph->nvtxs*sizeof(PetscInt));
1258: /* use -1 as a default value for which_dof array */
1259: for (n=0;n<graph->nvtxs;n++) graph->which_dof[n] = -1;
1260: PetscMemzero(graph->special_dof,graph->nvtxs*sizeof(PetscInt));
1261: /* zeroes first pointer to neighbour set */
1262: if (graph->nvtxs) {
1263: graph->neighbours_set[0] = 0;
1264: }
1265: /* zeroes workspace for values of ncc */
1266: graph->subset_ncc = 0;
1267: graph->subset_ref_node = 0;
1268: /* maxcount for cc */
1269: graph->maxcount = maxcount;
1270: return(0);
1271: }
1273: PetscErrorCode PCBDDCGraphDestroy(PCBDDCGraph* graph)
1274: {
1278: PCBDDCGraphResetCSR(*graph);
1279: PCBDDCGraphResetCoords(*graph);
1280: PCBDDCGraphReset(*graph);
1281: PetscFree(*graph);
1282: return(0);
1283: }
1285: PetscErrorCode PCBDDCGraphCreate(PCBDDCGraph *graph)
1286: {
1287: PCBDDCGraph new_graph;
1291: PetscNew(&new_graph);
1292: new_graph->custom_minimal_size = 1;
1293: new_graph->commsizelimit = 1;
1294: *graph = new_graph;
1295: return(0);
1296: }