Actual source code: bddcgraph.c
petsc-3.5.4 2015-05-23
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 PCBDDCGraphASCIIView(PCBDDCGraph graph, PetscInt verbosity_level, PetscViewer viewer)
8: {
9: PetscInt i,j,tabs;
10: PetscInt* queue_in_global_numbering;
14: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
15: PetscViewerASCIIGetTab(viewer,&tabs);
16: PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");
17: PetscViewerFlush(viewer);
18: PetscViewerASCIISynchronizedPrintf(viewer,"Local BDDC graph for subdomain %04d\n",PetscGlobalRank);
19: PetscViewerASCIISynchronizedPrintf(viewer,"Number of vertices %d\n",graph->nvtxs);
20: if (verbosity_level > 1) {
21: for (i=0;i<graph->nvtxs;i++) {
22: PetscViewerASCIISynchronizedPrintf(viewer,"%d:\n",i);
23: PetscViewerASCIISynchronizedPrintf(viewer," which_dof: %d\n",graph->which_dof[i]);
24: PetscViewerASCIISynchronizedPrintf(viewer," special_dof: %d\n",graph->special_dof[i]);
25: PetscViewerASCIISynchronizedPrintf(viewer," neighbours: %d\n",graph->count[i]);
26: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
27: if (graph->count[i]) {
28: PetscViewerASCIISynchronizedPrintf(viewer," set of neighbours:");
29: for (j=0;j<graph->count[i];j++) {
30: PetscViewerASCIISynchronizedPrintf(viewer," %d",graph->neighbours_set[i][j]);
31: }
32: PetscViewerASCIISynchronizedPrintf(viewer,"\n");
33: }
34: PetscViewerASCIISetTab(viewer,tabs);
35: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
36: if (graph->mirrors) {
37: PetscViewerASCIISynchronizedPrintf(viewer," mirrors: %d\n",graph->mirrors[i]);
38: if (graph->mirrors[i]) {
39: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
40: PetscViewerASCIISynchronizedPrintf(viewer," set of mirrors:");
41: for (j=0;j<graph->mirrors[i];j++) {
42: PetscViewerASCIISynchronizedPrintf(viewer," %d",graph->mirrors_set[i][j]);
43: }
44: PetscViewerASCIISynchronizedPrintf(viewer,"\n");
45: PetscViewerASCIISetTab(viewer,tabs);
46: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
47: }
48: }
49: if (verbosity_level > 2) {
50: if (graph->xadj && graph->adjncy) {
51: PetscViewerASCIISynchronizedPrintf(viewer," local adj list:");
52: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
53: for (j=graph->xadj[i];j<graph->xadj[i+1];j++) {
54: PetscViewerASCIISynchronizedPrintf(viewer," %d",graph->adjncy[j]);
55: }
56: PetscViewerASCIISynchronizedPrintf(viewer,"\n");
57: PetscViewerASCIISetTab(viewer,tabs);
58: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
59: }
60: }
61: PetscViewerASCIISynchronizedPrintf(viewer," interface subset id: %d\n",graph->subset[i]);
62: if (graph->subset[i] && graph->subset_ncc) {
63: PetscViewerASCIISynchronizedPrintf(viewer," ncc for subset: %d\n",graph->subset_ncc[graph->subset[i]-1]);
64: }
65: }
66: }
67: PetscViewerASCIISynchronizedPrintf(viewer,"Total number of connected components %d\n",graph->ncc);
68: PetscMalloc1(graph->cptr[graph->ncc],&queue_in_global_numbering);
69: ISLocalToGlobalMappingApply(graph->l2gmap,graph->cptr[graph->ncc],graph->queue,queue_in_global_numbering);
70: for (i=0;i<graph->ncc;i++) {
71: PetscInt node_num=graph->queue[graph->cptr[i]];
72: PetscBool printcc = PETSC_FALSE;
73: PetscViewerASCIISynchronizedPrintf(viewer," %d (neighs:",i);
74: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
75: for (j=0;j<graph->count[node_num];j++) {
76: PetscViewerASCIISynchronizedPrintf(viewer," %d",graph->neighbours_set[node_num][j]);
77: }
78: PetscViewerASCIISynchronizedPrintf(viewer,"):");
79: if (verbosity_level > 1) {
80: printcc = PETSC_TRUE;
81: } else if (graph->count[node_num] > 1 || (graph->count[node_num] == 1 && graph->special_dof[node_num] == PCBDDCGRAPH_NEUMANN_MARK)) {
82: printcc = PETSC_TRUE;
83: }
84: if (printcc) {
85: for (j=graph->cptr[i];j<graph->cptr[i+1];j++) {
86: PetscViewerASCIISynchronizedPrintf(viewer," %d (%d)",graph->queue[j],queue_in_global_numbering[j]);
87: }
88: }
89: PetscViewerASCIISynchronizedPrintf(viewer,"\n");
90: PetscViewerASCIISetTab(viewer,tabs);
91: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
92: }
93: PetscFree(queue_in_global_numbering);
94: if (graph->custom_minimal_size > 1 && verbosity_level > 1) {
95: PetscViewerASCIISynchronizedPrintf(viewer,"Custom minimal size %d\n",graph->custom_minimal_size);
96: }
97: PetscViewerFlush(viewer);
98: return(0);
99: }
103: PetscErrorCode PCBDDCGraphGetCandidatesIS(PCBDDCGraph graph, PetscBool use_faces, PetscBool use_edges, PetscBool use_vertices, PetscInt *n_faces, IS *FacesIS[], PetscInt *n_edges, IS *EdgesIS[], IS *VerticesIS)
104: {
105: IS *ISForFaces,*ISForEdges,ISForVertices;
106: PetscInt i,j,nfc,nec,nvc,*idx;
107: PetscBool twodim_flag;
111: /* loop on ccs to evalute number of faces, edges and vertices */
112: nfc = 0;
113: nec = 0;
114: nvc = 0;
115: twodim_flag = PETSC_FALSE;
116: for (i=0;i<graph->ncc;i++) {
117: PetscInt repdof = graph->queue[graph->cptr[i]];
118: if (graph->cptr[i+1]-graph->cptr[i] > graph->custom_minimal_size) {
119: if (graph->count[repdof] == 1 && graph->special_dof[repdof] != PCBDDCGRAPH_NEUMANN_MARK) {
120: nfc++;
121: } else { /* note that nec will be zero in 2d */
122: nec++;
123: }
124: } else {
125: if (graph->count[repdof] > 1 || (graph->count[repdof] == 1 && graph->special_dof[repdof] == PCBDDCGRAPH_NEUMANN_MARK)) {
126: nvc += graph->cptr[i+1]-graph->cptr[i];
127: }
128: }
129: }
130: j=0;
131: MPI_Allreduce(&nec,&j,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)graph->l2gmap));
132: if (!j) { /* we are in a 2D case -> no faces, only edges */
133: nec = nfc;
134: nfc = 0;
135: twodim_flag = PETSC_TRUE;
136: }
137: /* allocate IS arrays for faces, edges. Vertices need a single index set. */
138: ISForFaces = 0;
139: ISForEdges = 0;
140: ISForVertices = 0;
141: if (use_faces && nfc) {
142: PetscMalloc1(nfc,&ISForFaces);
143: }
144: if (use_edges && nec) {
145: PetscMalloc1(nec,&ISForEdges);
146: }
147: if (use_vertices && nvc) {
148: PetscMalloc1(nvc,&idx);
149: }
150: /* loop on ccs to compute index sets for faces and edges */
151: nfc = 0;
152: nec = 0;
153: for (i=0;i<graph->ncc;i++) {
154: PetscInt repdof = graph->queue[graph->cptr[i]];
155: if (graph->cptr[i+1]-graph->cptr[i] > graph->custom_minimal_size) {
156: if (graph->count[repdof] == 1 && graph->special_dof[repdof] != PCBDDCGRAPH_NEUMANN_MARK) {
157: if (twodim_flag) {
158: if (use_edges) {
159: ISCreateGeneral(PETSC_COMM_SELF,graph->cptr[i+1]-graph->cptr[i],&graph->queue[graph->cptr[i]],PETSC_COPY_VALUES,&ISForEdges[nec]);
160: nec++;
161: }
162: } else {
163: if (use_faces) {
164: ISCreateGeneral(PETSC_COMM_SELF,graph->cptr[i+1]-graph->cptr[i],&graph->queue[graph->cptr[i]],PETSC_COPY_VALUES,&ISForFaces[nfc]);
165: nfc++;
166: }
167: }
168: } else {
169: if (use_edges) {
170: ISCreateGeneral(PETSC_COMM_SELF,graph->cptr[i+1]-graph->cptr[i],&graph->queue[graph->cptr[i]],PETSC_COPY_VALUES,&ISForEdges[nec]);
171: nec++;
172: }
173: }
174: }
175: }
176: /* index set for vertices */
177: if (use_vertices && nvc) {
178: nvc = 0;
179: for (i=0;i<graph->ncc;i++) {
180: if (graph->cptr[i+1]-graph->cptr[i] <= graph->custom_minimal_size) {
181: PetscInt repdof = graph->queue[graph->cptr[i]];
182: if (graph->count[repdof] > 1 || (graph->count[repdof] == 1 && graph->special_dof[repdof] == PCBDDCGRAPH_NEUMANN_MARK)) {
183: for (j=graph->cptr[i];j<graph->cptr[i+1];j++) {
184: idx[nvc]=graph->queue[j];
185: nvc++;
186: }
187: }
188: }
189: }
190: /* sort vertex set (by local ordering) */
191: PetscSortInt(nvc,idx);
192: ISCreateGeneral(PETSC_COMM_SELF,nvc,idx,PETSC_OWN_POINTER,&ISForVertices);
193: }
194: /* get back info */
195: *n_faces = nfc;
196: *FacesIS = ISForFaces;
197: *n_edges = nec;
198: *EdgesIS = ISForEdges;
199: *VerticesIS = ISForVertices;
200: return(0);
201: }
205: PetscErrorCode PCBDDCGraphComputeConnectedComponents(PCBDDCGraph graph)
206: {
207: PetscInt i,adapt_interface,adapt_interface_reduced;
208: MPI_Comm interface_comm;
212: /* compute connected components locally */
213: PetscObjectGetComm((PetscObject)(graph->l2gmap),&interface_comm);
214: PCBDDCGraphComputeConnectedComponentsLocal(graph);
215: /* check consistency of connected components among neighbouring subdomains -> it adapt them in case it is needed */
216: adapt_interface = 0;
217: adapt_interface_reduced = 0;
218: for (i=0;i<graph->n_subsets;i++) {
219: /* We are not sure that on a given subset of the local interface,
220: with two connected components, the latters be the same among sharing subdomains */
221: if (graph->subset_ncc[i] > 1) {
222: adapt_interface=1;
223: break;
224: }
225: }
226: MPI_Allreduce(&adapt_interface,&adapt_interface_reduced,1,MPIU_INT,MPI_LOR,interface_comm);
228: if (graph->n_subsets && adapt_interface_reduced) {
229: MPI_Request *send_requests;
230: MPI_Request *recv_requests;
231: PetscInt *aux_new_xadj,*new_xadj,*new_adjncy,**temp_buffer;
232: PetscInt *old_xadj,*old_adjncy;
233: PetscInt j,k,s,sum_requests,buffer_size,size_of_recv,temp_buffer_size;
234: PetscMPIInt rank,neigh,tag,mpi_buffer_size;
235: PetscInt *cum_recv_counts,*subset_to_nodes_indices,*recv_buffer_subset,*nodes_to_temp_buffer_indices;
236: PetscInt *send_buffer,*recv_buffer,*queue_in_global_numbering,*sizes_of_sends,*add_to_subset;
237: PetscInt start_of_recv,start_of_send,size_of_send,global_subset_counter,ins_val;
238: PetscBool *subset_cc_adapt,same_set;
240: /* Retrict adjacency graph using information from previously computed connected components */
241: PetscMalloc1(graph->nvtxs,&aux_new_xadj);
242: for (i=0;i<graph->nvtxs;i++) {
243: aux_new_xadj[i]=1;
244: }
245: for (i=0;i<graph->ncc;i++) {
246: k = graph->cptr[i+1]-graph->cptr[i];
247: for (j=0;j<k;j++) {
248: aux_new_xadj[graph->queue[graph->cptr[i]+j]]=k;
249: }
250: }
251: j = 0;
252: for (i=0;i<graph->nvtxs;i++) {
253: j += aux_new_xadj[i];
254: }
255: PetscMalloc1((graph->nvtxs+1),&new_xadj);
256: PetscMalloc1(j,&new_adjncy);
257: new_xadj[0]=0;
258: for (i=0;i<graph->nvtxs;i++) {
259: new_xadj[i+1]=new_xadj[i]+aux_new_xadj[i];
260: if (aux_new_xadj[i]==1) {
261: new_adjncy[new_xadj[i]]=i;
262: }
263: }
264: PetscFree(aux_new_xadj);
265: for (i=0;i<graph->ncc;i++) {
266: k = graph->cptr[i+1]-graph->cptr[i];
267: for (j=0;j<k;j++) {
268: PetscMemcpy(&new_adjncy[new_xadj[graph->queue[graph->cptr[i]+j]]],&graph->queue[graph->cptr[i]],k*sizeof(PetscInt));
269: }
270: }
271: /* set temporarly new CSR into graph */
272: old_xadj = graph->xadj;
273: old_adjncy = graph->adjncy;
274: graph->xadj = new_xadj;
275: graph->adjncy = new_adjncy;
276: /* allocate some space */
277: MPI_Comm_rank(interface_comm,&rank);
278: PetscMalloc1((graph->n_subsets+1),&cum_recv_counts);
279: PetscMemzero(cum_recv_counts,(graph->n_subsets+1)*sizeof(*cum_recv_counts));
280: PetscMalloc1(graph->n_subsets,&subset_to_nodes_indices);
281: /* first count how many neighbours per connected component I will receive from */
282: cum_recv_counts[0]=0;
283: for (i=1;i<graph->n_subsets+1;i++) {
284: j = 0;
285: while (graph->subset[j] != i) {
286: j++;
287: }
288: subset_to_nodes_indices[i-1]=j;
289: /* We don't want sends/recvs_to/from_self -> here I don't count myself */
290: cum_recv_counts[i]=cum_recv_counts[i-1]+graph->count[j];
291: }
292: PetscMalloc1(2*cum_recv_counts[graph->n_subsets],&recv_buffer_subset);
293: PetscMalloc1(cum_recv_counts[graph->n_subsets],&send_requests);
294: PetscMalloc1(cum_recv_counts[graph->n_subsets],&recv_requests);
295: for (i=0;i<cum_recv_counts[graph->n_subsets];i++) {
296: send_requests[i]=MPI_REQUEST_NULL;
297: recv_requests[i]=MPI_REQUEST_NULL;
298: }
299: /* exchange with my neighbours the number of my connected components on the shared interface */
300: sum_requests = 0;
301: for (i=0;i<graph->n_subsets;i++) {
302: j = subset_to_nodes_indices[i];
303: PetscMPIIntCast(graph->subset_ref_node[i],&tag);
304: for (k=0;k<graph->count[j];k++) {
305: PetscMPIIntCast(graph->neighbours_set[j][k],&neigh);
306: MPI_Isend(&graph->subset_ncc[i],1,MPIU_INT,neigh,tag,interface_comm,&send_requests[sum_requests]);
307: MPI_Irecv(&recv_buffer_subset[sum_requests],1,MPIU_INT,neigh,tag,interface_comm,&recv_requests[sum_requests]);
308: sum_requests++;
309: }
310: }
311: MPI_Waitall(sum_requests,recv_requests,MPI_STATUSES_IGNORE);
312: MPI_Waitall(sum_requests,send_requests,MPI_STATUSES_IGNORE);
313: /* determine the connected component I need to adapt */
314: PetscMalloc1(graph->n_subsets,&subset_cc_adapt);
315: PetscMemzero(subset_cc_adapt,graph->n_subsets*sizeof(*subset_cc_adapt));
316: for (i=0;i<graph->n_subsets;i++) {
317: for (j=cum_recv_counts[i];j<cum_recv_counts[i+1];j++){
318: /* The first condition is natural (someone has a different number of ccs than me), the second one is just to be safe */
319: if (graph->subset_ncc[i] != recv_buffer_subset[j] || graph->subset_ncc[i] > 1) {
320: subset_cc_adapt[i] = PETSC_TRUE;
321: break;
322: }
323: }
324: }
325: buffer_size = 0;
326: for (i=0;i<graph->n_subsets;i++) {
327: if (subset_cc_adapt[i]) {
328: for (j=i;j<graph->ncc;j++) {
329: if (graph->subset[graph->queue[graph->cptr[j]]] == i+1) { /* WARNING -> subset values goes from 1 to graph->n_subsets included */
330: buffer_size += 1 + graph->cptr[j+1]-graph->cptr[j];
331: }
332: }
333: }
334: }
335: PetscMalloc1(buffer_size,&send_buffer);
336: /* now get from neighbours their ccs (in global numbering) and adapt them (in case it is needed) */
337: PetscMalloc1(graph->cptr[graph->ncc],&queue_in_global_numbering);
338: ISLocalToGlobalMappingApply(graph->l2gmap,graph->cptr[graph->ncc],graph->queue,queue_in_global_numbering);
339: /* determine how much data to send (size of each queue plus the global indices) and communicate it to neighbours */
340: PetscMalloc1(graph->n_subsets,&sizes_of_sends);
341: PetscMemzero(sizes_of_sends,graph->n_subsets*sizeof(*sizes_of_sends));
342: sum_requests = 0;
343: start_of_send = 0;
344: start_of_recv = cum_recv_counts[graph->n_subsets];
345: for (i=0;i<graph->n_subsets;i++) {
346: if (subset_cc_adapt[i]) {
347: size_of_send = 0;
348: for (j=i;j<graph->ncc;j++) {
349: if (graph->subset[graph->queue[graph->cptr[j]]] == i+1) { /* WARNING -> subset values goes from 1 to graph->n_subsets included */
350: send_buffer[start_of_send+size_of_send]=graph->cptr[j+1]-graph->cptr[j];
351: size_of_send += 1;
352: PetscMemcpy(&send_buffer[start_of_send+size_of_send],
353: &queue_in_global_numbering[graph->cptr[j]],
354: (graph->cptr[j+1]-graph->cptr[j])*sizeof(*send_buffer));
355: size_of_send = size_of_send+graph->cptr[j+1]-graph->cptr[j];
356: }
357: }
358: j = subset_to_nodes_indices[i];
359: sizes_of_sends[i] = size_of_send;
360: PetscMPIIntCast(graph->subset_ref_node[i]+1,&tag);
361: for (k=0;k<graph->count[j];k++) {
362: PetscMPIIntCast(graph->neighbours_set[j][k],&neigh);
363: MPI_Isend(&sizes_of_sends[i],1,MPIU_INT,neigh,tag,interface_comm,&send_requests[sum_requests]);
364: MPI_Irecv(&recv_buffer_subset[sum_requests+start_of_recv],1,MPIU_INT,neigh,tag,interface_comm,&recv_requests[sum_requests]);
365: sum_requests++;
366: }
367: start_of_send += size_of_send;
368: }
369: }
370: MPI_Waitall(sum_requests,send_requests,MPI_STATUSES_IGNORE);
371: MPI_Waitall(sum_requests,recv_requests,MPI_STATUSES_IGNORE);
372: buffer_size = 0;
373: for (k=0;k<sum_requests;k++) {
374: buffer_size += recv_buffer_subset[start_of_recv+k];
375: }
376: PetscMalloc1(buffer_size,&recv_buffer);
377: /* now exchange the data */
378: start_of_recv = 0;
379: start_of_send = 0;
380: sum_requests = 0;
381: for (i=0;i<graph->n_subsets;i++) {
382: if (subset_cc_adapt[i]) {
383: size_of_send = sizes_of_sends[i];
384: j = subset_to_nodes_indices[i];
385: PetscMPIIntCast(graph->subset_ref_node[i]+2,&tag);
386: for (k=0;k<graph->count[j];k++) {
387: PetscMPIIntCast(graph->neighbours_set[j][k],&neigh);
388: PetscMPIIntCast(size_of_send,&mpi_buffer_size);
389: MPI_Isend(&send_buffer[start_of_send],mpi_buffer_size,MPIU_INT,neigh,tag,interface_comm,&send_requests[sum_requests]);
390: size_of_recv = recv_buffer_subset[cum_recv_counts[graph->n_subsets]+sum_requests];
391: PetscMPIIntCast(size_of_recv,&mpi_buffer_size);
392: MPI_Irecv(&recv_buffer[start_of_recv],mpi_buffer_size,MPIU_INT,neigh,tag,interface_comm,&recv_requests[sum_requests]);
393: start_of_recv += size_of_recv;
394: sum_requests++;
395: }
396: start_of_send += size_of_send;
397: }
398: }
399: MPI_Waitall(sum_requests,recv_requests,MPI_STATUSES_IGNORE);
400: MPI_Waitall(sum_requests,send_requests,MPI_STATUSES_IGNORE);
401: for (j=0;j<buffer_size;) {
402: ISGlobalToLocalMappingApply(graph->l2gmap,IS_GTOLM_MASK,recv_buffer[j],&recv_buffer[j+1],&recv_buffer[j],&recv_buffer[j+1]);
403: /* we need to adapt the output of GlobalToLocal mapping if there are mirrored nodes */
404: if (graph->mirrors) {
405: PetscBool mirrored_found=PETSC_FALSE;
406: for (k=0;k<recv_buffer[j];k++) {
407: if (graph->mirrors[recv_buffer[j+k+1]]) {
408: mirrored_found=PETSC_TRUE;
409: recv_buffer[j+k+1]=graph->mirrors_set[recv_buffer[j+k+1]][0];
410: }
411: }
412: if (mirrored_found) {
413: PetscSortInt(recv_buffer[j],&recv_buffer[j+1]);
414: k=0;
415: while (k<recv_buffer[j]) {
416: for (s=1;s<graph->mirrors[recv_buffer[j+1+k]];s++) {
417: recv_buffer[j+1+k+s] = graph->mirrors_set[recv_buffer[j+1+k]][s];
418: }
419: k+=graph->mirrors[recv_buffer[j+1+k]]+s;
420: }
421: }
422: }
423: k = recv_buffer[j]+1;
424: j += k;
425: }
426: sum_requests = cum_recv_counts[graph->n_subsets];
427: start_of_recv = 0;
428: PetscMalloc1(graph->nvtxs,&nodes_to_temp_buffer_indices);
429: global_subset_counter = 0;
430: for (i=0;i<graph->n_subsets;i++) {
431: if (subset_cc_adapt[i]) {
432: temp_buffer_size = 0;
433: /* find nodes on the shared interface we need to adapt */
434: for (j=0;j<graph->nvtxs;j++) {
435: if (graph->subset[j]==i+1) {
436: nodes_to_temp_buffer_indices[j] = temp_buffer_size;
437: temp_buffer_size++;
438: } else {
439: nodes_to_temp_buffer_indices[j] = -1;
440: }
441: }
442: /* allocate some temporary space */
443: PetscMalloc1(temp_buffer_size,&temp_buffer);
444: PetscMalloc1(temp_buffer_size*(cum_recv_counts[i+1]-cum_recv_counts[i]),&temp_buffer[0]);
445: PetscMemzero(temp_buffer[0],temp_buffer_size*(cum_recv_counts[i+1]-cum_recv_counts[i])*sizeof(PetscInt));
446: for (j=1;j<temp_buffer_size;j++) {
447: temp_buffer[j] = temp_buffer[j-1]+cum_recv_counts[i+1]-cum_recv_counts[i];
448: }
449: /* analyze contributions from neighbouring subdomains for i-th conn comp
450: temp buffer structure:
451: supposing part of the interface has dimension 5 (for example with global dofs 0,1,2,3,4)
452: 3 neighs procs with structured connected components:
453: neigh 0: [0 1 4], [2 3]; (2 connected components)
454: neigh 1: [0 1], [2 3 4]; (2 connected components)
455: neigh 2: [0 4], [1], [2 3]; (3 connected components)
456: tempbuffer (row-oriented) will be filled as:
457: [ 0, 0, 0;
458: 0, 0, 1;
459: 1, 1, 2;
460: 1, 1, 2;
461: 0, 1, 0; ];
462: This way we can simply find intersections of ccs among neighs.
463: For the example above, the graph->subset array will be modified to reproduce the following 4 connected components [0], [1], [2 3], [4];
464: */
465: for (j=0;j<cum_recv_counts[i+1]-cum_recv_counts[i];j++) {
466: ins_val = 0;
467: size_of_recv = recv_buffer_subset[sum_requests]; /* total size of recv from neighs */
468: for (buffer_size=0;buffer_size<size_of_recv;) { /* loop until all data from neighs has been taken into account */
469: for (k=1;k<recv_buffer[buffer_size+start_of_recv]+1;k++) { /* filling properly temp_buffer using data from a single recv */
470: temp_buffer[nodes_to_temp_buffer_indices[recv_buffer[start_of_recv+buffer_size+k]]][j] = ins_val;
471: }
472: buffer_size += k;
473: ins_val++;
474: }
475: start_of_recv += size_of_recv;
476: sum_requests++;
477: }
478: PetscMalloc1(temp_buffer_size,&add_to_subset);
479: PetscMemzero(add_to_subset,temp_buffer_size*sizeof(*add_to_subset));
480: for (j=0;j<temp_buffer_size;j++) {
481: if (!add_to_subset[j]) { /* found a new cc */
482: global_subset_counter++;
483: add_to_subset[j] = global_subset_counter;
484: for (k=j+1;k<temp_buffer_size;k++) { /* check for other nodes in new cc */
485: same_set = PETSC_TRUE;
486: for (s=0;s<cum_recv_counts[i+1]-cum_recv_counts[i];s++) {
487: if (temp_buffer[j][s]!=temp_buffer[k][s]) {
488: same_set = PETSC_FALSE;
489: break;
490: }
491: }
492: if (same_set) {
493: add_to_subset[k] = global_subset_counter;
494: }
495: }
496: }
497: }
498: /* insert new data in subset array */
499: temp_buffer_size = 0;
500: for (j=0;j<graph->nvtxs;j++) {
501: if (graph->subset[j]==i+1) {
502: graph->subset[j] = graph->n_subsets+add_to_subset[temp_buffer_size];
503: temp_buffer_size++;
504: }
505: }
506: PetscFree(temp_buffer[0]);
507: PetscFree(temp_buffer);
508: PetscFree(add_to_subset);
509: }
510: }
511: PetscFree(nodes_to_temp_buffer_indices);
512: PetscFree(sizes_of_sends);
513: PetscFree(send_requests);
514: PetscFree(recv_requests);
515: PetscFree(recv_buffer);
516: PetscFree(recv_buffer_subset);
517: PetscFree(send_buffer);
518: PetscFree(cum_recv_counts);
519: PetscFree(subset_to_nodes_indices);
520: PetscFree(subset_cc_adapt);
521: /* We are ready to find for connected components consistent among neighbouring subdomains */
522: if (global_subset_counter) {
523: PetscBTMemzero(graph->nvtxs,graph->touched);
524: global_subset_counter = 0;
525: for (i=0;i<graph->nvtxs;i++) {
526: if (graph->subset[i] && !PetscBTLookup(graph->touched,i)) {
527: global_subset_counter++;
528: for (j=i+1;j<graph->nvtxs;j++) {
529: if (!PetscBTLookup(graph->touched,j) && graph->subset[j]==graph->subset[i]) {
530: graph->subset[j] = global_subset_counter;
531: PetscBTSet(graph->touched,j);
532: }
533: }
534: graph->subset[i] = global_subset_counter;
535: PetscBTSet(graph->touched,i);
536: }
537: }
538: /* refine connected components locally */
539: PCBDDCGraphComputeConnectedComponentsLocal(graph);
540: }
541: /* restore original CSR graph of dofs */
542: PetscFree(new_xadj);
543: PetscFree(new_adjncy);
544: graph->xadj = old_xadj;
545: graph->adjncy = old_adjncy;
546: PetscFree(queue_in_global_numbering);
547: }
548: return(0);
549: }
551: /* The following code has been adapted from function IsConnectedSubdomain contained
552: in source file contig.c of METIS library (version 5.0.1)
553: It finds connected components for each subset */
556: PetscErrorCode PCBDDCGraphComputeConnectedComponentsLocal(PCBDDCGraph graph)
557: {
558: PetscInt i,j,k,first,last,nleft,ncc,pid,cum_queue,n,ncc_pid;
559: PetscInt *queue_global;
563: /* quiet return if no csr info is available */
564: if (!graph->xadj || !graph->adjncy) {
565: return(0);
566: }
568: /* reset any previous search of connected components */
569: PetscBTMemzero(graph->nvtxs,graph->touched);
570: graph->n_subsets = 0;
571: for (i=0;i<graph->nvtxs;i++) {
572: if (graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK || !graph->count[i]) {
573: PetscBTSet(graph->touched,i);
574: graph->subset[i] = 0;
575: }
576: graph->n_subsets = PetscMax(graph->n_subsets,graph->subset[i]);
577: }
578: PetscFree(graph->subset_ncc);
579: PetscMalloc1(graph->n_subsets,&graph->subset_ncc);
580: PetscMemzero(graph->subset_ncc,graph->n_subsets*sizeof(*graph->subset_ncc));
581: PetscMemzero(graph->cptr,(graph->nvtxs+1)*sizeof(*graph->cptr));
582: PetscMemzero(graph->queue,graph->nvtxs*sizeof(*graph->queue));
584: /* begin search for connected components */
585: cum_queue = 0;
586: ncc = 0;
587: for (n=0;n<graph->n_subsets;n++) {
588: pid = n+1; /* partition labeled by 0 is discarded */
589: nleft = 0;
590: for (i=0;i<graph->nvtxs;i++) {
591: if (graph->subset[i] == pid) {
592: nleft++;
593: }
594: }
595: for (i=0; i<graph->nvtxs; i++) {
596: if (graph->subset[i] == pid) {
597: break;
598: }
599: }
600: PetscBTSet(graph->touched,i);
601: graph->queue[cum_queue] = i;
602: first = 0;
603: last = 1;
604: graph->cptr[ncc] = cum_queue;
605: ncc_pid = 0;
606: while (first != nleft) {
607: if (first == last) {
608: graph->cptr[++ncc] = first+cum_queue;
609: ncc_pid++;
610: for (i=0; i<graph->nvtxs; i++) { /* TODO-> use a while! */
611: if (graph->subset[i] == pid && !PetscBTLookup(graph->touched,i)) {
612: break;
613: }
614: }
615: graph->queue[cum_queue+last] = i;
616: last++;
617: PetscBTSet(graph->touched,i);
618: }
619: i = graph->queue[cum_queue+first];
620: first++;
621: for (j=graph->xadj[i];j<graph->xadj[i+1];j++) {
622: k = graph->adjncy[j];
623: if (graph->subset[k] == pid && !PetscBTLookup(graph->touched,k)) {
624: graph->queue[cum_queue+last] = k;
625: last++;
626: PetscBTSet(graph->touched,k);
627: }
628: }
629: }
630: graph->cptr[++ncc] = first+cum_queue;
631: ncc_pid++;
632: cum_queue = graph->cptr[ncc];
633: graph->subset_ncc[n] = ncc_pid;
634: }
635: graph->ncc = ncc;
636: /* For consistency among neighbours, I need to sort (by global ordering) each connected component */
637: PetscMalloc1(graph->cptr[graph->ncc],&queue_global);
638: ISLocalToGlobalMappingApply(graph->l2gmap,graph->cptr[graph->ncc],graph->queue,queue_global);
639: for (i=0;i<graph->ncc;i++) {
640: PetscSortIntWithArray(graph->cptr[i+1]-graph->cptr[i],&queue_global[graph->cptr[i]],&graph->queue[graph->cptr[i]]);
641: }
642: PetscFree(queue_global);
643: return(0);
644: }
648: PetscErrorCode PCBDDCGraphSetUp(PCBDDCGraph graph, PetscInt custom_minimal_size, IS neumann_is, IS dirichlet_is, PetscInt n_ISForDofs, IS ISForDofs[], IS custom_primal_vertices)
649: {
650: VecScatter scatter_ctx;
651: Vec local_vec,local_vec2,global_vec;
652: IS to,from;
653: MPI_Comm comm;
654: PetscScalar *array,*array2;
655: const PetscInt *is_indices;
656: PetscInt n_neigh,*neigh,*n_shared,**shared,*queue_global;
657: PetscInt i,j,k,s,total_counts,nodes_touched,is_size;
659: PetscBool same_set,mirrors_found;
662: PetscObjectGetComm((PetscObject)(graph->l2gmap),&comm);
663: /* custom_minimal_size */
664: graph->custom_minimal_size = PetscMax(graph->custom_minimal_size,custom_minimal_size);
665: /* get info l2gmap and allocate work vectors */
666: ISLocalToGlobalMappingGetInfo(graph->l2gmap,&n_neigh,&neigh,&n_shared,&shared);
667: ISLocalToGlobalMappingGetIndices(graph->l2gmap,&is_indices);
668: j = 0;
669: for (i=0;i<graph->nvtxs;i++) {
670: j = PetscMax(j,is_indices[i]);
671: }
672: MPI_Allreduce(&j,&i,1,MPIU_INT,MPI_MAX,comm);
673: i++;
674: VecCreate(PETSC_COMM_SELF,&local_vec);
675: VecSetSizes(local_vec,PETSC_DECIDE,graph->nvtxs);
676: VecSetType(local_vec,VECSTANDARD);
677: VecDuplicate(local_vec,&local_vec2);
678: VecCreate(comm,&global_vec);
679: VecSetSizes(global_vec,PETSC_DECIDE,i);
680: VecSetType(global_vec,VECSTANDARD);
681: ISCreateStride(PETSC_COMM_SELF,graph->nvtxs,0,1,&to);
682: ISLocalToGlobalMappingApplyIS(graph->l2gmap,to,&from);
683: VecScatterCreate(global_vec,from,local_vec,to,&scatter_ctx);
685: /* get local periodic nodes */
686: mirrors_found = PETSC_FALSE;
687: if (graph->nvtxs && n_neigh) {
688: for (i=0; i<n_shared[0]; i++) graph->count[shared[0][i]] += 1;
689: for (i=0; i<n_shared[0]; i++) {
690: if (graph->count[shared[0][i]] > 1) {
691: mirrors_found = PETSC_TRUE;
692: break;
693: }
694: }
695: }
696: /* compute local mirrors (if any) */
697: if (mirrors_found) {
698: PetscInt *local_indices,*global_indices;
699: /* get arrays of local and global indices */
700: PetscMalloc1(graph->nvtxs,&local_indices);
701: ISGetIndices(to,(const PetscInt**)&is_indices);
702: PetscMemcpy(local_indices,is_indices,graph->nvtxs*sizeof(PetscInt));
703: ISRestoreIndices(to,(const PetscInt**)&is_indices);
704: PetscMalloc1(graph->nvtxs,&global_indices);
705: ISGetIndices(from,(const PetscInt**)&is_indices);
706: PetscMemcpy(global_indices,is_indices,graph->nvtxs*sizeof(PetscInt));
707: ISRestoreIndices(from,(const PetscInt**)&is_indices);
708: /* allocate space for mirrors */
709: PetscMalloc2(graph->nvtxs,&graph->mirrors,graph->nvtxs,&graph->mirrors_set);
710: PetscMemzero(graph->mirrors,graph->nvtxs*sizeof(PetscInt));
711: graph->mirrors_set[0] = 0;
713: k=0;
714: for (i=0;i<n_shared[0];i++) {
715: j=shared[0][i];
716: if (graph->count[j] > 1) {
717: graph->mirrors[j]++;
718: k++;
719: }
720: }
721: /* allocate space for set of mirrors */
722: PetscMalloc1(k,&graph->mirrors_set[0]);
723: for (i=1;i<graph->nvtxs;i++)
724: graph->mirrors_set[i]=graph->mirrors_set[i-1]+graph->mirrors[i-1];
726: /* fill arrays */
727: PetscMemzero(graph->mirrors,graph->nvtxs*sizeof(PetscInt));
728: for (j=0;j<n_shared[0];j++) {
729: i=shared[0][j];
730: if (graph->count[i] > 1)
731: graph->mirrors_set[i][graph->mirrors[i]++]=global_indices[i];
732: }
733: PetscSortIntWithArray(graph->nvtxs,global_indices,local_indices);
734: for (i=0;i<graph->nvtxs;i++) {
735: if (graph->mirrors[i] > 0) {
736: PetscFindInt(graph->mirrors_set[i][0],graph->nvtxs,global_indices,&k);
737: j = global_indices[k];
738: while ( k > 0 && global_indices[k-1] == j) k--;
739: for (j=0;j<graph->mirrors[i];j++) {
740: graph->mirrors_set[i][j]=local_indices[k+j];
741: }
742: PetscSortInt(graph->mirrors[i],graph->mirrors_set[i]);
743: }
744: }
745: PetscFree(local_indices);
746: PetscFree(global_indices);
747: }
748: PetscMemzero(graph->count,graph->nvtxs*sizeof(*graph->count));
749: ISDestroy(&to);
750: ISDestroy(&from);
752: /* Count total number of neigh per node */
753: k=0;
754: for (i=1;i<n_neigh;i++) {
755: k += n_shared[i];
756: for (j=0;j<n_shared[i];j++) {
757: graph->count[shared[i][j]] += 1;
758: }
759: }
760: /* Allocate space for storing the set of neighbours for each node */
761: if (graph->nvtxs) {
762: PetscMalloc1(k,&graph->neighbours_set[0]);
763: }
764: for (i=1;i<graph->nvtxs;i++) { /* dont count myself */
765: graph->neighbours_set[i]=graph->neighbours_set[i-1]+graph->count[i-1];
766: }
767: /* Get information for sharing subdomains */
768: PetscMemzero(graph->count,graph->nvtxs*sizeof(*graph->count));
769: for (i=1;i<n_neigh;i++) { /* dont count myself */
770: s = n_shared[i];
771: for (j=0;j<s;j++) {
772: k = shared[i][j];
773: graph->neighbours_set[k][graph->count[k]] = neigh[i];
774: graph->count[k] += 1;
775: }
776: }
777: /* sort set of sharing subdomains */
778: for (i=0;i<graph->nvtxs;i++) {
779: PetscSortRemoveDupsInt(&graph->count[i],graph->neighbours_set[i]);
780: }
781: /*
782: Get info for dofs splitting
783: User can specify only a subset; an additional field is considered as a complementary set
784: */
785: for (i=0;i<graph->nvtxs;i++) {
786: graph->which_dof[i] = n_ISForDofs; /* by default a dof belongs to the complement set */
787: }
788: for (i=0;i<n_ISForDofs;i++) {
789: ISGetLocalSize(ISForDofs[i],&is_size);
790: ISGetIndices(ISForDofs[i],(const PetscInt**)&is_indices);
791: for (j=0;j<is_size;j++) {
792: if (is_indices[j] > -1 && is_indices[j] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */
793: graph->which_dof[is_indices[j]] = i;
794: }
795: }
796: ISRestoreIndices(ISForDofs[i],(const PetscInt**)&is_indices);
797: }
798: /* Take into account Neumann nodes */
799: VecSet(local_vec,0.0);
800: VecSet(local_vec2,0.0);
801: if (neumann_is) {
802: VecGetArray(local_vec,&array);
803: ISGetLocalSize(neumann_is,&is_size);
804: ISGetIndices(neumann_is,(const PetscInt**)&is_indices);
805: for (i=0;i<is_size;i++) {
806: if (is_indices[i] > -1 && is_indices[i] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */
807: array[is_indices[i]] = 1.0;
808: }
809: }
810: ISRestoreIndices(neumann_is,(const PetscInt**)&is_indices);
811: VecRestoreArray(local_vec,&array);
812: }
813: /* Neumann nodes: impose consistency among neighbours */
814: VecSet(global_vec,0.0);
815: VecScatterBegin(scatter_ctx,local_vec,global_vec,ADD_VALUES,SCATTER_REVERSE);
816: VecScatterEnd(scatter_ctx,local_vec,global_vec,ADD_VALUES,SCATTER_REVERSE);
817: VecScatterBegin(scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_FORWARD);
818: VecScatterEnd(scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_FORWARD);
819: VecGetArray(local_vec,&array);
820: for (i=0;i<graph->nvtxs;i++) {
821: if (PetscRealPart(array[i]) > 0.1) {
822: graph->special_dof[i] = PCBDDCGRAPH_NEUMANN_MARK;
823: }
824: }
825: VecRestoreArray(local_vec,&array);
826: /* Take into account Dirichlet nodes */
827: VecSet(local_vec2,0.0);
828: if (dirichlet_is) {
829: VecGetArray(local_vec,&array);
830: VecGetArray(local_vec2,&array2);
831: ISGetLocalSize(dirichlet_is,&is_size);
832: ISGetIndices(dirichlet_is,(const PetscInt**)&is_indices);
833: for (i=0;i<is_size;i++){
834: if (is_indices[i] > -1 && is_indices[i] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */
835: k = is_indices[i];
836: if (graph->count[k] && !PetscBTLookup(graph->touched,k)) {
837: if (PetscRealPart(array[k]) > 0.1) {
838: 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);
839: }
840: array2[k] = 1.0;
841: }
842: }
843: }
844: ISRestoreIndices(dirichlet_is,(const PetscInt**)&is_indices);
845: VecRestoreArray(local_vec,&array);
846: VecRestoreArray(local_vec2,&array2);
847: }
848: /* Dirichlet nodes: impose consistency among neighbours */
849: VecSet(global_vec,0.0);
850: VecScatterBegin(scatter_ctx,local_vec2,global_vec,ADD_VALUES,SCATTER_REVERSE);
851: VecScatterEnd(scatter_ctx,local_vec2,global_vec,ADD_VALUES,SCATTER_REVERSE);
852: VecScatterBegin(scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_FORWARD);
853: VecScatterEnd(scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_FORWARD);
854: VecGetArray(local_vec,&array);
855: for (i=0;i<graph->nvtxs;i++) {
856: if (PetscRealPart(array[i]) > 0.1) {
857: PetscBTSet(graph->touched,i);
858: graph->subset[i] = 0; /* dirichlet nodes treated as internal -> is it ok? */
859: graph->special_dof[i] = PCBDDCGRAPH_DIRICHLET_MARK;
860: }
861: }
862: VecRestoreArray(local_vec,&array);
864: /* mark local periodic nodes (if any) and adapt CSR graph (if any) */
865: if (graph->mirrors) {
866: for (i=0;i<graph->nvtxs;i++)
867: if (graph->mirrors[i])
868: graph->special_dof[i] = PCBDDCGRAPH_LOCAL_PERIODIC_MARK;
870: if (graph->xadj && graph->adjncy) {
871: PetscInt *new_xadj,*new_adjncy;
872: /* sort CSR graph */
873: for (i=0;i<graph->nvtxs;i++)
874: PetscSortInt(graph->xadj[i+1]-graph->xadj[i],&graph->adjncy[graph->xadj[i]]);
876: /* adapt local CSR graph in case of local periodicity */
877: k=0;
878: for (i=0;i<graph->nvtxs;i++)
879: for (j=graph->xadj[i];j<graph->xadj[i+1];j++)
880: k += graph->mirrors[graph->adjncy[j]];
882: PetscMalloc1((graph->nvtxs+1),&new_xadj);
883: PetscMalloc1((k+graph->xadj[graph->nvtxs]),&new_adjncy);
884: new_xadj[0]=0;
885: for (i=0;i<graph->nvtxs;i++) {
886: k = graph->xadj[i+1]-graph->xadj[i];
887: PetscMemcpy(&new_adjncy[new_xadj[i]],&graph->adjncy[graph->xadj[i]],k*sizeof(PetscInt));
888: new_xadj[i+1]=new_xadj[i]+k;
889: for (j=graph->xadj[i];j<graph->xadj[i+1];j++) {
890: k = graph->mirrors[graph->adjncy[j]];
891: PetscMemcpy(&new_adjncy[new_xadj[i+1]],graph->mirrors_set[graph->adjncy[j]],k*sizeof(PetscInt));
892: new_xadj[i+1]+=k;
893: }
894: k = new_xadj[i+1]-new_xadj[i];
895: PetscSortRemoveDupsInt(&k,&new_adjncy[new_xadj[i]]);
896: new_xadj[i+1]=new_xadj[i]+k;
897: }
898: /* set new CSR into graph */
899: PetscFree(graph->xadj);
900: PetscFree(graph->adjncy);
901: graph->xadj = new_xadj;
902: graph->adjncy = new_adjncy;
903: }
904: }
906: /* mark special nodes -> each will become a single node equivalence class */
907: if (custom_primal_vertices) {
908: ISGetLocalSize(custom_primal_vertices,&is_size);
909: ISGetIndices(custom_primal_vertices,(const PetscInt**)&is_indices);
910: for (i=0;i<is_size;i++) {
911: graph->special_dof[is_indices[i]] = PCBDDCGRAPH_SPECIAL_MARK-i;
912: }
913: ISRestoreIndices(custom_primal_vertices,(const PetscInt**)&is_indices);
914: }
915: /* mark interior nodes as touched and belonging to partition number 0 */
916: for (i=0;i<graph->nvtxs;i++) {
917: if (!graph->count[i]) {
918: PetscBTSet(graph->touched,i);
919: graph->subset[i] = 0;
920: }
921: }
922: /* init graph structure and compute default subsets */
923: nodes_touched=0;
924: for (i=0;i<graph->nvtxs;i++) {
925: if (PetscBTLookup(graph->touched,i)) {
926: nodes_touched++;
927: }
928: }
929: i = 0;
930: graph->ncc = 0;
931: total_counts = 0;
932: while (nodes_touched<graph->nvtxs) {
933: /* find first untouched node in local ordering */
934: while (PetscBTLookup(graph->touched,i)) {
935: i++;
936: }
937: PetscBTSet(graph->touched,i);
938: graph->subset[i] = graph->ncc+1;
939: graph->cptr[graph->ncc] = total_counts;
940: graph->queue[total_counts] = i;
941: total_counts++;
942: nodes_touched++;
943: /* now find all other nodes having the same set of sharing subdomains */
944: for (j=i+1;j<graph->nvtxs;j++) {
945: /* check for same number of sharing subdomains, dof number and same special mark */
946: 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]) {
947: /* check for same set of sharing subdomains */
948: same_set=PETSC_TRUE;
949: for (k=0;k<graph->count[j];k++){
950: if (graph->neighbours_set[i][k]!=graph->neighbours_set[j][k]) {
951: same_set=PETSC_FALSE;
952: }
953: }
954: /* I found a friend of mine */
955: if (same_set) {
956: graph->subset[j]=graph->ncc+1;
957: PetscBTSet(graph->touched,j);
958: nodes_touched++;
959: graph->queue[total_counts] = j;
960: total_counts++;
961: }
962: }
963: }
964: graph->ncc++;
965: }
966: /* set default number of subsets (at this point no info on csr graph has been taken into account, so n_subsets = ncc */
967: graph->n_subsets = graph->ncc;
968: PetscMalloc1(graph->n_subsets,&graph->subset_ncc);
969: for (i=0;i<graph->n_subsets;i++) {
970: graph->subset_ncc[i] = 1;
971: }
972: /* final pointer */
973: graph->cptr[graph->ncc] = total_counts;
974: /* free memory allocated by ISLocalToGlobalMappingGetInfo */
975: ISLocalToGlobalMappingRestoreInfo(graph->l2gmap,&n_neigh,&neigh,&n_shared,&shared);
976: /* get a reference node (min index in global ordering) for each subset */
977: PetscMalloc1(graph->ncc,&graph->subset_ref_node);
978: PetscMalloc1(graph->cptr[graph->ncc],&queue_global);
979: ISLocalToGlobalMappingApply(graph->l2gmap,graph->cptr[graph->ncc],graph->queue,queue_global);
980: for (i=0;i<graph->ncc;i++) {
981: PetscInt minval = queue_global[graph->cptr[i]];
982: for (j=graph->cptr[i]+1;j<graph->cptr[i+1];j++) {
983: minval = PetscMin(minval,queue_global[j]);
984: }
985: graph->subset_ref_node[i] = minval;
986: }
987: PetscFree(queue_global);
988: /* free objects */
989: VecDestroy(&local_vec);
990: VecDestroy(&local_vec2);
991: VecDestroy(&global_vec);
992: VecScatterDestroy(&scatter_ctx);
993: return(0);
994: }
998: PetscErrorCode PCBDDCGraphResetCSR(PCBDDCGraph graph)
999: {
1003: PetscFree(graph->xadj);
1004: PetscFree(graph->adjncy);
1005: graph->nvtxs_csr = 0;
1006: return(0);
1007: }
1011: PetscErrorCode PCBDDCGraphReset(PCBDDCGraph graph)
1012: {
1016: ISLocalToGlobalMappingDestroy(&graph->l2gmap);
1017: PetscFree(graph->subset_ncc);
1018: PetscFree(graph->subset_ref_node);
1019: if (graph->nvtxs) {
1020: PetscFree(graph->neighbours_set[0]);
1021: }
1022: PetscBTDestroy(&graph->touched);
1023: PetscFree7(graph->count,
1024: graph->neighbours_set,
1025: graph->subset,
1026: graph->which_dof,
1027: graph->cptr,
1028: graph->queue,
1029: graph->special_dof);
1030: if (graph->mirrors) {
1031: PetscFree(graph->mirrors_set[0]);
1032: }
1033: PetscFree2(graph->mirrors,graph->mirrors_set);
1034: graph->nvtxs = 0;
1035: graph->n_subsets = 0;
1036: graph->custom_minimal_size = 1;
1037: return(0);
1038: }
1042: PetscErrorCode PCBDDCGraphInit(PCBDDCGraph graph, ISLocalToGlobalMapping l2gmap)
1043: {
1044: PetscInt n;
1050: /* raise an error if already allocated */
1051: if (graph->nvtxs) {
1052: SETERRQ(PetscObjectComm((PetscObject)l2gmap),PETSC_ERR_PLIB,"BDDCGraph already initialized");
1053: }
1054: /* set number of vertices */
1055: PetscObjectReference((PetscObject)l2gmap);
1056: graph->l2gmap = l2gmap;
1057: ISLocalToGlobalMappingGetSize(l2gmap,&n);
1058: graph->nvtxs = n;
1059: /* allocate used space */
1060: PetscBTCreate(graph->nvtxs,&graph->touched);
1061: PetscMalloc7(graph->nvtxs,&graph->count,
1062: graph->nvtxs,&graph->neighbours_set,
1063: graph->nvtxs,&graph->subset,
1064: graph->nvtxs,&graph->which_dof,
1065: graph->nvtxs+1,&graph->cptr,
1066: graph->nvtxs,&graph->queue,
1067: graph->nvtxs,&graph->special_dof);
1068: /* zeroes memory */
1069: PetscMemzero(graph->count,graph->nvtxs*sizeof(PetscInt));
1070: PetscMemzero(graph->subset,graph->nvtxs*sizeof(PetscInt));
1071: /* use -1 as a default value for which_dof array */
1072: for (n=0;n<graph->nvtxs;n++) graph->which_dof[n] = -1;
1073: PetscMemzero(graph->cptr,(graph->nvtxs+1)*sizeof(PetscInt));
1074: PetscMemzero(graph->queue,graph->nvtxs*sizeof(PetscInt));
1075: PetscMemzero(graph->special_dof,graph->nvtxs*sizeof(PetscInt));
1076: /* zeroes first pointer to neighbour set */
1077: if (graph->nvtxs) {
1078: graph->neighbours_set[0] = 0;
1079: }
1080: /* zeroes workspace for values of ncc */
1081: graph->subset_ncc = 0;
1082: graph->subset_ref_node = 0;
1083: return(0);
1084: }
1088: PetscErrorCode PCBDDCGraphDestroy(PCBDDCGraph* graph)
1089: {
1093: PCBDDCGraphReset(*graph);
1094: PetscFree(*graph);
1095: return(0);
1096: }
1100: PetscErrorCode PCBDDCGraphCreate(PCBDDCGraph *graph)
1101: {
1102: PCBDDCGraph new_graph;
1106: PetscMalloc(sizeof(*new_graph),&new_graph);
1107: /* local to global mapping of dofs */
1108: new_graph->l2gmap = 0;
1109: /* vertex size */
1110: new_graph->nvtxs = 0;
1111: new_graph->n_subsets = 0;
1112: new_graph->custom_minimal_size = 1;
1113: /* zeroes ponters */
1114: new_graph->mirrors = 0;
1115: new_graph->mirrors_set = 0;
1116: new_graph->neighbours_set = 0;
1117: new_graph->subset = 0;
1118: new_graph->which_dof = 0;
1119: new_graph->special_dof = 0;
1120: new_graph->cptr = 0;
1121: new_graph->queue = 0;
1122: new_graph->count = 0;
1123: new_graph->subset_ncc = 0;
1124: new_graph->subset_ref_node = 0;
1125: new_graph->touched = 0;
1126: /* zeroes pointers to csr graph of local nodes connectivity (optional data) */
1127: new_graph->nvtxs_csr = 0;
1128: new_graph->xadj = 0;
1129: new_graph->adjncy = 0;
1130: *graph = new_graph;
1131: return(0);
1132: }