Actual source code: bddcgraph.c
petsc-3.6.4 2016-04-12
1: #include <petsc/private/petscimpl.h>
2: #include <../src/ksp/pc/impls/bddc/bddcprivate.h>
3: #include <../src/ksp/pc/impls/bddc/bddcstructs.h>
7: PetscErrorCode PCBDDCGraphGetDirichletDofsB(PCBDDCGraph graph, IS* dirdofs)
8: {
12: if (graph->dirdofsB) {
13: PetscObjectReference((PetscObject)graph->dirdofsB);
14: } else if (graph->has_dirichlet) {
15: PetscInt i,size;
16: PetscInt *dirdofs_idxs;
18: size = 0;
19: for (i=0;i<graph->nvtxs;i++) {
20: if (graph->count[i] && graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK) size++;
21: }
23: PetscMalloc1(size,&dirdofs_idxs);
24: size = 0;
25: for (i=0;i<graph->nvtxs;i++) {
26: if (graph->count[i] && graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK) dirdofs_idxs[size++] = i;
27: }
28: ISCreateGeneral(PETSC_COMM_SELF,size,dirdofs_idxs,PETSC_OWN_POINTER,&graph->dirdofsB);
29: PetscObjectReference((PetscObject)graph->dirdofsB);
30: }
31: *dirdofs = graph->dirdofsB;
32: return(0);
33: }
37: PetscErrorCode PCBDDCGraphGetDirichletDofs(PCBDDCGraph graph, IS* dirdofs)
38: {
42: if (graph->dirdofs) {
43: PetscObjectReference((PetscObject)graph->dirdofs);
44: } else if (graph->has_dirichlet) {
45: PetscInt i,size;
46: PetscInt *dirdofs_idxs;
48: size = 0;
49: for (i=0;i<graph->nvtxs;i++) {
50: if (graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK) size++;
51: }
53: PetscMalloc1(size,&dirdofs_idxs);
54: size = 0;
55: for (i=0;i<graph->nvtxs;i++) {
56: if (graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK) dirdofs_idxs[size++] = i;
57: }
58: ISCreateGeneral(PetscObjectComm((PetscObject)graph->l2gmap),size,dirdofs_idxs,PETSC_OWN_POINTER,&graph->dirdofs);
59: PetscObjectReference((PetscObject)graph->dirdofs);
60: }
61: *dirdofs = graph->dirdofs;
62: return(0);
63: }
67: PetscErrorCode PCBDDCGraphASCIIView(PCBDDCGraph graph, PetscInt verbosity_level, PetscViewer viewer)
68: {
69: PetscInt i,j,tabs;
70: PetscInt* queue_in_global_numbering;
74: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
75: PetscViewerASCIIGetTab(viewer,&tabs);
76: PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");
77: PetscViewerFlush(viewer);
78: PetscViewerASCIISynchronizedPrintf(viewer,"Local BDDC graph for subdomain %04d\n",PetscGlobalRank);
79: PetscViewerASCIISynchronizedPrintf(viewer,"Number of vertices %d\n",graph->nvtxs);
80: if (verbosity_level > 1) {
81: for (i=0;i<graph->nvtxs;i++) {
82: PetscViewerASCIISynchronizedPrintf(viewer,"%d:\n",i);
83: PetscViewerASCIISynchronizedPrintf(viewer," which_dof: %d\n",graph->which_dof[i]);
84: PetscViewerASCIISynchronizedPrintf(viewer," special_dof: %d\n",graph->special_dof[i]);
85: PetscViewerASCIISynchronizedPrintf(viewer," neighbours: %d\n",graph->count[i]);
86: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
87: if (graph->count[i]) {
88: PetscViewerASCIISynchronizedPrintf(viewer," set of neighbours:");
89: for (j=0;j<graph->count[i];j++) {
90: PetscViewerASCIISynchronizedPrintf(viewer," %d",graph->neighbours_set[i][j]);
91: }
92: PetscViewerASCIISynchronizedPrintf(viewer,"\n");
93: }
94: PetscViewerASCIISetTab(viewer,tabs);
95: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
96: if (graph->mirrors) {
97: PetscViewerASCIISynchronizedPrintf(viewer," mirrors: %d\n",graph->mirrors[i]);
98: if (graph->mirrors[i]) {
99: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
100: PetscViewerASCIISynchronizedPrintf(viewer," set of mirrors:");
101: for (j=0;j<graph->mirrors[i];j++) {
102: PetscViewerASCIISynchronizedPrintf(viewer," %d",graph->mirrors_set[i][j]);
103: }
104: PetscViewerASCIISynchronizedPrintf(viewer,"\n");
105: PetscViewerASCIISetTab(viewer,tabs);
106: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
107: }
108: }
109: if (verbosity_level > 2) {
110: if (graph->xadj && graph->adjncy) {
111: PetscViewerASCIISynchronizedPrintf(viewer," local adj list:");
112: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
113: for (j=graph->xadj[i];j<graph->xadj[i+1];j++) {
114: PetscViewerASCIISynchronizedPrintf(viewer," %d",graph->adjncy[j]);
115: }
116: PetscViewerASCIISynchronizedPrintf(viewer,"\n");
117: PetscViewerASCIISetTab(viewer,tabs);
118: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
119: } else {
120: PetscViewerASCIISynchronizedPrintf(viewer," no adj info\n");
121: }
122: }
123: PetscViewerASCIISynchronizedPrintf(viewer," interface subset id: %d\n",graph->subset[i]);
124: if (graph->subset[i] && graph->subset_ncc) {
125: PetscViewerASCIISynchronizedPrintf(viewer," ncc for subset: %d\n",graph->subset_ncc[graph->subset[i]-1]);
126: }
127: }
128: }
129: PetscViewerASCIISynchronizedPrintf(viewer,"Total number of connected components %d\n",graph->ncc);
130: PetscMalloc1(graph->cptr[graph->ncc],&queue_in_global_numbering);
131: ISLocalToGlobalMappingApply(graph->l2gmap,graph->cptr[graph->ncc],graph->queue,queue_in_global_numbering);
132: for (i=0;i<graph->ncc;i++) {
133: PetscInt node_num=graph->queue[graph->cptr[i]];
134: PetscBool printcc = PETSC_FALSE;
135: PetscViewerASCIISynchronizedPrintf(viewer," %d (neighs:",i);
136: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
137: for (j=0;j<graph->count[node_num];j++) {
138: PetscViewerASCIISynchronizedPrintf(viewer," %d",graph->neighbours_set[node_num][j]);
139: }
140: PetscViewerASCIISynchronizedPrintf(viewer,"):");
141: if (verbosity_level > 1 || graph->twodim) {
142: printcc = PETSC_TRUE;
143: } else if (graph->count[node_num] > 1 || (graph->count[node_num] == 1 && graph->special_dof[node_num] == PCBDDCGRAPH_NEUMANN_MARK)) {
144: printcc = PETSC_TRUE;
145: }
146: if (printcc) {
147: for (j=graph->cptr[i];j<graph->cptr[i+1];j++) {
148: PetscViewerASCIISynchronizedPrintf(viewer," %d (%d)",graph->queue[j],queue_in_global_numbering[j]);
149: }
150: }
151: PetscViewerASCIISynchronizedPrintf(viewer,"\n");
152: PetscViewerASCIISetTab(viewer,tabs);
153: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
154: }
155: PetscFree(queue_in_global_numbering);
156: if (graph->custom_minimal_size > 1 && verbosity_level > 1) {
157: PetscViewerASCIISynchronizedPrintf(viewer,"Custom minimal size %d\n",graph->custom_minimal_size);
158: }
159: PetscViewerFlush(viewer);
160: return(0);
161: }
165: PetscErrorCode PCBDDCGraphGetCandidatesIS(PCBDDCGraph graph, PetscInt *n_faces, IS *FacesIS[], PetscInt *n_edges, IS *EdgesIS[], IS *VerticesIS)
166: {
167: IS *ISForFaces,*ISForEdges,ISForVertices;
168: PetscInt i,nfc,nec,nvc,*idx;
172: /* loop on ccs to evalute number of faces, edges and vertices */
173: nfc = 0;
174: nec = 0;
175: nvc = 0;
176: for (i=0;i<graph->ncc;i++) {
177: PetscInt repdof = graph->queue[graph->cptr[i]];
178: if (graph->cptr[i+1]-graph->cptr[i] > graph->custom_minimal_size) {
179: if (graph->count[repdof] == 1 && graph->special_dof[repdof] != PCBDDCGRAPH_NEUMANN_MARK) {
180: nfc++;
181: } else { /* note that nec will be zero in 2d */
182: nec++;
183: }
184: } else {
185: nvc += graph->cptr[i+1]-graph->cptr[i];
186: }
187: }
189: /* check if we are in 2D or 3D */
190: if (graph->twodim) { /* we are in a 2D case -> edges are shared by 2 subregions and faces don't exist */
191: nec = nfc;
192: nfc = 0;
193: }
195: /* allocate IS arrays for faces, edges. Vertices need a single index set. */
196: if (FacesIS) {
197: PetscMalloc1(nfc,&ISForFaces);
198: }
199: if (EdgesIS) {
200: PetscMalloc1(nec,&ISForEdges);
201: }
202: if (VerticesIS) {
203: PetscMalloc1(nvc,&idx);
204: }
206: /* loop on ccs to compute index sets for faces and edges */
207: if (!graph->queue_sorted) {
208: PetscInt *queue_global;
210: PetscMalloc1(graph->cptr[graph->ncc],&queue_global);
211: ISLocalToGlobalMappingApply(graph->l2gmap,graph->cptr[graph->ncc],graph->queue,queue_global);
212: for (i=0;i<graph->ncc;i++) {
213: PetscSortIntWithArray(graph->cptr[i+1]-graph->cptr[i],&queue_global[graph->cptr[i]],&graph->queue[graph->cptr[i]]);
214: }
215: PetscFree(queue_global);
216: graph->queue_sorted = PETSC_TRUE;
217: }
218: nfc = 0;
219: nec = 0;
220: for (i=0;i<graph->ncc;i++) {
221: PetscInt repdof = graph->queue[graph->cptr[i]];
222: if (graph->cptr[i+1]-graph->cptr[i] > graph->custom_minimal_size) {
223: if (graph->count[repdof] == 1 && graph->special_dof[repdof] != PCBDDCGRAPH_NEUMANN_MARK) {
224: if (graph->twodim) {
225: if (EdgesIS) {
226: ISCreateGeneral(PETSC_COMM_SELF,graph->cptr[i+1]-graph->cptr[i],&graph->queue[graph->cptr[i]],PETSC_COPY_VALUES,&ISForEdges[nec]);
227: }
228: nec++;
229: } else {
230: if (FacesIS) {
231: ISCreateGeneral(PETSC_COMM_SELF,graph->cptr[i+1]-graph->cptr[i],&graph->queue[graph->cptr[i]],PETSC_COPY_VALUES,&ISForFaces[nfc]);
232: }
233: nfc++;
234: }
235: } else {
236: if (EdgesIS) {
237: ISCreateGeneral(PETSC_COMM_SELF,graph->cptr[i+1]-graph->cptr[i],&graph->queue[graph->cptr[i]],PETSC_COPY_VALUES,&ISForEdges[nec]);
238: }
239: nec++;
240: }
241: }
242: }
243: /* index set for vertices */
244: if (VerticesIS) {
245: nvc = 0;
246: for (i=0;i<graph->ncc;i++) {
247: if (graph->cptr[i+1]-graph->cptr[i] <= graph->custom_minimal_size) {
248: PetscInt j;
250: for (j=graph->cptr[i];j<graph->cptr[i+1];j++) {
251: idx[nvc]=graph->queue[j];
252: nvc++;
253: }
254: }
255: }
256: /* sort vertex set (by local ordering) */
257: PetscSortInt(nvc,idx);
258: ISCreateGeneral(PETSC_COMM_SELF,nvc,idx,PETSC_OWN_POINTER,&ISForVertices);
259: }
260: /* get back info */
261: if (n_faces) *n_faces = nfc;
262: if (FacesIS) {
263: *FacesIS = ISForFaces;
264: }
265: if (n_edges) *n_edges = nec;
266: if (EdgesIS) {
267: *EdgesIS = ISForEdges;
268: }
269: if (VerticesIS) {
270: *VerticesIS = ISForVertices;
271: }
272: return(0);
273: }
277: PetscErrorCode PCBDDCGraphComputeConnectedComponents(PCBDDCGraph graph)
278: {
279: PetscBool adapt_interface_reduced;
280: MPI_Comm interface_comm;
281: PetscMPIInt size;
282: PetscInt i;
283: PetscBool twodim;
287: /* compute connected components locally */
288: PetscObjectGetComm((PetscObject)(graph->l2gmap),&interface_comm);
289: PCBDDCGraphComputeConnectedComponentsLocal(graph);
290: /* check consistency of connected components among neighbouring subdomains -> it adapt them in case it is needed */
291: MPI_Comm_size(interface_comm,&size);
292: adapt_interface_reduced = PETSC_FALSE;
293: if (size > 1) {
294: PetscInt i;
295: PetscBool adapt_interface = PETSC_FALSE;
296: for (i=0;i<graph->n_subsets;i++) {
297: /* We are not sure that on a given subset of the local interface,
298: with two connected components, the latters be the same among sharing subdomains */
299: if (graph->subset_ncc[i] > 1) {
300: adapt_interface = PETSC_TRUE;
301: break;
302: }
303: }
304: MPI_Allreduce(&adapt_interface,&adapt_interface_reduced,1,MPIU_BOOL,MPI_LOR,interface_comm);
305: }
307: if (graph->n_subsets && adapt_interface_reduced) {
308: /* old csr stuff */
309: PetscInt *aux_new_xadj,*new_xadj,*new_adjncy;
310: PetscInt *old_xadj,*old_adjncy;
311: PetscBT subset_cc_adapt;
312: /* MPI */
313: MPI_Request *send_requests,*recv_requests;
314: PetscInt *send_buffer,*recv_buffer;
315: PetscInt sum_requests,start_of_recv,start_of_send;
316: PetscInt *cum_recv_counts;
317: /* others */
318: PetscInt *labels;
319: PetscInt global_subset_counter;
320: PetscInt j,k,s;
322: /* Retrict adjacency graph using information from previously computed connected components */
323: PetscMalloc1(graph->nvtxs,&aux_new_xadj);
324: for (i=0;i<graph->nvtxs;i++) {
325: aux_new_xadj[i]=1;
326: }
327: for (i=0;i<graph->ncc;i++) {
328: k = graph->cptr[i+1]-graph->cptr[i];
329: for (j=0;j<k;j++) {
330: aux_new_xadj[graph->queue[graph->cptr[i]+j]]=k;
331: }
332: }
333: j = 0;
334: for (i=0;i<graph->nvtxs;i++) {
335: j += aux_new_xadj[i];
336: }
337: PetscMalloc1(graph->nvtxs+1,&new_xadj);
338: PetscMalloc1(j,&new_adjncy);
339: new_xadj[0]=0;
340: for (i=0;i<graph->nvtxs;i++) {
341: new_xadj[i+1]=new_xadj[i]+aux_new_xadj[i];
342: if (aux_new_xadj[i]==1) {
343: new_adjncy[new_xadj[i]]=i;
344: }
345: }
346: PetscFree(aux_new_xadj);
347: PetscMalloc1(graph->nvtxs,&labels);
348: PetscMemzero(labels,graph->nvtxs*sizeof(*labels));
349: for (i=0;i<graph->ncc;i++) {
350: k = graph->cptr[i+1]-graph->cptr[i];
351: for (j=0;j<k;j++) {
352: PetscMemcpy(&new_adjncy[new_xadj[graph->queue[graph->cptr[i]+j]]],&graph->queue[graph->cptr[i]],k*sizeof(PetscInt));
353: labels[graph->queue[graph->cptr[i]+j]] = i;
354: }
355: }
356: /* set temporarly new CSR into graph */
357: old_xadj = graph->xadj;
358: old_adjncy = graph->adjncy;
359: graph->xadj = new_xadj;
360: graph->adjncy = new_adjncy;
361: /* allocate some space */
362: PetscMalloc1(graph->n_subsets+1,&cum_recv_counts);
363: PetscMemzero(cum_recv_counts,(graph->n_subsets+1)*sizeof(*cum_recv_counts));
364: /* first count how many neighbours per connected component I will receive from */
365: cum_recv_counts[0]=0;
366: for (i=0;i<graph->n_subsets;i++) {
367: cum_recv_counts[i+1]=cum_recv_counts[i]+graph->count[graph->subsets[i][0]];
368: }
369: PetscMalloc1(cum_recv_counts[graph->n_subsets],&recv_buffer);
370: PetscMalloc2(cum_recv_counts[graph->n_subsets],&send_requests,cum_recv_counts[graph->n_subsets],&recv_requests);
371: for (i=0;i<cum_recv_counts[graph->n_subsets];i++) {
372: send_requests[i]=MPI_REQUEST_NULL;
373: recv_requests[i]=MPI_REQUEST_NULL;
374: }
375: /* exchange with my neighbours the number of my connected components on the subset of interface */
376: sum_requests = 0;
377: for (i=0;i<graph->n_subsets;i++) {
378: PetscMPIInt neigh,tag;
380: j = graph->subsets[i][0];
381: PetscMPIIntCast(2*graph->subset_ref_node[i],&tag);
382: for (k=0;k<graph->count[j];k++) {
383: PetscMPIIntCast(graph->neighbours_set[j][k],&neigh);
384: MPI_Isend(&graph->subset_ncc[i],1,MPIU_INT,neigh,tag,interface_comm,&send_requests[sum_requests]);
385: MPI_Irecv(&recv_buffer[sum_requests],1,MPIU_INT,neigh,tag,interface_comm,&recv_requests[sum_requests]);
386: sum_requests++;
387: }
388: }
389: MPI_Waitall(sum_requests,recv_requests,MPI_STATUSES_IGNORE);
390: MPI_Waitall(sum_requests,send_requests,MPI_STATUSES_IGNORE);
391: /* determine the subsets I have to adapt */
392: PetscBTCreate(graph->n_subsets,&subset_cc_adapt);
393: PetscBTMemzero(graph->n_subsets,subset_cc_adapt);
394: for (i=0;i<graph->n_subsets;i++) {
395: for (j=cum_recv_counts[i];j<cum_recv_counts[i+1];j++){
396: /* adapt if more than one cc is present */
397: if (graph->subset_ncc[i] > 1 || recv_buffer[j] > 1) {
398: PetscBTSet(subset_cc_adapt,i);
399: break;
400: }
401: }
402: }
403: PetscFree(recv_buffer);
404: /* determine send/recv buffers sizes */
405: j = 0;
406: for (i=0;i<graph->n_subsets;i++) {
407: if (PetscBTLookup(subset_cc_adapt,i)) {
408: j += graph->subsets_size[i];
409: }
410: }
411: k = 0;
412: for (i=0;i<graph->n_subsets;i++) {
413: if (PetscBTLookup(subset_cc_adapt,i)) {
414: k += (cum_recv_counts[i+1]-cum_recv_counts[i])*graph->subsets_size[i];
415: }
416: }
417: PetscMalloc2(j,&send_buffer,k,&recv_buffer);
419: /* fill send buffer */
420: j = 0;
421: for (i=0;i<graph->n_subsets;i++) {
422: if (PetscBTLookup(subset_cc_adapt,i)) {
423: for (k=0;k<graph->subsets_size[i];k++) {
424: send_buffer[j++] = labels[graph->subsets[i][k]];
425: }
426: }
427: }
428: PetscFree(labels);
430: /* now exchange the data */
431: start_of_recv = 0;
432: start_of_send = 0;
433: sum_requests = 0;
434: for (i=0;i<graph->n_subsets;i++) {
435: if (PetscBTLookup(subset_cc_adapt,i)) {
436: PetscMPIInt neigh,tag;
437: PetscInt size_of_send = graph->subsets_size[i];
439: j = graph->subsets[i][0];
440: PetscMPIIntCast(2*graph->subset_ref_node[i]+1,&tag);
441: for (k=0;k<graph->count[j];k++) {
442: PetscMPIIntCast(graph->neighbours_set[j][k],&neigh);
443: MPI_Isend(&send_buffer[start_of_send],size_of_send,MPIU_INT,neigh,tag,interface_comm,&send_requests[sum_requests]);
444: MPI_Irecv(&recv_buffer[start_of_recv],size_of_send,MPIU_INT,neigh,tag,interface_comm,&recv_requests[sum_requests]);
445: start_of_recv += size_of_send;
446: sum_requests++;
447: }
448: start_of_send += size_of_send;
449: }
450: }
451: MPI_Waitall(sum_requests,recv_requests,MPI_STATUSES_IGNORE);
452: MPI_Waitall(sum_requests,send_requests,MPI_STATUSES_IGNORE);
453: PetscFree2(send_requests,recv_requests);
455: start_of_recv = 0;
456: global_subset_counter = 0;
457: for (i=0;i<graph->n_subsets;i++) {
458: if (PetscBTLookup(subset_cc_adapt,i)) {
459: PetscInt **temp_buffer,temp_buffer_size,*private_label;
461: /* allocate some temporary space */
462: temp_buffer_size = graph->subsets_size[i];
463: PetscMalloc1(temp_buffer_size,&temp_buffer);
464: PetscMalloc1(temp_buffer_size*(cum_recv_counts[i+1]-cum_recv_counts[i]),&temp_buffer[0]);
465: PetscMemzero(temp_buffer[0],temp_buffer_size*(cum_recv_counts[i+1]-cum_recv_counts[i])*sizeof(PetscInt));
466: for (j=1;j<temp_buffer_size;j++) {
467: temp_buffer[j] = temp_buffer[j-1]+cum_recv_counts[i+1]-cum_recv_counts[i];
468: }
469: /* analyze contributions from neighbouring subdomains for i-th subset
470: temp buffer structure:
471: supposing part of the interface has dimension 5
472: e.g with global dofs 0,1,2,3,4, locally ordered on the current process as 0,4,3,1,2
473: 3 neighs procs having connected components:
474: neigh 0: [0 1 4], [2 3], labels [4,7] (2 connected components)
475: neigh 1: [0 1], [2 3 4], labels [3 2] (2 connected components)
476: neigh 2: [0 4], [1], [2 3], labels [1 5 6] (3 connected components)
477: tempbuffer (row-oriented) will be filled as:
478: [ 4, 3, 1;
479: 4, 2, 1;
480: 7, 2, 6;
481: 4, 3, 5;
482: 7, 2, 6; ];
483: This way we can simply find intersections of ccs among neighs.
484: The graph->subset array will need to be modified. The output for the example is [0], [1], [2 3], [4];
485: */
486: for (j=0;j<cum_recv_counts[i+1]-cum_recv_counts[i];j++) {
487: for (k=0;k<temp_buffer_size;k++) {
488: temp_buffer[k][j] = recv_buffer[start_of_recv+k];
489: }
490: start_of_recv += temp_buffer_size;
491: }
492: PetscMalloc1(temp_buffer_size,&private_label);
493: PetscMemzero(private_label,temp_buffer_size*sizeof(*private_label));
494: for (j=0;j<temp_buffer_size;j++) {
495: if (!private_label[j]) { /* found a new cc */
496: PetscBool same_set;
498: global_subset_counter++;
499: private_label[j] = global_subset_counter;
500: for (k=j+1;k<temp_buffer_size;k++) { /* check for other nodes in new cc */
501: same_set = PETSC_TRUE;
502: for (s=0;s<cum_recv_counts[i+1]-cum_recv_counts[i];s++) {
503: if (temp_buffer[j][s] != temp_buffer[k][s]) {
504: same_set = PETSC_FALSE;
505: break;
506: }
507: }
508: if (same_set) {
509: private_label[k] = global_subset_counter;
510: }
511: }
512: }
513: }
514: /* insert private labels in graph structure */
515: for (j=0;j<graph->subsets_size[i];j++) {
516: graph->subset[graph->subsets[i][j]] = graph->n_subsets+private_label[j];
517: }
518: PetscFree(private_label);
519: PetscFree(temp_buffer[0]);
520: PetscFree(temp_buffer);
521: }
522: }
523: PetscFree2(send_buffer,recv_buffer);
524: PetscFree(cum_recv_counts);
525: PetscBTDestroy(&subset_cc_adapt);
526: /* We are ready to find for connected components which are consistent among neighbouring subdomains */
527: if (global_subset_counter) {
528: PetscBTMemzero(graph->nvtxs,graph->touched);
529: global_subset_counter = 0;
530: for (i=0;i<graph->nvtxs;i++) {
531: if (graph->subset[i] && !PetscBTLookup(graph->touched,i)) {
532: global_subset_counter++;
533: for (j=i+1;j<graph->nvtxs;j++) {
534: if (!PetscBTLookup(graph->touched,j) && graph->subset[j]==graph->subset[i]) {
535: graph->subset[j] = global_subset_counter;
536: PetscBTSet(graph->touched,j);
537: }
538: }
539: graph->subset[i] = global_subset_counter;
540: PetscBTSet(graph->touched,i);
541: }
542: }
543: /* refine connected components locally */
544: PCBDDCGraphComputeConnectedComponentsLocal(graph);
545: }
546: /* restore original CSR graph of dofs */
547: PetscFree(new_xadj);
548: PetscFree(new_adjncy);
549: graph->xadj = old_xadj;
550: graph->adjncy = old_adjncy;
551: }
553: /* Determine if we are in 2D or 3D */
554: twodim = PETSC_TRUE;
555: for (i=0;i<graph->ncc;i++) {
556: PetscInt repdof = graph->queue[graph->cptr[i]];
557: if (graph->cptr[i+1]-graph->cptr[i] > graph->custom_minimal_size) {
558: if (graph->count[repdof] > 1 || graph->special_dof[repdof] == PCBDDCGRAPH_NEUMANN_MARK) {
559: twodim = PETSC_FALSE;
560: break;
561: }
562: }
563: }
564: MPI_Allreduce(&twodim,&graph->twodim,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)graph->l2gmap));
565: return(0);
566: }
568: /* The following code has been adapted from function IsConnectedSubdomain contained
569: in source file contig.c of METIS library (version 5.0.1)
570: It finds connected components for each subset */
573: PetscErrorCode PCBDDCGraphComputeConnectedComponentsLocal(PCBDDCGraph graph)
574: {
575: PetscInt i,j,k,first,last,nleft,ncc,pid,cum_queue,n,ncc_pid;
576: PetscMPIInt size;
580: /* quiet return if no csr info is available */
581: if (!graph->xadj || !graph->adjncy) {
582: return(0);
583: }
585: /* reset any previous search of connected components */
586: PetscBTMemzero(graph->nvtxs,graph->touched);
587: graph->n_subsets = 0;
588: MPI_Comm_size(PetscObjectComm((PetscObject)graph->l2gmap),&size);
589: if (size == 1) {
590: if (graph->nvtxs) {
591: graph->n_subsets = 1;
592: for (i=0;i<graph->nvtxs;i++) {
593: graph->subset[i] = 1;
594: }
595: }
596: } else {
597: for (i=0;i<graph->nvtxs;i++) {
598: if (graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK || !graph->count[i]) {
599: PetscBTSet(graph->touched,i);
600: graph->subset[i] = 0;
601: }
602: graph->n_subsets = PetscMax(graph->n_subsets,graph->subset[i]);
603: }
604: }
605: PetscFree(graph->subset_ncc);
606: PetscMalloc1(graph->n_subsets,&graph->subset_ncc);
608: /* begin search for connected components */
609: cum_queue = 0;
610: ncc = 0;
611: for (n=0;n<graph->n_subsets;n++) {
612: pid = n+1; /* partition labeled by 0 is discarded */
613: nleft = 0;
614: for (i=0;i<graph->nvtxs;i++) {
615: if (graph->subset[i] == pid) {
616: nleft++;
617: }
618: }
619: for (i=0; i<graph->nvtxs; i++) {
620: if (graph->subset[i] == pid) {
621: break;
622: }
623: }
624: PetscBTSet(graph->touched,i);
625: graph->queue[cum_queue] = i;
626: first = 0;
627: last = 1;
628: graph->cptr[ncc] = cum_queue;
629: ncc_pid = 0;
630: while (first != nleft) {
631: if (first == last) {
632: graph->cptr[++ncc] = first+cum_queue;
633: ncc_pid++;
634: for (i=0; i<graph->nvtxs; i++) { /* TODO-> use a while! */
635: if (graph->subset[i] == pid && !PetscBTLookup(graph->touched,i)) {
636: break;
637: }
638: }
639: graph->queue[cum_queue+last] = i;
640: last++;
641: PetscBTSet(graph->touched,i);
642: }
643: i = graph->queue[cum_queue+first];
644: first++;
645: for (j=graph->xadj[i];j<graph->xadj[i+1];j++) {
646: k = graph->adjncy[j];
647: if (graph->subset[k] == pid && !PetscBTLookup(graph->touched,k)) {
648: graph->queue[cum_queue+last] = k;
649: last++;
650: PetscBTSet(graph->touched,k);
651: }
652: }
653: }
654: graph->cptr[++ncc] = first+cum_queue;
655: ncc_pid++;
656: cum_queue = graph->cptr[ncc];
657: graph->subset_ncc[n] = ncc_pid;
658: }
659: graph->ncc = ncc;
660: graph->queue_sorted = PETSC_FALSE;
661: return(0);
662: }
666: PetscErrorCode PCBDDCGraphSetUp(PCBDDCGraph graph, PetscInt custom_minimal_size, IS neumann_is, IS dirichlet_is, PetscInt n_ISForDofs, IS ISForDofs[], IS custom_primal_vertices)
667: {
668: VecScatter scatter_ctx;
669: Vec local_vec,local_vec2,global_vec;
670: IS to,from,subset,subset_n;
671: MPI_Comm comm;
672: PetscScalar *array,*array2;
673: const PetscInt *is_indices;
674: PetscInt n_neigh,*neigh,*n_shared,**shared,*queue_global;
675: PetscInt i,j,k,s,total_counts,nodes_touched,is_size;
676: PetscMPIInt size;
677: PetscBool same_set,mirrors_found;
681: graph->has_dirichlet = PETSC_FALSE;
682: if (dirichlet_is) {
684: graph->has_dirichlet = PETSC_TRUE;
685: }
686: PetscObjectGetComm((PetscObject)(graph->l2gmap),&comm);
687: /* custom_minimal_size */
688: graph->custom_minimal_size = PetscMax(graph->custom_minimal_size,custom_minimal_size);
689: /* get info l2gmap and allocate work vectors */
690: ISLocalToGlobalMappingGetInfo(graph->l2gmap,&n_neigh,&neigh,&n_shared,&shared);
691: /* check if we have any local periodic nodes (periodic BCs) */
692: mirrors_found = PETSC_FALSE;
693: if (graph->nvtxs && n_neigh) {
694: for (i=0; i<n_shared[0]; i++) graph->count[shared[0][i]] += 1;
695: for (i=0; i<n_shared[0]; i++) {
696: if (graph->count[shared[0][i]] > 1) {
697: mirrors_found = PETSC_TRUE;
698: break;
699: }
700: }
701: }
702: /* create some workspace objects */
703: local_vec = NULL;
704: local_vec2 = NULL;
705: global_vec = NULL;
706: to = NULL;
707: from = NULL;
708: scatter_ctx = NULL;
709: if (n_ISForDofs || dirichlet_is || neumann_is || custom_primal_vertices) {
710: VecCreate(PETSC_COMM_SELF,&local_vec);
711: VecSetSizes(local_vec,PETSC_DECIDE,graph->nvtxs);
712: VecSetType(local_vec,VECSTANDARD);
713: VecDuplicate(local_vec,&local_vec2);
714: VecCreate(comm,&global_vec);
715: VecSetSizes(global_vec,PETSC_DECIDE,graph->nvtxs_global);
716: VecSetType(global_vec,VECSTANDARD);
717: ISCreateStride(PETSC_COMM_SELF,graph->nvtxs,0,1,&to);
718: ISLocalToGlobalMappingApplyIS(graph->l2gmap,to,&from);
719: VecScatterCreate(global_vec,from,local_vec,to,&scatter_ctx);
720: } else if (mirrors_found) {
721: ISCreateStride(PETSC_COMM_SELF,graph->nvtxs,0,1,&to);
722: ISLocalToGlobalMappingApplyIS(graph->l2gmap,to,&from);
723: }
724: /* compute local mirrors (if any) */
725: if (mirrors_found) {
726: PetscInt *local_indices,*global_indices;
727: /* get arrays of local and global indices */
728: PetscMalloc1(graph->nvtxs,&local_indices);
729: ISGetIndices(to,(const PetscInt**)&is_indices);
730: PetscMemcpy(local_indices,is_indices,graph->nvtxs*sizeof(PetscInt));
731: ISRestoreIndices(to,(const PetscInt**)&is_indices);
732: PetscMalloc1(graph->nvtxs,&global_indices);
733: ISGetIndices(from,(const PetscInt**)&is_indices);
734: PetscMemcpy(global_indices,is_indices,graph->nvtxs*sizeof(PetscInt));
735: ISRestoreIndices(from,(const PetscInt**)&is_indices);
736: /* allocate space for mirrors */
737: PetscMalloc2(graph->nvtxs,&graph->mirrors,graph->nvtxs,&graph->mirrors_set);
738: PetscMemzero(graph->mirrors,graph->nvtxs*sizeof(PetscInt));
739: graph->mirrors_set[0] = 0;
741: k=0;
742: for (i=0;i<n_shared[0];i++) {
743: j=shared[0][i];
744: if (graph->count[j] > 1) {
745: graph->mirrors[j]++;
746: k++;
747: }
748: }
749: /* allocate space for set of mirrors */
750: PetscMalloc1(k,&graph->mirrors_set[0]);
751: for (i=1;i<graph->nvtxs;i++)
752: graph->mirrors_set[i]=graph->mirrors_set[i-1]+graph->mirrors[i-1];
754: /* fill arrays */
755: PetscMemzero(graph->mirrors,graph->nvtxs*sizeof(PetscInt));
756: for (j=0;j<n_shared[0];j++) {
757: i=shared[0][j];
758: if (graph->count[i] > 1)
759: graph->mirrors_set[i][graph->mirrors[i]++]=global_indices[i];
760: }
761: PetscSortIntWithArray(graph->nvtxs,global_indices,local_indices);
762: for (i=0;i<graph->nvtxs;i++) {
763: if (graph->mirrors[i] > 0) {
764: PetscFindInt(graph->mirrors_set[i][0],graph->nvtxs,global_indices,&k);
765: j = global_indices[k];
766: while ( k > 0 && global_indices[k-1] == j) k--;
767: for (j=0;j<graph->mirrors[i];j++) {
768: graph->mirrors_set[i][j]=local_indices[k+j];
769: }
770: PetscSortInt(graph->mirrors[i],graph->mirrors_set[i]);
771: }
772: }
773: PetscFree(local_indices);
774: PetscFree(global_indices);
775: }
776: PetscMemzero(graph->count,graph->nvtxs*sizeof(*graph->count));
777: ISDestroy(&to);
778: ISDestroy(&from);
780: /* Count total number of neigh per node */
781: k=0;
782: for (i=1;i<n_neigh;i++) {
783: k += n_shared[i];
784: for (j=0;j<n_shared[i];j++) {
785: graph->count[shared[i][j]] += 1;
786: }
787: }
788: /* Allocate space for storing the set of neighbours for each node */
789: if (graph->nvtxs) {
790: PetscMalloc1(k,&graph->neighbours_set[0]);
791: }
792: for (i=1;i<graph->nvtxs;i++) { /* dont count myself */
793: graph->neighbours_set[i]=graph->neighbours_set[i-1]+graph->count[i-1];
794: }
795: /* Get information for sharing subdomains */
796: PetscMemzero(graph->count,graph->nvtxs*sizeof(*graph->count));
797: for (i=1;i<n_neigh;i++) { /* dont count myself */
798: s = n_shared[i];
799: for (j=0;j<s;j++) {
800: k = shared[i][j];
801: graph->neighbours_set[k][graph->count[k]] = neigh[i];
802: graph->count[k] += 1;
803: }
804: }
805: /* sort set of sharing subdomains */
806: for (i=0;i<graph->nvtxs;i++) {
807: PetscSortRemoveDupsInt(&graph->count[i],graph->neighbours_set[i]);
808: }
809: /* free memory allocated by ISLocalToGlobalMappingGetInfo */
810: ISLocalToGlobalMappingRestoreInfo(graph->l2gmap,&n_neigh,&neigh,&n_shared,&shared);
812: /*
813: Get info for dofs splitting
814: User can specify just a subset; an additional field is considered as a complementary field
815: */
816: for (i=0;i<graph->nvtxs;i++) {
817: graph->which_dof[i] = n_ISForDofs; /* by default a dof belongs to the complement set */
818: }
819: if (n_ISForDofs) {
820: VecSet(local_vec,-1.0);
821: }
822: for (i=0;i<n_ISForDofs;i++) {
823: VecGetArray(local_vec,&array);
824: ISGetLocalSize(ISForDofs[i],&is_size);
825: ISGetIndices(ISForDofs[i],(const PetscInt**)&is_indices);
826: for (j=0;j<is_size;j++) {
827: if (is_indices[j] > -1 && is_indices[j] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */
828: graph->which_dof[is_indices[j]] = i;
829: array[is_indices[j]] = 1.*i;
830: }
831: }
832: ISRestoreIndices(ISForDofs[i],(const PetscInt**)&is_indices);
833: VecRestoreArray(local_vec,&array);
834: }
835: /* Check consistency among neighbours */
836: if (n_ISForDofs) {
837: VecScatterBegin(scatter_ctx,local_vec,global_vec,INSERT_VALUES,SCATTER_REVERSE);
838: VecScatterEnd(scatter_ctx,local_vec,global_vec,INSERT_VALUES,SCATTER_REVERSE);
839: VecScatterBegin(scatter_ctx,global_vec,local_vec2,INSERT_VALUES,SCATTER_FORWARD);
840: VecScatterEnd(scatter_ctx,global_vec,local_vec2,INSERT_VALUES,SCATTER_FORWARD);
841: VecGetArray(local_vec,&array);
842: VecGetArray(local_vec2,&array2);
843: for (i=0;i<graph->nvtxs;i++){
844: PetscInt field1,field2;
846: field1 = (PetscInt)PetscRealPart(array[i]);
847: field2 = (PetscInt)PetscRealPart(array2[i]);
848: if (field1 != field2) {
849: SETERRQ3(comm,PETSC_ERR_USER,"Local node %d have been assigned two different field ids %d and %d at the same time\n",i,field1,field2);
850: }
851: }
852: VecRestoreArray(local_vec,&array);
853: VecRestoreArray(local_vec2,&array2);
854: }
855: if (neumann_is || dirichlet_is) {
856: /* Take into account Neumann nodes */
857: VecSet(local_vec,0.0);
858: VecSet(local_vec2,0.0);
859: if (neumann_is) {
860: VecGetArray(local_vec,&array);
861: ISGetLocalSize(neumann_is,&is_size);
862: ISGetIndices(neumann_is,(const PetscInt**)&is_indices);
863: for (i=0;i<is_size;i++) {
864: if (is_indices[i] > -1 && is_indices[i] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */
865: array[is_indices[i]] = 1.0;
866: }
867: }
868: ISRestoreIndices(neumann_is,(const PetscInt**)&is_indices);
869: VecRestoreArray(local_vec,&array);
870: }
871: /* Neumann nodes: impose consistency among neighbours */
872: VecSet(global_vec,0.0);
873: VecScatterBegin(scatter_ctx,local_vec,global_vec,ADD_VALUES,SCATTER_REVERSE);
874: VecScatterEnd(scatter_ctx,local_vec,global_vec,ADD_VALUES,SCATTER_REVERSE);
875: VecScatterBegin(scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_FORWARD);
876: VecScatterEnd(scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_FORWARD);
877: VecGetArray(local_vec,&array);
878: for (i=0;i<graph->nvtxs;i++) {
879: if (PetscRealPart(array[i]) > 0.1) {
880: graph->special_dof[i] = PCBDDCGRAPH_NEUMANN_MARK;
881: }
882: }
883: VecRestoreArray(local_vec,&array);
884: /* Take into account Dirichlet nodes */
885: VecSet(local_vec2,0.0);
886: if (dirichlet_is) {
887: VecGetArray(local_vec,&array);
888: VecGetArray(local_vec2,&array2);
889: ISGetLocalSize(dirichlet_is,&is_size);
890: ISGetIndices(dirichlet_is,(const PetscInt**)&is_indices);
891: for (i=0;i<is_size;i++){
892: if (is_indices[i] > -1 && is_indices[i] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */
893: k = is_indices[i];
894: if (graph->count[k] && !PetscBTLookup(graph->touched,k)) {
895: if (PetscRealPart(array[k]) > 0.1) {
896: SETERRQ1(comm,PETSC_ERR_USER,"BDDC cannot have boundary nodes which are marked Neumann and Dirichlet at the same time! Local node %d is wrong\n",k);
897: }
898: array2[k] = 1.0;
899: }
900: }
901: }
902: ISRestoreIndices(dirichlet_is,(const PetscInt**)&is_indices);
903: VecRestoreArray(local_vec,&array);
904: VecRestoreArray(local_vec2,&array2);
905: }
906: /* Dirichlet nodes: impose consistency among neighbours */
907: VecSet(global_vec,0.0);
908: VecScatterBegin(scatter_ctx,local_vec2,global_vec,ADD_VALUES,SCATTER_REVERSE);
909: VecScatterEnd(scatter_ctx,local_vec2,global_vec,ADD_VALUES,SCATTER_REVERSE);
910: VecScatterBegin(scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_FORWARD);
911: VecScatterEnd(scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_FORWARD);
912: VecGetArray(local_vec,&array);
913: for (i=0;i<graph->nvtxs;i++) {
914: if (PetscRealPart(array[i]) > 0.1) {
915: PetscBTSet(graph->touched,i);
916: graph->subset[i] = 0; /* dirichlet nodes treated as internal -> is it ok? */
917: graph->special_dof[i] = PCBDDCGRAPH_DIRICHLET_MARK;
918: }
919: }
920: VecRestoreArray(local_vec,&array);
921: }
922: /* mark local periodic nodes (if any) and adapt CSR graph (if any) */
923: if (graph->mirrors) {
924: for (i=0;i<graph->nvtxs;i++)
925: if (graph->mirrors[i])
926: graph->special_dof[i] = PCBDDCGRAPH_LOCAL_PERIODIC_MARK;
928: if (graph->xadj && graph->adjncy) {
929: PetscInt *new_xadj,*new_adjncy;
930: /* sort CSR graph */
931: for (i=0;i<graph->nvtxs;i++)
932: PetscSortInt(graph->xadj[i+1]-graph->xadj[i],&graph->adjncy[graph->xadj[i]]);
934: /* adapt local CSR graph in case of local periodicity */
935: k=0;
936: for (i=0;i<graph->nvtxs;i++)
937: for (j=graph->xadj[i];j<graph->xadj[i+1];j++)
938: k += graph->mirrors[graph->adjncy[j]];
940: PetscMalloc1(graph->nvtxs+1,&new_xadj);
941: PetscMalloc1(k+graph->xadj[graph->nvtxs],&new_adjncy);
942: new_xadj[0]=0;
943: for (i=0;i<graph->nvtxs;i++) {
944: k = graph->xadj[i+1]-graph->xadj[i];
945: PetscMemcpy(&new_adjncy[new_xadj[i]],&graph->adjncy[graph->xadj[i]],k*sizeof(PetscInt));
946: new_xadj[i+1]=new_xadj[i]+k;
947: for (j=graph->xadj[i];j<graph->xadj[i+1];j++) {
948: k = graph->mirrors[graph->adjncy[j]];
949: PetscMemcpy(&new_adjncy[new_xadj[i+1]],graph->mirrors_set[graph->adjncy[j]],k*sizeof(PetscInt));
950: new_xadj[i+1]+=k;
951: }
952: k = new_xadj[i+1]-new_xadj[i];
953: PetscSortRemoveDupsInt(&k,&new_adjncy[new_xadj[i]]);
954: new_xadj[i+1]=new_xadj[i]+k;
955: }
956: /* set new CSR into graph */
957: PetscFree(graph->xadj);
958: PetscFree(graph->adjncy);
959: graph->xadj = new_xadj;
960: graph->adjncy = new_adjncy;
961: }
962: }
964: /* mark special nodes (if any) -> each will become a single node equivalence class */
965: if (custom_primal_vertices) {
966: VecSet(local_vec,0.0);
967: VecGetArray(local_vec,&array);
968: ISGetLocalSize(custom_primal_vertices,&is_size);
969: ISGetIndices(custom_primal_vertices,(const PetscInt**)&is_indices);
970: for (i=0;i<is_size;i++){
971: if (is_indices[i] > -1 && is_indices[i] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */
972: array[is_indices[i]] = 1.0;
973: }
974: }
975: ISRestoreIndices(custom_primal_vertices,(const PetscInt**)&is_indices);
976: VecRestoreArray(local_vec,&array);
977: /* special nodes: impose consistency among neighbours */
978: VecSet(global_vec,0.0);
979: VecScatterBegin(scatter_ctx,local_vec,global_vec,ADD_VALUES,SCATTER_REVERSE);
980: VecScatterEnd(scatter_ctx,local_vec,global_vec,ADD_VALUES,SCATTER_REVERSE);
981: VecScatterBegin(scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_FORWARD);
982: VecScatterEnd(scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_FORWARD);
983: VecGetArray(local_vec,&array);
984: j = 0;
985: for (i=0;i<graph->nvtxs;i++) {
986: if (PetscRealPart(array[i]) > 0.1 && graph->special_dof[i] != PCBDDCGRAPH_DIRICHLET_MARK) {
987: graph->special_dof[i] = PCBDDCGRAPH_SPECIAL_MARK-j;
988: j++;
989: }
990: }
991: VecRestoreArray(local_vec,&array);
992: }
994: /* mark interior nodes as touched and belonging to partition number 0 */
995: for (i=0;i<graph->nvtxs;i++) {
996: if (!graph->count[i]) {
997: PetscBTSet(graph->touched,i);
998: graph->subset[i] = 0;
999: }
1000: }
1001: /* init graph structure and compute default subsets */
1002: nodes_touched=0;
1003: for (i=0;i<graph->nvtxs;i++) {
1004: if (PetscBTLookup(graph->touched,i)) {
1005: nodes_touched++;
1006: }
1007: }
1008: i = 0;
1009: graph->ncc = 0;
1010: total_counts = 0;
1012: /* allocated space for queues */
1013: MPI_Comm_size(comm,&size);
1014: if (size == 1) {
1015: PetscMalloc2(graph->nvtxs+1,&graph->cptr,graph->nvtxs,&graph->queue);
1016: } else {
1017: PetscInt nused = graph->nvtxs - nodes_touched;
1018: PetscMalloc2(nused+1,&graph->cptr,nused,&graph->queue);
1019: }
1021: while (nodes_touched<graph->nvtxs) {
1022: /* find first untouched node in local ordering */
1023: while (PetscBTLookup(graph->touched,i)) {
1024: i++;
1025: }
1026: PetscBTSet(graph->touched,i);
1027: graph->subset[i] = graph->ncc+1;
1028: graph->cptr[graph->ncc] = total_counts;
1029: graph->queue[total_counts] = i;
1030: total_counts++;
1031: nodes_touched++;
1032: /* now find all other nodes having the same set of sharing subdomains */
1033: for (j=i+1;j<graph->nvtxs;j++) {
1034: /* check for same number of sharing subdomains, dof number and same special mark */
1035: 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]) {
1036: /* check for same set of sharing subdomains */
1037: same_set=PETSC_TRUE;
1038: for (k=0;k<graph->count[j];k++){
1039: if (graph->neighbours_set[i][k]!=graph->neighbours_set[j][k]) {
1040: same_set=PETSC_FALSE;
1041: }
1042: }
1043: /* I found a friend of mine */
1044: if (same_set) {
1045: graph->subset[j]=graph->ncc+1;
1046: PetscBTSet(graph->touched,j);
1047: nodes_touched++;
1048: graph->queue[total_counts] = j;
1049: total_counts++;
1050: }
1051: }
1052: }
1053: graph->ncc++;
1054: }
1055: /* set default number of subsets (at this point no info on csr graph has been taken into account, so n_subsets = ncc */
1056: graph->n_subsets = graph->ncc;
1057: PetscMalloc1(graph->n_subsets,&graph->subset_ncc);
1058: for (i=0;i<graph->n_subsets;i++) {
1059: graph->subset_ncc[i] = 1;
1060: }
1061: /* final pointer */
1062: graph->cptr[graph->ncc] = total_counts;
1064: /* save information on subsets (needed if we have to adapt the connected components later) */
1065: /* For consistency reasons (among neighbours), I need to sort (by global ordering) each subset */
1066: /* Get a reference node (min index in global ordering) for each subset for tagging messages */
1067: PetscMalloc1(graph->ncc,&graph->subset_ref_node);
1068: PetscMalloc1(graph->cptr[graph->ncc],&queue_global);
1069: ISLocalToGlobalMappingApply(graph->l2gmap,graph->cptr[graph->ncc],graph->queue,queue_global);
1070: for (j=0;j<graph->ncc;j++) {
1071: PetscSortIntWithArray(graph->cptr[j+1]-graph->cptr[j],&queue_global[graph->cptr[j]],&graph->queue[graph->cptr[j]]);
1072: graph->subset_ref_node[j] = graph->queue[graph->cptr[j]];
1073: }
1074: PetscFree(queue_global);
1075: graph->queue_sorted = PETSC_TRUE;
1076: if (graph->ncc) {
1077: PetscMalloc2(graph->ncc,&graph->subsets_size,graph->ncc,&graph->subsets);
1078: PetscMalloc1(graph->cptr[graph->ncc],&graph->subsets[0]);
1079: PetscMemzero(graph->subsets[0],graph->cptr[graph->ncc]*sizeof(PetscInt));
1080: for (j=1;j<graph->ncc;j++) {
1081: graph->subsets_size[j-1] = graph->cptr[j] - graph->cptr[j-1];
1082: graph->subsets[j] = graph->subsets[j-1] + graph->subsets_size[j-1];
1083: }
1084: graph->subsets_size[graph->ncc-1] = graph->cptr[graph->ncc] - graph->cptr[graph->ncc-1];
1085: for (j=0;j<graph->ncc;j++) {
1086: PetscMemcpy(graph->subsets[j],&graph->queue[graph->cptr[j]],graph->subsets_size[j]*sizeof(PetscInt));
1087: }
1088: }
1089: /* renumber reference nodes */
1090: ISCreateGeneral(PetscObjectComm((PetscObject)(graph->l2gmap)),graph->ncc,graph->subset_ref_node,PETSC_COPY_VALUES,&subset_n);
1091: ISLocalToGlobalMappingApplyIS(graph->l2gmap,subset_n,&subset);
1092: ISDestroy(&subset_n);
1093: PCBDDCSubsetNumbering(subset,NULL,NULL,&subset_n);
1094: ISDestroy(&subset);
1095: ISGetLocalSize(subset_n,&k);
1096: if (k != graph->ncc) {
1097: SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Invalid size of new subset! %d != %d",k,graph->ncc);
1098: }
1099: ISGetIndices(subset_n,&is_indices);
1100: PetscMemcpy(graph->subset_ref_node,is_indices,graph->ncc*sizeof(PetscInt));
1101: ISRestoreIndices(subset_n,&is_indices);
1102: ISDestroy(&subset_n);
1104: /* free workspace */
1105: VecDestroy(&local_vec);
1106: VecDestroy(&local_vec2);
1107: VecDestroy(&global_vec);
1108: VecScatterDestroy(&scatter_ctx);
1109: return(0);
1110: }
1114: PetscErrorCode PCBDDCGraphResetCSR(PCBDDCGraph graph)
1115: {
1119: PetscFree(graph->xadj);
1120: PetscFree(graph->adjncy);
1121: graph->nvtxs_csr = 0;
1122: return(0);
1123: }
1127: PetscErrorCode PCBDDCGraphReset(PCBDDCGraph graph)
1128: {
1132: ISLocalToGlobalMappingDestroy(&graph->l2gmap);
1133: PetscFree(graph->subset_ncc);
1134: PetscFree(graph->subset_ref_node);
1135: if (graph->nvtxs) {
1136: PetscFree(graph->neighbours_set[0]);
1137: }
1138: PetscBTDestroy(&graph->touched);
1139: PetscFree5(graph->count,
1140: graph->neighbours_set,
1141: graph->subset,
1142: graph->which_dof,
1143: graph->special_dof);
1144: PetscFree2(graph->cptr,graph->queue);
1145: if (graph->mirrors) {
1146: PetscFree(graph->mirrors_set[0]);
1147: }
1148: PetscFree2(graph->mirrors,graph->mirrors_set);
1149: if (graph->subsets) {
1150: PetscFree(graph->subsets[0]);
1151: }
1152: PetscFree2(graph->subsets_size,graph->subsets);
1153: ISDestroy(&graph->dirdofs);
1154: ISDestroy(&graph->dirdofsB);
1155: graph->has_dirichlet = PETSC_FALSE;
1156: graph->nvtxs = 0;
1157: graph->nvtxs_global = 0;
1158: graph->n_subsets = 0;
1159: graph->custom_minimal_size = 1;
1160: return(0);
1161: }
1165: PetscErrorCode PCBDDCGraphInit(PCBDDCGraph graph, ISLocalToGlobalMapping l2gmap, PetscInt N)
1166: {
1167: PetscInt n;
1174: /* raise an error if already allocated */
1175: if (graph->nvtxs_global) {
1176: SETERRQ(PetscObjectComm((PetscObject)l2gmap),PETSC_ERR_PLIB,"BDDCGraph already initialized");
1177: }
1178: /* set number of vertices */
1179: PetscObjectReference((PetscObject)l2gmap);
1180: graph->l2gmap = l2gmap;
1181: ISLocalToGlobalMappingGetSize(l2gmap,&n);
1182: graph->nvtxs = n;
1183: graph->nvtxs_global = N;
1184: /* allocate used space */
1185: PetscBTCreate(graph->nvtxs,&graph->touched);
1186: PetscMalloc5(graph->nvtxs,&graph->count,
1187: graph->nvtxs,&graph->neighbours_set,
1188: graph->nvtxs,&graph->subset,
1189: graph->nvtxs,&graph->which_dof,
1190: graph->nvtxs,&graph->special_dof);
1191: /* zeroes memory */
1192: PetscMemzero(graph->count,graph->nvtxs*sizeof(PetscInt));
1193: PetscMemzero(graph->subset,graph->nvtxs*sizeof(PetscInt));
1194: /* use -1 as a default value for which_dof array */
1195: for (n=0;n<graph->nvtxs;n++) graph->which_dof[n] = -1;
1196: PetscMemzero(graph->special_dof,graph->nvtxs*sizeof(PetscInt));
1197: /* zeroes first pointer to neighbour set */
1198: if (graph->nvtxs) {
1199: graph->neighbours_set[0] = 0;
1200: }
1201: /* zeroes workspace for values of ncc */
1202: graph->subset_ncc = 0;
1203: graph->subset_ref_node = 0;
1204: return(0);
1205: }
1209: PetscErrorCode PCBDDCGraphDestroy(PCBDDCGraph* graph)
1210: {
1214: PCBDDCGraphReset(*graph);
1215: PetscFree(*graph);
1216: return(0);
1217: }
1221: PetscErrorCode PCBDDCGraphCreate(PCBDDCGraph *graph)
1222: {
1223: PCBDDCGraph new_graph;
1227: PetscNew(&new_graph);
1228: new_graph->custom_minimal_size = 1;
1229: *graph = new_graph;
1230: return(0);
1231: }