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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 206 } 207 PetscCall(PetscCalloc1(graph->ncc, &mark)); 208 /* loop on ccs to evaluate 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(PETSC_SUCCESS); 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 PetscCheck(buffer_size >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected buffer_size %" PetscInt_FMT " >= 0", buffer_size); 599 /* compute pointers */ 600 for (j = 1; j < buffer_size; j++) refine_buffer[j] = refine_buffer[j - 1] + sharingprocs; 601 /* analyze contributions from subdomains that share the i-th subset 602 The structure of refine_buffer is suitable to find intersections of ccs among sharingprocs. 603 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) 604 sharing procs connected components: 605 neigh 0: [0 1 4], [2 3], labels [4,7] (2 connected components) 606 neigh 1: [0 1], [2 3 4], labels [3 2] (2 connected components) 607 neigh 2: [0 4], [1], [2 3], labels [1 5 6] (3 connected components) 608 refine_buffer will be filled as: 609 [ 4, 3, 1; 610 4, 2, 1; 611 7, 2, 6; 612 4, 3, 5; 613 7, 2, 6; ]; 614 The connected components in local ordering are [0], [1], [2 3], [4] */ 615 /* fill temp_buffer */ 616 for (k = 0; k < buffer_size; k++) refine_buffer[k][0] = labels[graph->subset_idxs[i][k]]; 617 for (j = 0; j < sharingprocs - 1; j++) { 618 for (k = 0; k < buffer_size; k++) refine_buffer[k][j + 1] = recv_buffer[start_of_recv + k]; 619 start_of_recv += buffer_size; 620 } 621 PetscCall(PetscArrayzero(private_labels, buffer_size)); 622 for (j = 0; j < buffer_size; j++) { 623 if (!private_labels[j]) { /* found a new cc */ 624 PetscBool same_set; 625 626 graph->cptr[ncc] = cum_queue; 627 ncc++; 628 subset_counter++; 629 private_labels[j] = subset_counter; 630 graph->queue[cum_queue++] = graph->subset_idxs[i][j]; 631 for (k = j + 1; k < buffer_size; k++) { /* check for other nodes in new cc */ 632 same_set = PETSC_TRUE; 633 for (s = 0; s < sharingprocs; s++) { 634 if (refine_buffer[j][s] != refine_buffer[k][s]) { 635 same_set = PETSC_FALSE; 636 break; 637 } 638 } 639 if (same_set) { 640 private_labels[k] = subset_counter; 641 graph->queue[cum_queue++] = graph->subset_idxs[i][k]; 642 } 643 } 644 } 645 } 646 graph->cptr[ncc] = cum_queue; 647 graph->subset_ncc[i] = subset_counter; 648 graph->queue_sorted = PETSC_FALSE; 649 } else { /* this subset does not need to be adapted */ 650 PetscCall(PetscArraycpy(graph->queue + cum_queue, graph->subset_idxs[i], graph->subset_size[i])); 651 ncc++; 652 cum_queue += graph->subset_size[i]; 653 graph->cptr[ncc] = cum_queue; 654 } 655 } 656 graph->cptr[ncc] = cum_queue; 657 graph->ncc = ncc; 658 if (mss) { 659 PetscCall(PetscFree2(refine_buffer[0], private_labels)); 660 PetscCall(PetscFree(refine_buffer)); 661 } 662 PetscCall(PetscFree(labels)); 663 PetscCallMPI(MPI_Waitall(sum_requests, send_requests, MPI_STATUSES_IGNORE)); 664 PetscCall(PetscFree2(send_requests, recv_requests)); 665 PetscCall(PetscFree2(send_buffer, recv_buffer)); 666 PetscCall(PetscFree(cum_recv_counts)); 667 PetscCall(PetscBTDestroy(&subset_cc_adapt)); 668 } 669 PetscCall(PetscBTDestroy(&cornerp)); 670 671 /* Determine if we are in 2D or 3D */ 672 if (!graph->twodimset) { 673 PetscBool twodim = PETSC_TRUE; 674 for (i = 0; i < graph->ncc; i++) { 675 PetscInt repdof = graph->queue[graph->cptr[i]]; 676 PetscInt ccsize = graph->cptr[i + 1] - graph->cptr[i]; 677 if (graph->count[repdof] > 1 && ccsize > graph->custom_minimal_size) { 678 twodim = PETSC_FALSE; 679 break; 680 } 681 } 682 PetscCall(MPIU_Allreduce(&twodim, &graph->twodim, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)graph->l2gmap))); 683 graph->twodimset = PETSC_TRUE; 684 } 685 PetscFunctionReturn(PETSC_SUCCESS); 686 } 687 688 static inline PetscErrorCode PCBDDCGraphComputeCC_Private(PCBDDCGraph graph, PetscInt pid, PetscInt *queue_tip, PetscInt n_prev, PetscInt *n_added) 689 { 690 PetscInt i, j, n; 691 PetscInt *xadj = graph->xadj, *adjncy = graph->adjncy; 692 PetscBT touched = graph->touched; 693 PetscBool havecsr = (PetscBool)(!!xadj); 694 PetscBool havesubs = (PetscBool)(!!graph->n_local_subs); 695 696 PetscFunctionBegin; 697 n = 0; 698 if (havecsr && !havesubs) { 699 for (i = -n_prev; i < 0; i++) { 700 PetscInt start_dof = queue_tip[i]; 701 /* we assume that if a dof has a size 1 adjacency list and the corresponding entry is negative, it is connected to all dofs */ 702 if (xadj[start_dof + 1] - xadj[start_dof] == 1 && adjncy[xadj[start_dof]] < 0) { 703 for (j = 0; j < graph->subset_size[pid - 1]; j++) { /* pid \in [1,graph->n_subsets] */ 704 PetscInt dof = graph->subset_idxs[pid - 1][j]; 705 if (!PetscBTLookup(touched, dof) && graph->subset[dof] == pid) { 706 PetscCall(PetscBTSet(touched, dof)); 707 queue_tip[n] = dof; 708 n++; 709 } 710 } 711 } else { 712 for (j = xadj[start_dof]; j < xadj[start_dof + 1]; j++) { 713 PetscInt dof = adjncy[j]; 714 if (!PetscBTLookup(touched, dof) && graph->subset[dof] == pid) { 715 PetscCall(PetscBTSet(touched, dof)); 716 queue_tip[n] = dof; 717 n++; 718 } 719 } 720 } 721 } 722 } else if (havecsr && havesubs) { 723 PetscInt sid = graph->local_subs[queue_tip[-n_prev]]; 724 for (i = -n_prev; i < 0; i++) { 725 PetscInt start_dof = queue_tip[i]; 726 /* 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 */ 727 if (xadj[start_dof + 1] - xadj[start_dof] == 1 && adjncy[xadj[start_dof]] < 0) { 728 for (j = 0; j < graph->subset_size[pid - 1]; j++) { /* pid \in [1,graph->n_subsets] */ 729 PetscInt dof = graph->subset_idxs[pid - 1][j]; 730 if (!PetscBTLookup(touched, dof) && graph->subset[dof] == pid && graph->local_subs[dof] == sid) { 731 PetscCall(PetscBTSet(touched, dof)); 732 queue_tip[n] = dof; 733 n++; 734 } 735 } 736 } else { 737 for (j = xadj[start_dof]; j < xadj[start_dof + 1]; j++) { 738 PetscInt dof = adjncy[j]; 739 if (!PetscBTLookup(touched, dof) && graph->subset[dof] == pid && graph->local_subs[dof] == sid) { 740 PetscCall(PetscBTSet(touched, dof)); 741 queue_tip[n] = dof; 742 n++; 743 } 744 } 745 } 746 } 747 } else if (havesubs) { /* sub info only */ 748 PetscInt sid = graph->local_subs[queue_tip[-n_prev]]; 749 for (j = 0; j < graph->subset_size[pid - 1]; j++) { /* pid \in [1,graph->n_subsets] */ 750 PetscInt dof = graph->subset_idxs[pid - 1][j]; 751 if (!PetscBTLookup(touched, dof) && graph->subset[dof] == pid && graph->local_subs[dof] == sid) { 752 PetscCall(PetscBTSet(touched, dof)); 753 queue_tip[n] = dof; 754 n++; 755 } 756 } 757 } else { 758 for (j = 0; j < graph->subset_size[pid - 1]; j++) { /* pid \in [1,graph->n_subsets] */ 759 PetscInt dof = graph->subset_idxs[pid - 1][j]; 760 if (!PetscBTLookup(touched, dof) && graph->subset[dof] == pid) { 761 PetscCall(PetscBTSet(touched, dof)); 762 queue_tip[n] = dof; 763 n++; 764 } 765 } 766 } 767 *n_added = n; 768 PetscFunctionReturn(PETSC_SUCCESS); 769 } 770 771 PetscErrorCode PCBDDCGraphComputeConnectedComponentsLocal(PCBDDCGraph graph) 772 { 773 PetscInt ncc, cum_queue, n; 774 PetscMPIInt commsize; 775 776 PetscFunctionBegin; 777 PetscCheck(graph->setupcalled, PetscObjectComm((PetscObject)graph->l2gmap), PETSC_ERR_ORDER, "PCBDDCGraphSetUp should be called first"); 778 /* quiet return if there isn't any local info */ 779 if (!graph->xadj && !graph->n_local_subs) PetscFunctionReturn(PETSC_SUCCESS); 780 781 /* reset any previous search of connected components */ 782 PetscCall(PetscBTMemzero(graph->nvtxs, graph->touched)); 783 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)graph->l2gmap), &commsize)); 784 if (commsize > graph->commsizelimit) { 785 PetscInt i; 786 for (i = 0; i < graph->nvtxs; i++) { 787 if (graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK || !graph->count[i]) PetscCall(PetscBTSet(graph->touched, i)); 788 } 789 } 790 791 /* begin search for connected components */ 792 cum_queue = 0; 793 ncc = 0; 794 for (n = 0; n < graph->n_subsets; n++) { 795 PetscInt pid = n + 1; /* partition labeled by 0 is discarded */ 796 PetscInt found = 0, prev = 0, first = 0, ncc_pid = 0; 797 while (found != graph->subset_size[n]) { 798 PetscInt added = 0; 799 if (!prev) { /* search for new starting dof */ 800 while (PetscBTLookup(graph->touched, graph->subset_idxs[n][first])) first++; 801 PetscCall(PetscBTSet(graph->touched, graph->subset_idxs[n][first])); 802 graph->queue[cum_queue] = graph->subset_idxs[n][first]; 803 graph->cptr[ncc] = cum_queue; 804 prev = 1; 805 cum_queue++; 806 found++; 807 ncc_pid++; 808 ncc++; 809 } 810 PetscCall(PCBDDCGraphComputeCC_Private(graph, pid, graph->queue + cum_queue, prev, &added)); 811 if (!added) { 812 graph->subset_ncc[n] = ncc_pid; 813 graph->cptr[ncc] = cum_queue; 814 } 815 prev = added; 816 found += added; 817 cum_queue += added; 818 if (added && found == graph->subset_size[n]) { 819 graph->subset_ncc[n] = ncc_pid; 820 graph->cptr[ncc] = cum_queue; 821 } 822 } 823 } 824 graph->ncc = ncc; 825 graph->queue_sorted = PETSC_FALSE; 826 PetscFunctionReturn(PETSC_SUCCESS); 827 } 828 829 PetscErrorCode PCBDDCGraphSetUp(PCBDDCGraph graph, PetscInt custom_minimal_size, IS neumann_is, IS dirichlet_is, PetscInt n_ISForDofs, IS ISForDofs[], IS custom_primal_vertices) 830 { 831 IS subset, subset_n; 832 MPI_Comm comm; 833 const PetscInt *is_indices; 834 PetscInt n_neigh, *neigh, *n_shared, **shared, *queue_global; 835 PetscInt i, j, k, s, total_counts, nodes_touched, is_size; 836 PetscMPIInt commsize; 837 PetscBool same_set, mirrors_found; 838 839 PetscFunctionBegin; 840 PetscValidLogicalCollectiveInt(graph->l2gmap, custom_minimal_size, 2); 841 if (neumann_is) { 842 PetscValidHeaderSpecific(neumann_is, IS_CLASSID, 3); 843 PetscCheckSameComm(graph->l2gmap, 1, neumann_is, 3); 844 } 845 graph->has_dirichlet = PETSC_FALSE; 846 if (dirichlet_is) { 847 PetscValidHeaderSpecific(dirichlet_is, IS_CLASSID, 4); 848 PetscCheckSameComm(graph->l2gmap, 1, dirichlet_is, 4); 849 graph->has_dirichlet = PETSC_TRUE; 850 } 851 PetscValidLogicalCollectiveInt(graph->l2gmap, n_ISForDofs, 5); 852 for (i = 0; i < n_ISForDofs; i++) { 853 PetscValidHeaderSpecific(ISForDofs[i], IS_CLASSID, 6); 854 PetscCheckSameComm(graph->l2gmap, 1, ISForDofs[i], 6); 855 } 856 if (custom_primal_vertices) { 857 PetscValidHeaderSpecific(custom_primal_vertices, IS_CLASSID, 7); 858 PetscCheckSameComm(graph->l2gmap, 1, custom_primal_vertices, 7); 859 } 860 PetscCall(PetscObjectGetComm((PetscObject)graph->l2gmap, &comm)); 861 PetscCallMPI(MPI_Comm_size(comm, &commsize)); 862 863 /* custom_minimal_size */ 864 graph->custom_minimal_size = custom_minimal_size; 865 /* get info l2gmap and allocate work vectors */ 866 PetscCall(ISLocalToGlobalMappingGetInfo(graph->l2gmap, &n_neigh, &neigh, &n_shared, &shared)); 867 /* check if we have any local periodic nodes (periodic BCs) */ 868 mirrors_found = PETSC_FALSE; 869 if (graph->nvtxs && n_neigh) { 870 for (i = 0; i < n_shared[0]; i++) graph->count[shared[0][i]] += 1; 871 for (i = 0; i < n_shared[0]; i++) { 872 if (graph->count[shared[0][i]] > 1) { 873 mirrors_found = PETSC_TRUE; 874 break; 875 } 876 } 877 } 878 /* compute local mirrors (if any) */ 879 if (mirrors_found) { 880 IS to, from; 881 PetscInt *local_indices, *global_indices; 882 883 PetscCall(ISCreateStride(PETSC_COMM_SELF, graph->nvtxs, 0, 1, &to)); 884 PetscCall(ISLocalToGlobalMappingApplyIS(graph->l2gmap, to, &from)); 885 /* get arrays of local and global indices */ 886 PetscCall(PetscMalloc1(graph->nvtxs, &local_indices)); 887 PetscCall(ISGetIndices(to, (const PetscInt **)&is_indices)); 888 PetscCall(PetscArraycpy(local_indices, is_indices, graph->nvtxs)); 889 PetscCall(ISRestoreIndices(to, (const PetscInt **)&is_indices)); 890 PetscCall(PetscMalloc1(graph->nvtxs, &global_indices)); 891 PetscCall(ISGetIndices(from, (const PetscInt **)&is_indices)); 892 PetscCall(PetscArraycpy(global_indices, is_indices, graph->nvtxs)); 893 PetscCall(ISRestoreIndices(from, (const PetscInt **)&is_indices)); 894 /* allocate space for mirrors */ 895 PetscCall(PetscMalloc2(graph->nvtxs, &graph->mirrors, graph->nvtxs, &graph->mirrors_set)); 896 PetscCall(PetscArrayzero(graph->mirrors, graph->nvtxs)); 897 graph->mirrors_set[0] = NULL; 898 899 k = 0; 900 for (i = 0; i < n_shared[0]; i++) { 901 j = shared[0][i]; 902 if (graph->count[j] > 1) { 903 graph->mirrors[j]++; 904 k++; 905 } 906 } 907 /* allocate space for set of mirrors */ 908 PetscCall(PetscMalloc1(k, &graph->mirrors_set[0])); 909 for (i = 1; i < graph->nvtxs; i++) graph->mirrors_set[i] = graph->mirrors_set[i - 1] + graph->mirrors[i - 1]; 910 911 /* fill arrays */ 912 PetscCall(PetscArrayzero(graph->mirrors, graph->nvtxs)); 913 for (j = 0; j < n_shared[0]; j++) { 914 i = shared[0][j]; 915 if (graph->count[i] > 1) graph->mirrors_set[i][graph->mirrors[i]++] = global_indices[i]; 916 } 917 PetscCall(PetscSortIntWithArray(graph->nvtxs, global_indices, local_indices)); 918 for (i = 0; i < graph->nvtxs; i++) { 919 if (graph->mirrors[i] > 0) { 920 PetscCall(PetscFindInt(graph->mirrors_set[i][0], graph->nvtxs, global_indices, &k)); 921 j = global_indices[k]; 922 while (k > 0 && global_indices[k - 1] == j) k--; 923 for (j = 0; j < graph->mirrors[i]; j++) graph->mirrors_set[i][j] = local_indices[k + j]; 924 PetscCall(PetscSortInt(graph->mirrors[i], graph->mirrors_set[i])); 925 } 926 } 927 PetscCall(PetscFree(local_indices)); 928 PetscCall(PetscFree(global_indices)); 929 PetscCall(ISDestroy(&to)); 930 PetscCall(ISDestroy(&from)); 931 } 932 PetscCall(PetscArrayzero(graph->count, graph->nvtxs)); 933 934 /* Count total number of neigh per node */ 935 k = 0; 936 for (i = 1; i < n_neigh; i++) { 937 k += n_shared[i]; 938 for (j = 0; j < n_shared[i]; j++) graph->count[shared[i][j]] += 1; 939 } 940 /* Allocate space for storing the set of neighbours for each node */ 941 if (graph->nvtxs) PetscCall(PetscMalloc1(k, &graph->neighbours_set[0])); 942 for (i = 1; i < graph->nvtxs; i++) { /* dont count myself */ 943 graph->neighbours_set[i] = PetscSafePointerPlusOffset(graph->neighbours_set[i - 1], graph->count[i - 1]); 944 } 945 /* Get information for sharing subdomains */ 946 PetscCall(PetscArrayzero(graph->count, graph->nvtxs)); 947 for (i = 1; i < n_neigh; i++) { /* dont count myself */ 948 s = n_shared[i]; 949 for (j = 0; j < s; j++) { 950 k = shared[i][j]; 951 graph->neighbours_set[k][graph->count[k]] = neigh[i]; 952 graph->count[k] += 1; 953 } 954 } 955 /* sort set of sharing subdomains */ 956 for (i = 0; i < graph->nvtxs; i++) PetscCall(PetscSortRemoveDupsInt(&graph->count[i], graph->neighbours_set[i])); 957 /* free memory allocated by ISLocalToGlobalMappingGetInfo */ 958 PetscCall(ISLocalToGlobalMappingRestoreInfo(graph->l2gmap, &n_neigh, &neigh, &n_shared, &shared)); 959 960 /* 961 Get info for dofs splitting 962 User can specify just a subset; an additional field is considered as a complementary field 963 */ 964 for (i = 0, k = 0; i < n_ISForDofs; i++) { 965 PetscInt bs; 966 967 PetscCall(ISGetBlockSize(ISForDofs[i], &bs)); 968 k += bs; 969 } 970 for (i = 0; i < graph->nvtxs; i++) graph->which_dof[i] = k; /* by default a dof belongs to the complement set */ 971 for (i = 0, k = 0; i < n_ISForDofs; i++) { 972 PetscInt bs; 973 974 PetscCall(ISGetLocalSize(ISForDofs[i], &is_size)); 975 PetscCall(ISGetBlockSize(ISForDofs[i], &bs)); 976 PetscCall(ISGetIndices(ISForDofs[i], (const PetscInt **)&is_indices)); 977 for (j = 0; j < is_size / bs; j++) { 978 PetscInt b; 979 980 for (b = 0; b < bs; b++) { 981 PetscInt jj = bs * j + b; 982 983 if (is_indices[jj] > -1 && is_indices[jj] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */ 984 graph->which_dof[is_indices[jj]] = k + b; 985 } 986 } 987 } 988 PetscCall(ISRestoreIndices(ISForDofs[i], (const PetscInt **)&is_indices)); 989 k += bs; 990 } 991 992 /* Take into account Neumann nodes */ 993 if (neumann_is) { 994 PetscCall(ISGetLocalSize(neumann_is, &is_size)); 995 PetscCall(ISGetIndices(neumann_is, (const PetscInt **)&is_indices)); 996 for (i = 0; i < is_size; i++) { 997 if (is_indices[i] > -1 && is_indices[i] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */ 998 graph->special_dof[is_indices[i]] = PCBDDCGRAPH_NEUMANN_MARK; 999 } 1000 } 1001 PetscCall(ISRestoreIndices(neumann_is, (const PetscInt **)&is_indices)); 1002 } 1003 /* Take into account Dirichlet nodes (they overwrite any neumann boundary mark previously set) */ 1004 if (dirichlet_is) { 1005 PetscCall(ISGetLocalSize(dirichlet_is, &is_size)); 1006 PetscCall(ISGetIndices(dirichlet_is, (const PetscInt **)&is_indices)); 1007 for (i = 0; i < is_size; i++) { 1008 if (is_indices[i] > -1 && is_indices[i] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */ 1009 if (commsize > graph->commsizelimit) { /* dirichlet nodes treated as internal */ 1010 PetscCall(PetscBTSet(graph->touched, is_indices[i])); 1011 graph->subset[is_indices[i]] = 0; 1012 } 1013 graph->special_dof[is_indices[i]] = PCBDDCGRAPH_DIRICHLET_MARK; 1014 } 1015 } 1016 PetscCall(ISRestoreIndices(dirichlet_is, (const PetscInt **)&is_indices)); 1017 } 1018 /* mark local periodic nodes (if any) and adapt CSR graph (if any) */ 1019 if (graph->mirrors) { 1020 for (i = 0; i < graph->nvtxs; i++) 1021 if (graph->mirrors[i]) graph->special_dof[i] = PCBDDCGRAPH_LOCAL_PERIODIC_MARK; 1022 1023 if (graph->xadj) { 1024 PetscInt *new_xadj, *new_adjncy; 1025 /* sort CSR graph */ 1026 for (i = 0; i < graph->nvtxs; i++) PetscCall(PetscSortInt(graph->xadj[i + 1] - graph->xadj[i], &graph->adjncy[graph->xadj[i]])); 1027 /* adapt local CSR graph in case of local periodicity */ 1028 k = 0; 1029 for (i = 0; i < graph->nvtxs; i++) 1030 for (j = graph->xadj[i]; j < graph->xadj[i + 1]; j++) k += graph->mirrors[graph->adjncy[j]]; 1031 1032 PetscCall(PetscMalloc1(graph->nvtxs + 1, &new_xadj)); 1033 PetscCall(PetscMalloc1(k + graph->xadj[graph->nvtxs], &new_adjncy)); 1034 new_xadj[0] = 0; 1035 for (i = 0; i < graph->nvtxs; i++) { 1036 k = graph->xadj[i + 1] - graph->xadj[i]; 1037 PetscCall(PetscArraycpy(&new_adjncy[new_xadj[i]], &graph->adjncy[graph->xadj[i]], k)); 1038 new_xadj[i + 1] = new_xadj[i] + k; 1039 for (j = graph->xadj[i]; j < graph->xadj[i + 1]; j++) { 1040 k = graph->mirrors[graph->adjncy[j]]; 1041 PetscCall(PetscArraycpy(&new_adjncy[new_xadj[i + 1]], graph->mirrors_set[graph->adjncy[j]], k)); 1042 new_xadj[i + 1] += k; 1043 } 1044 k = new_xadj[i + 1] - new_xadj[i]; 1045 PetscCall(PetscSortRemoveDupsInt(&k, &new_adjncy[new_xadj[i]])); 1046 new_xadj[i + 1] = new_xadj[i] + k; 1047 } 1048 /* set new CSR into graph */ 1049 PetscCall(PetscFree(graph->xadj)); 1050 PetscCall(PetscFree(graph->adjncy)); 1051 graph->xadj = new_xadj; 1052 graph->adjncy = new_adjncy; 1053 } 1054 } 1055 1056 /* mark special nodes (if any) -> each will become a single node equivalence class */ 1057 if (custom_primal_vertices) { 1058 PetscCall(ISGetLocalSize(custom_primal_vertices, &is_size)); 1059 PetscCall(ISGetIndices(custom_primal_vertices, (const PetscInt **)&is_indices)); 1060 for (i = 0, j = 0; i < is_size; i++) { 1061 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 */ 1062 graph->special_dof[is_indices[i]] = PCBDDCGRAPH_SPECIAL_MARK - j; 1063 j++; 1064 } 1065 } 1066 PetscCall(ISRestoreIndices(custom_primal_vertices, (const PetscInt **)&is_indices)); 1067 } 1068 1069 /* mark interior nodes (if commsize > graph->commsizelimit) as touched and belonging to partition number 0 */ 1070 if (commsize > graph->commsizelimit) { 1071 for (i = 0; i < graph->nvtxs; i++) { 1072 if (!graph->count[i]) { 1073 PetscCall(PetscBTSet(graph->touched, i)); 1074 graph->subset[i] = 0; 1075 } 1076 } 1077 } 1078 1079 /* init graph structure and compute default subsets */ 1080 nodes_touched = 0; 1081 for (i = 0; i < graph->nvtxs; i++) { 1082 if (PetscBTLookup(graph->touched, i)) nodes_touched++; 1083 } 1084 i = 0; 1085 graph->ncc = 0; 1086 total_counts = 0; 1087 1088 /* allocated space for queues */ 1089 if (commsize == graph->commsizelimit) { 1090 PetscCall(PetscMalloc2(graph->nvtxs + 1, &graph->cptr, graph->nvtxs, &graph->queue)); 1091 } else { 1092 PetscInt nused = graph->nvtxs - nodes_touched; 1093 PetscCall(PetscMalloc2(nused + 1, &graph->cptr, nused, &graph->queue)); 1094 } 1095 1096 while (nodes_touched < graph->nvtxs) { 1097 /* find first untouched node in local ordering */ 1098 while (PetscBTLookup(graph->touched, i)) i++; 1099 PetscCall(PetscBTSet(graph->touched, i)); 1100 graph->subset[i] = graph->ncc + 1; 1101 graph->cptr[graph->ncc] = total_counts; 1102 graph->queue[total_counts] = i; 1103 total_counts++; 1104 nodes_touched++; 1105 /* now find all other nodes having the same set of sharing subdomains */ 1106 for (j = i + 1; j < graph->nvtxs; j++) { 1107 /* check for same number of sharing subdomains, dof number and same special mark */ 1108 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]) { 1109 /* check for same set of sharing subdomains */ 1110 same_set = PETSC_TRUE; 1111 for (k = 0; k < graph->count[j]; k++) { 1112 if (graph->neighbours_set[i][k] != graph->neighbours_set[j][k]) same_set = PETSC_FALSE; 1113 } 1114 /* I have found a friend of mine */ 1115 if (same_set) { 1116 PetscCall(PetscBTSet(graph->touched, j)); 1117 graph->subset[j] = graph->ncc + 1; 1118 nodes_touched++; 1119 graph->queue[total_counts] = j; 1120 total_counts++; 1121 } 1122 } 1123 } 1124 graph->ncc++; 1125 } 1126 /* 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 */ 1127 graph->n_subsets = graph->ncc; 1128 PetscCall(PetscMalloc1(graph->n_subsets, &graph->subset_ncc)); 1129 for (i = 0; i < graph->n_subsets; i++) graph->subset_ncc[i] = 1; 1130 /* final pointer */ 1131 graph->cptr[graph->ncc] = total_counts; 1132 1133 /* For consistency reasons (among neighbours), I need to sort (by global ordering) each connected component */ 1134 /* Get a reference node (min index in global ordering) for each subset for tagging messages */ 1135 PetscCall(PetscMalloc1(graph->ncc, &graph->subset_ref_node)); 1136 PetscCall(PetscMalloc1(graph->cptr[graph->ncc], &queue_global)); 1137 PetscCall(ISLocalToGlobalMappingApply(graph->l2gmap, graph->cptr[graph->ncc], graph->queue, queue_global)); 1138 for (j = 0; j < graph->ncc; j++) { 1139 PetscCall(PetscSortIntWithArray(graph->cptr[j + 1] - graph->cptr[j], &queue_global[graph->cptr[j]], &graph->queue[graph->cptr[j]])); 1140 graph->subset_ref_node[j] = graph->queue[graph->cptr[j]]; 1141 } 1142 PetscCall(PetscFree(queue_global)); 1143 graph->queue_sorted = PETSC_TRUE; 1144 1145 /* save information on subsets (needed when analyzing the connected components) */ 1146 if (graph->ncc) { 1147 PetscCall(PetscMalloc2(graph->ncc, &graph->subset_size, graph->ncc, &graph->subset_idxs)); 1148 PetscCall(PetscMalloc1(graph->cptr[graph->ncc], &graph->subset_idxs[0])); 1149 PetscCall(PetscArrayzero(graph->subset_idxs[0], graph->cptr[graph->ncc])); 1150 for (j = 1; j < graph->ncc; j++) { 1151 graph->subset_size[j - 1] = graph->cptr[j] - graph->cptr[j - 1]; 1152 graph->subset_idxs[j] = graph->subset_idxs[j - 1] + graph->subset_size[j - 1]; 1153 } 1154 graph->subset_size[graph->ncc - 1] = graph->cptr[graph->ncc] - graph->cptr[graph->ncc - 1]; 1155 PetscCall(PetscArraycpy(graph->subset_idxs[0], graph->queue, graph->cptr[graph->ncc])); 1156 } 1157 1158 /* renumber reference nodes */ 1159 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)graph->l2gmap), graph->ncc, graph->subset_ref_node, PETSC_COPY_VALUES, &subset_n)); 1160 PetscCall(ISLocalToGlobalMappingApplyIS(graph->l2gmap, subset_n, &subset)); 1161 PetscCall(ISDestroy(&subset_n)); 1162 PetscCall(ISRenumber(subset, NULL, NULL, &subset_n)); 1163 PetscCall(ISDestroy(&subset)); 1164 PetscCall(ISGetLocalSize(subset_n, &k)); 1165 PetscCheck(k == graph->ncc, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid size of new subset! %" PetscInt_FMT " != %" PetscInt_FMT, k, graph->ncc); 1166 PetscCall(ISGetIndices(subset_n, &is_indices)); 1167 PetscCall(PetscArraycpy(graph->subset_ref_node, is_indices, graph->ncc)); 1168 PetscCall(ISRestoreIndices(subset_n, &is_indices)); 1169 PetscCall(ISDestroy(&subset_n)); 1170 1171 /* free workspace */ 1172 graph->setupcalled = PETSC_TRUE; 1173 PetscFunctionReturn(PETSC_SUCCESS); 1174 } 1175 1176 PetscErrorCode PCBDDCGraphResetCoords(PCBDDCGraph graph) 1177 { 1178 PetscFunctionBegin; 1179 if (!graph) PetscFunctionReturn(PETSC_SUCCESS); 1180 PetscCall(PetscFree(graph->coords)); 1181 graph->cdim = 0; 1182 graph->cnloc = 0; 1183 graph->cloc = PETSC_FALSE; 1184 PetscFunctionReturn(PETSC_SUCCESS); 1185 } 1186 1187 PetscErrorCode PCBDDCGraphResetCSR(PCBDDCGraph graph) 1188 { 1189 PetscFunctionBegin; 1190 if (!graph) PetscFunctionReturn(PETSC_SUCCESS); 1191 if (graph->freecsr) { 1192 PetscCall(PetscFree(graph->xadj)); 1193 PetscCall(PetscFree(graph->adjncy)); 1194 } else { 1195 graph->xadj = NULL; 1196 graph->adjncy = NULL; 1197 } 1198 graph->freecsr = PETSC_FALSE; 1199 graph->nvtxs_csr = 0; 1200 PetscFunctionReturn(PETSC_SUCCESS); 1201 } 1202 1203 PetscErrorCode PCBDDCGraphReset(PCBDDCGraph graph) 1204 { 1205 PetscFunctionBegin; 1206 if (!graph) PetscFunctionReturn(PETSC_SUCCESS); 1207 PetscCall(ISLocalToGlobalMappingDestroy(&graph->l2gmap)); 1208 PetscCall(PetscFree(graph->subset_ncc)); 1209 PetscCall(PetscFree(graph->subset_ref_node)); 1210 if (graph->nvtxs) PetscCall(PetscFree(graph->neighbours_set[0])); 1211 PetscCall(PetscBTDestroy(&graph->touched)); 1212 PetscCall(PetscFree5(graph->count, graph->neighbours_set, graph->subset, graph->which_dof, graph->special_dof)); 1213 PetscCall(PetscFree2(graph->cptr, graph->queue)); 1214 if (graph->mirrors) PetscCall(PetscFree(graph->mirrors_set[0])); 1215 PetscCall(PetscFree2(graph->mirrors, graph->mirrors_set)); 1216 if (graph->subset_idxs) PetscCall(PetscFree(graph->subset_idxs[0])); 1217 PetscCall(PetscFree2(graph->subset_size, graph->subset_idxs)); 1218 PetscCall(ISDestroy(&graph->dirdofs)); 1219 PetscCall(ISDestroy(&graph->dirdofsB)); 1220 if (graph->n_local_subs) PetscCall(PetscFree(graph->local_subs)); 1221 graph->has_dirichlet = PETSC_FALSE; 1222 graph->twodimset = PETSC_FALSE; 1223 graph->twodim = PETSC_FALSE; 1224 graph->nvtxs = 0; 1225 graph->nvtxs_global = 0; 1226 graph->n_subsets = 0; 1227 graph->custom_minimal_size = 1; 1228 graph->n_local_subs = 0; 1229 graph->maxcount = PETSC_MAX_INT; 1230 graph->setupcalled = PETSC_FALSE; 1231 PetscFunctionReturn(PETSC_SUCCESS); 1232 } 1233 1234 PetscErrorCode PCBDDCGraphInit(PCBDDCGraph graph, ISLocalToGlobalMapping l2gmap, PetscInt N, PetscInt maxcount) 1235 { 1236 PetscInt n; 1237 1238 PetscFunctionBegin; 1239 PetscAssertPointer(graph, 1); 1240 PetscValidHeaderSpecific(l2gmap, IS_LTOGM_CLASSID, 2); 1241 PetscValidLogicalCollectiveInt(l2gmap, N, 3); 1242 PetscValidLogicalCollectiveInt(l2gmap, maxcount, 4); 1243 /* raise an error if already allocated */ 1244 PetscCheck(!graph->nvtxs_global, PetscObjectComm((PetscObject)l2gmap), PETSC_ERR_PLIB, "BDDCGraph already initialized"); 1245 /* set number of vertices */ 1246 PetscCall(PetscObjectReference((PetscObject)l2gmap)); 1247 graph->l2gmap = l2gmap; 1248 PetscCall(ISLocalToGlobalMappingGetSize(l2gmap, &n)); 1249 graph->nvtxs = n; 1250 graph->nvtxs_global = N; 1251 /* allocate used space */ 1252 PetscCall(PetscBTCreate(graph->nvtxs, &graph->touched)); 1253 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)); 1254 /* zeroes memory */ 1255 PetscCall(PetscArrayzero(graph->count, graph->nvtxs)); 1256 PetscCall(PetscArrayzero(graph->subset, graph->nvtxs)); 1257 /* use -1 as a default value for which_dof array */ 1258 for (n = 0; n < graph->nvtxs; n++) graph->which_dof[n] = -1; 1259 PetscCall(PetscArrayzero(graph->special_dof, graph->nvtxs)); 1260 /* zeroes first pointer to neighbour set */ 1261 if (graph->nvtxs) graph->neighbours_set[0] = NULL; 1262 /* zeroes workspace for values of ncc */ 1263 graph->subset_ncc = NULL; 1264 graph->subset_ref_node = NULL; 1265 /* maxcount for cc */ 1266 graph->maxcount = maxcount; 1267 PetscFunctionReturn(PETSC_SUCCESS); 1268 } 1269 1270 PetscErrorCode PCBDDCGraphDestroy(PCBDDCGraph *graph) 1271 { 1272 PetscFunctionBegin; 1273 PetscCall(PCBDDCGraphResetCSR(*graph)); 1274 PetscCall(PCBDDCGraphResetCoords(*graph)); 1275 PetscCall(PCBDDCGraphReset(*graph)); 1276 PetscCall(PetscFree(*graph)); 1277 PetscFunctionReturn(PETSC_SUCCESS); 1278 } 1279 1280 PetscErrorCode PCBDDCGraphCreate(PCBDDCGraph *graph) 1281 { 1282 PCBDDCGraph new_graph; 1283 1284 PetscFunctionBegin; 1285 PetscCall(PetscNew(&new_graph)); 1286 new_graph->custom_minimal_size = 1; 1287 new_graph->commsizelimit = 1; 1288 *graph = new_graph; 1289 PetscFunctionReturn(PETSC_SUCCESS); 1290 } 1291