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