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