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