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