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