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