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