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