1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/isimpl.h> 3 #include <petscsf.h> 4 #include <petscds.h> 5 6 /* get adjacencies due to point-to-point constraints that can't be found with DMPlexGetAdjacency() */ 7 static PetscErrorCode DMPlexComputeAnchorAdjacencies(DM dm, PetscBool useCone, PetscBool useClosure, PetscSection *anchorSectionAdj, PetscInt *anchorAdj[]) 8 { 9 PetscInt pStart, pEnd; 10 PetscSection section, sectionGlobal, adjSec, aSec; 11 IS aIS; 12 13 PetscFunctionBegin; 14 PetscCall(DMGetLocalSection(dm, §ion)); 15 PetscCall(DMGetGlobalSection(dm, §ionGlobal)); 16 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section), &adjSec)); 17 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 18 PetscCall(PetscSectionSetChart(adjSec, pStart, pEnd)); 19 20 PetscCall(DMPlexGetAnchors(dm, &aSec, &aIS)); 21 if (aSec) { 22 const PetscInt *anchors; 23 PetscInt p, q, a, aSize, *offsets, aStart, aEnd, *inverse, iSize, *adj, adjSize; 24 PetscInt *tmpAdjP = NULL, *tmpAdjQ = NULL; 25 PetscSection inverseSec; 26 27 /* invert the constraint-to-anchor map */ 28 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)aSec), &inverseSec)); 29 PetscCall(PetscSectionSetChart(inverseSec, pStart, pEnd)); 30 PetscCall(ISGetLocalSize(aIS, &aSize)); 31 PetscCall(ISGetIndices(aIS, &anchors)); 32 33 for (p = 0; p < aSize; p++) { 34 PetscInt a = anchors[p]; 35 36 PetscCall(PetscSectionAddDof(inverseSec, a, 1)); 37 } 38 PetscCall(PetscSectionSetUp(inverseSec)); 39 PetscCall(PetscSectionGetStorageSize(inverseSec, &iSize)); 40 PetscCall(PetscMalloc1(iSize, &inverse)); 41 PetscCall(PetscCalloc1(pEnd - pStart, &offsets)); 42 PetscCall(PetscSectionGetChart(aSec, &aStart, &aEnd)); 43 for (p = aStart; p < aEnd; p++) { 44 PetscInt dof, off; 45 46 PetscCall(PetscSectionGetDof(aSec, p, &dof)); 47 PetscCall(PetscSectionGetOffset(aSec, p, &off)); 48 49 for (q = 0; q < dof; q++) { 50 PetscInt iOff; 51 52 a = anchors[off + q]; 53 PetscCall(PetscSectionGetOffset(inverseSec, a, &iOff)); 54 inverse[iOff + offsets[a - pStart]++] = p; 55 } 56 } 57 PetscCall(ISRestoreIndices(aIS, &anchors)); 58 PetscCall(PetscFree(offsets)); 59 60 /* construct anchorAdj and adjSec 61 * 62 * loop over anchors: 63 * construct anchor adjacency 64 * loop over constrained: 65 * construct constrained adjacency 66 * if not in anchor adjacency, add to dofs 67 * setup adjSec, allocate anchorAdj 68 * loop over anchors: 69 * construct anchor adjacency 70 * loop over constrained: 71 * construct constrained adjacency 72 * if not in anchor adjacency 73 * if not already in list, put in list 74 * sort, unique, reduce dof count 75 * optional: compactify 76 */ 77 for (p = pStart; p < pEnd; p++) { 78 PetscInt iDof, iOff, i, r, s, numAdjP = PETSC_DETERMINE; 79 80 PetscCall(PetscSectionGetDof(inverseSec, p, &iDof)); 81 if (!iDof) continue; 82 PetscCall(PetscSectionGetOffset(inverseSec, p, &iOff)); 83 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, PETSC_TRUE, &numAdjP, &tmpAdjP)); 84 for (i = 0; i < iDof; i++) { 85 PetscInt iNew = 0, qAdj, qAdjDof, qAdjCDof, numAdjQ = PETSC_DETERMINE; 86 87 q = inverse[iOff + i]; 88 PetscCall(DMPlexGetAdjacency_Internal(dm, q, useCone, useClosure, PETSC_TRUE, &numAdjQ, &tmpAdjQ)); 89 for (r = 0; r < numAdjQ; r++) { 90 qAdj = tmpAdjQ[r]; 91 if ((qAdj < pStart) || (qAdj >= pEnd)) continue; 92 for (s = 0; s < numAdjP; s++) { 93 if (qAdj == tmpAdjP[s]) break; 94 } 95 if (s < numAdjP) continue; 96 PetscCall(PetscSectionGetDof(section, qAdj, &qAdjDof)); 97 PetscCall(PetscSectionGetConstraintDof(section, qAdj, &qAdjCDof)); 98 iNew += qAdjDof - qAdjCDof; 99 } 100 PetscCall(PetscSectionAddDof(adjSec, p, iNew)); 101 } 102 } 103 104 PetscCall(PetscSectionSetUp(adjSec)); 105 PetscCall(PetscSectionGetStorageSize(adjSec, &adjSize)); 106 PetscCall(PetscMalloc1(adjSize, &adj)); 107 108 for (p = pStart; p < pEnd; p++) { 109 PetscInt iDof, iOff, i, r, s, aOff, aOffOrig, aDof, numAdjP = PETSC_DETERMINE; 110 111 PetscCall(PetscSectionGetDof(inverseSec, p, &iDof)); 112 if (!iDof) continue; 113 PetscCall(PetscSectionGetOffset(inverseSec, p, &iOff)); 114 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, PETSC_TRUE, &numAdjP, &tmpAdjP)); 115 PetscCall(PetscSectionGetDof(adjSec, p, &aDof)); 116 PetscCall(PetscSectionGetOffset(adjSec, p, &aOff)); 117 aOffOrig = aOff; 118 for (i = 0; i < iDof; i++) { 119 PetscInt qAdj, qAdjDof, qAdjCDof, qAdjOff, nd, numAdjQ = PETSC_DETERMINE; 120 121 q = inverse[iOff + i]; 122 PetscCall(DMPlexGetAdjacency_Internal(dm, q, useCone, useClosure, PETSC_TRUE, &numAdjQ, &tmpAdjQ)); 123 for (r = 0; r < numAdjQ; r++) { 124 qAdj = tmpAdjQ[r]; 125 if ((qAdj < pStart) || (qAdj >= pEnd)) continue; 126 for (s = 0; s < numAdjP; s++) { 127 if (qAdj == tmpAdjP[s]) break; 128 } 129 if (s < numAdjP) continue; 130 PetscCall(PetscSectionGetDof(section, qAdj, &qAdjDof)); 131 PetscCall(PetscSectionGetConstraintDof(section, qAdj, &qAdjCDof)); 132 PetscCall(PetscSectionGetOffset(sectionGlobal, qAdj, &qAdjOff)); 133 for (nd = 0; nd < qAdjDof - qAdjCDof; ++nd) adj[aOff++] = (qAdjOff < 0 ? -(qAdjOff + 1) : qAdjOff) + nd; 134 } 135 } 136 PetscCall(PetscSortRemoveDupsInt(&aDof, PetscSafePointerPlusOffset(adj, aOffOrig))); 137 PetscCall(PetscSectionSetDof(adjSec, p, aDof)); 138 } 139 *anchorAdj = adj; 140 141 /* clean up */ 142 PetscCall(PetscSectionDestroy(&inverseSec)); 143 PetscCall(PetscFree(inverse)); 144 PetscCall(PetscFree(tmpAdjP)); 145 PetscCall(PetscFree(tmpAdjQ)); 146 } else { 147 *anchorAdj = NULL; 148 PetscCall(PetscSectionSetUp(adjSec)); 149 } 150 *anchorSectionAdj = adjSec; 151 PetscFunctionReturn(PETSC_SUCCESS); 152 } 153 154 // Determine if any of the local adjancencies match a leaf and root of the pointSF. 155 // When using isoperiodic boundary conditions, it is possible for periodic (leaf) and donor (root) pairs to be on the same rank. 156 // This check is done to ensure the adjancency in these cases is only counted for one of the mesh points rather than both. 157 static inline PetscErrorCode AdjancencyContainsLeafRootPair(PetscInt num_pairs, const PetscInt leaves[], const PetscInt roots[], PetscInt numAdj, const PetscInt tmpAdj[], PetscInt *num_leaves_found, PetscInt leaves_found[]) 158 { 159 PetscInt root_index = -1, leaf_, num_roots_found = 0; 160 161 PetscFunctionBeginUser; 162 *num_leaves_found = 0; 163 for (PetscInt q = 0; q < numAdj; q++) { 164 const PetscInt padj = tmpAdj[q]; 165 PetscCall(PetscFindInt(padj, num_pairs, roots, &root_index)); 166 if (root_index >= 0) { 167 leaves_found[num_roots_found] = root_index; // Initially use leaves_found to store pair indices 168 num_roots_found++; 169 break; 170 } 171 } 172 if (num_roots_found == 0) PetscFunctionReturn(PETSC_SUCCESS); 173 174 for (PetscInt i = 0; i < num_roots_found; i++) { 175 leaf_ = leaves[root_index]; 176 for (PetscInt q = 0; q < numAdj; q++) { 177 if (tmpAdj[q] == leaf_) { 178 leaves_found[*num_leaves_found] = leaf_; 179 (*num_leaves_found)++; 180 continue; 181 } 182 } 183 } 184 185 PetscCall(PetscIntSortSemiOrdered(*num_leaves_found, leaves_found)); 186 PetscFunctionReturn(PETSC_SUCCESS); 187 } 188 189 static PetscErrorCode DMPlexCreateAdjacencySection_Static(DM dm, PetscInt bs, PetscSF sfDof, PetscBool useCone, PetscBool useClosure, PetscBool useAnchors, PetscSection *sA, PetscInt **colIdx) 190 { 191 MPI_Comm comm; 192 PetscMPIInt myrank; 193 PetscBool doCommLocal, doComm, debug = PETSC_FALSE; 194 PetscSF sf, sfAdj; 195 PetscSection section, sectionGlobal, leafSectionAdj, rootSectionAdj, sectionAdj, anchorSectionAdj; 196 PetscInt nroots, nleaves, l, p, r; 197 const PetscInt *leaves; 198 const PetscSFNode *remotes; 199 PetscInt dim, pStart, pEnd, numDof, globalOffStart, globalOffEnd, numCols; 200 PetscInt *tmpAdj = NULL, *adj, *rootAdj, *anchorAdj = NULL, *cols, *remoteOffsets, *rootsMyRankPair = NULL, *leavesMyRankPair = NULL; 201 PetscInt adjSize, numMyRankPair = 0; 202 203 PetscFunctionBegin; 204 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 205 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_view_preallocation", &debug, NULL)); 206 PetscCallMPI(MPI_Comm_rank(comm, &myrank)); 207 PetscCall(DMGetDimension(dm, &dim)); 208 PetscCall(DMGetIsoperiodicPointSF_Internal(dm, &sf)); 209 PetscCall(DMGetLocalSection(dm, §ion)); 210 PetscCall(DMGetGlobalSection(dm, §ionGlobal)); 211 PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL)); 212 doCommLocal = nroots >= 0 ? PETSC_TRUE : PETSC_FALSE; 213 PetscCallMPI(MPIU_Allreduce(&doCommLocal, &doComm, 1, MPI_C_BOOL, MPI_LAND, comm)); 214 /* Create section for dof adjacency (dof ==> # adj dof) */ 215 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 216 PetscCall(PetscSectionGetStorageSize(section, &numDof)); 217 PetscCall(PetscSectionCreate(comm, &leafSectionAdj)); 218 PetscCall(PetscSectionSetChart(leafSectionAdj, 0, numDof)); 219 PetscCall(PetscSectionCreate(comm, &rootSectionAdj)); 220 PetscCall(PetscSectionSetChart(rootSectionAdj, 0, numDof)); 221 PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, &remotes)); 222 223 // Store leaf-root pairs if remote.rank is the current rank 224 if (nleaves >= 0) PetscCall(PetscMalloc2(nleaves, &rootsMyRankPair, nleaves, &leavesMyRankPair)); 225 for (PetscInt l = 0; l < nleaves; l++) { 226 if (remotes[l].rank == myrank) { 227 rootsMyRankPair[numMyRankPair] = remotes[l].index; 228 leavesMyRankPair[numMyRankPair] = leaves[l]; 229 numMyRankPair++; 230 } 231 } 232 PetscCall(PetscIntSortSemiOrdered(numMyRankPair, rootsMyRankPair)); 233 PetscCall(PetscIntSortSemiOrdered(numMyRankPair, leavesMyRankPair)); 234 if (debug) { 235 PetscCall(PetscPrintf(comm, "Roots on the same rank:\n")); 236 PetscCall(PetscIntView(numMyRankPair, rootsMyRankPair, NULL)); 237 } 238 /* 239 section - maps points to (# dofs, local dofs) 240 sectionGlobal - maps points to (# dofs, global dofs) 241 leafSectionAdj - maps unowned local dofs to # adj dofs 242 rootSectionAdj - maps owned local dofs to # adj dofs 243 adj - adj global dofs indexed by leafSectionAdj 244 rootAdj - adj global dofs indexed by rootSectionAdj 245 sf - describes shared points across procs 246 sfDof - describes shared dofs across procs 247 sfAdj - describes shared adjacent dofs across procs 248 ** The bootstrapping process involves six rounds with similar structure of visiting neighbors of each point. 249 (0). If there are point-to-point constraints, add the adjacencies of constrained points to anchors in anchorAdj 250 (This is done in DMPlexComputeAnchorAdjacencies()) 251 1. Visit unowned points on interface, count adjacencies placing in leafSectionAdj 252 Reduce those counts to rootSectionAdj (now redundantly counting some interface points) 253 2. Visit owned points on interface, count adjacencies placing in rootSectionAdj 254 Create sfAdj connecting rootSectionAdj and leafSectionAdj 255 3. Visit unowned points on interface, write adjacencies to adj 256 Gather adj to rootAdj (note that there is redundancy in rootAdj when multiple procs find the same adjacencies) 257 4. Visit owned points on interface, write adjacencies to rootAdj 258 Remove redundancy in rootAdj 259 ** The last two traversals use transitive closure 260 5. Visit all owned points in the subdomain, count dofs for each point (sectionAdj) 261 Allocate memory addressed by sectionAdj (cols) 262 6. Visit all owned points in the subdomain, insert dof adjacencies into cols 263 ** Knowing all the column adjacencies, check ownership and sum into dnz and onz 264 */ 265 PetscCall(DMPlexComputeAnchorAdjacencies(dm, useCone, useClosure, &anchorSectionAdj, &anchorAdj)); 266 for (l = 0; l < nleaves; ++l) { 267 PetscInt dof, off, d, q, anDof; 268 PetscInt p = leaves[l], numAdj = PETSC_DETERMINE; 269 270 if ((p < pStart) || (p >= pEnd)) continue; 271 PetscCall(PetscSectionGetDof(section, p, &dof)); 272 PetscCall(PetscSectionGetOffset(section, p, &off)); 273 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 274 for (q = 0; q < numAdj; ++q) { 275 const PetscInt padj = tmpAdj[q]; 276 PetscInt ndof, ncdof; 277 278 if ((padj < pStart) || (padj >= pEnd)) continue; 279 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 280 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 281 for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(leafSectionAdj, d, ndof - ncdof)); 282 } 283 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 284 if (anDof) { 285 for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(leafSectionAdj, d, anDof)); 286 } 287 } 288 PetscCall(PetscSectionSetUp(leafSectionAdj)); 289 if (debug) { 290 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Leaves:\n")); 291 PetscCall(PetscSectionView(leafSectionAdj, NULL)); 292 } 293 /* Get maximum remote adjacency sizes for owned dofs on interface (roots) */ 294 if (doComm) { 295 PetscCall(PetscSFReduceBegin(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM)); 296 PetscCall(PetscSFReduceEnd(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM)); 297 PetscCall(PetscSectionInvalidateMaxDof_Internal(rootSectionAdj)); 298 } 299 if (debug) { 300 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots:\n")); 301 PetscCall(PetscSectionView(rootSectionAdj, NULL)); 302 } 303 /* Add in local adjacency sizes for owned dofs on interface (roots) */ 304 for (p = pStart; p < pEnd; ++p) { 305 PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof; 306 307 PetscCall(PetscSectionGetDof(section, p, &dof)); 308 PetscCall(PetscSectionGetOffset(section, p, &off)); 309 if (!dof) continue; 310 PetscCall(PetscSectionGetDof(rootSectionAdj, off, &adof)); 311 if (adof <= 0) continue; 312 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 313 for (q = 0; q < numAdj; ++q) { 314 const PetscInt padj = tmpAdj[q]; 315 PetscInt ndof, ncdof; 316 317 if ((padj < pStart) || (padj >= pEnd)) continue; 318 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 319 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 320 for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(rootSectionAdj, d, ndof - ncdof)); 321 } 322 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 323 if (anDof) { 324 for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(rootSectionAdj, d, anDof)); 325 } 326 } 327 PetscCall(PetscSectionSetUp(rootSectionAdj)); 328 if (debug) { 329 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots after local additions:\n")); 330 PetscCall(PetscSectionView(rootSectionAdj, NULL)); 331 } 332 /* Create adj SF based on dof SF */ 333 PetscCall(PetscSFCreateRemoteOffsets(sfDof, rootSectionAdj, leafSectionAdj, &remoteOffsets)); 334 PetscCall(PetscSFCreateSectionSF(sfDof, rootSectionAdj, remoteOffsets, leafSectionAdj, &sfAdj)); 335 PetscCall(PetscFree(remoteOffsets)); 336 if (debug) { 337 PetscCall(PetscPrintf(comm, "Adjacency SF for Preallocation:\n")); 338 PetscCall(PetscSFView(sfAdj, NULL)); 339 } 340 /* Create leaf adjacency */ 341 PetscCall(PetscSectionSetUp(leafSectionAdj)); 342 PetscCall(PetscSectionGetStorageSize(leafSectionAdj, &adjSize)); 343 PetscCall(PetscCalloc1(adjSize, &adj)); 344 for (l = 0; l < nleaves; ++l) { 345 PetscInt dof, off, d, q, anDof, anOff; 346 PetscInt p = leaves[l], numAdj = PETSC_DETERMINE; 347 348 if ((p < pStart) || (p >= pEnd)) continue; 349 PetscCall(PetscSectionGetDof(section, p, &dof)); 350 PetscCall(PetscSectionGetOffset(section, p, &off)); 351 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 352 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 353 PetscCall(PetscSectionGetOffset(anchorSectionAdj, p, &anOff)); 354 for (d = off; d < off + dof; ++d) { 355 PetscInt aoff, i = 0; 356 357 PetscCall(PetscSectionGetOffset(leafSectionAdj, d, &aoff)); 358 for (q = 0; q < numAdj; ++q) { 359 const PetscInt padj = tmpAdj[q]; 360 PetscInt ndof, ncdof, ngoff, nd; 361 362 if ((padj < pStart) || (padj >= pEnd)) continue; 363 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 364 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 365 PetscCall(PetscSectionGetOffset(sectionGlobal, padj, &ngoff)); 366 for (nd = 0; nd < ndof - ncdof; ++nd) { 367 adj[aoff + i] = (ngoff < 0 ? -(ngoff + 1) : ngoff) + nd; 368 ++i; 369 } 370 } 371 for (q = 0; q < anDof; q++) { 372 adj[aoff + i] = anchorAdj[anOff + q]; 373 ++i; 374 } 375 } 376 } 377 /* Debugging */ 378 if (debug) { 379 PetscCall(PetscPrintf(comm, "Leaf adjacency indices\n")); 380 PetscCall(PetscSectionArrayView(leafSectionAdj, adj, PETSC_INT, NULL)); 381 } 382 /* Gather adjacent indices to root */ 383 PetscCall(PetscSectionGetStorageSize(rootSectionAdj, &adjSize)); 384 PetscCall(PetscMalloc1(adjSize, &rootAdj)); 385 for (r = 0; r < adjSize; ++r) rootAdj[r] = -1; 386 if (doComm) { 387 const PetscInt *indegree; 388 PetscInt *remoteadj, radjsize = 0; 389 390 PetscCall(PetscSFComputeDegreeBegin(sfAdj, &indegree)); 391 PetscCall(PetscSFComputeDegreeEnd(sfAdj, &indegree)); 392 for (p = 0; p < adjSize; ++p) radjsize += indegree[p]; 393 PetscCall(PetscMalloc1(radjsize, &remoteadj)); 394 PetscCall(PetscSFGatherBegin(sfAdj, MPIU_INT, adj, remoteadj)); 395 PetscCall(PetscSFGatherEnd(sfAdj, MPIU_INT, adj, remoteadj)); 396 for (p = 0, l = 0, r = 0; p < adjSize; ++p, l = PetscMax(p, l + indegree[p - 1])) { 397 PetscInt s; 398 for (s = 0; s < indegree[p]; ++s, ++r) rootAdj[l + s] = remoteadj[r]; 399 } 400 PetscCheck(r == radjsize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %" PetscInt_FMT " != %" PetscInt_FMT, r, radjsize); 401 PetscCheck(l == adjSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %" PetscInt_FMT " != %" PetscInt_FMT, l, adjSize); 402 PetscCall(PetscFree(remoteadj)); 403 } 404 PetscCall(PetscSFDestroy(&sfAdj)); 405 PetscCall(PetscFree(adj)); 406 /* Debugging */ 407 if (debug) { 408 PetscCall(PetscPrintf(comm, "Root adjacency indices after gather\n")); 409 PetscCall(PetscSectionArrayView(rootSectionAdj, rootAdj, PETSC_INT, NULL)); 410 } 411 /* Add in local adjacency indices for owned dofs on interface (roots) */ 412 for (p = pStart; p < pEnd; ++p) { 413 PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof, anOff; 414 415 PetscCall(PetscSectionGetDof(section, p, &dof)); 416 PetscCall(PetscSectionGetOffset(section, p, &off)); 417 if (!dof) continue; 418 PetscCall(PetscSectionGetDof(rootSectionAdj, off, &adof)); 419 if (adof <= 0) continue; 420 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 421 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 422 PetscCall(PetscSectionGetOffset(anchorSectionAdj, p, &anOff)); 423 for (d = off; d < off + dof; ++d) { 424 PetscInt adof, aoff, i; 425 426 PetscCall(PetscSectionGetDof(rootSectionAdj, d, &adof)); 427 PetscCall(PetscSectionGetOffset(rootSectionAdj, d, &aoff)); 428 i = adof - 1; 429 for (q = 0; q < anDof; q++) { 430 rootAdj[aoff + i] = anchorAdj[anOff + q]; 431 --i; 432 } 433 for (q = 0; q < numAdj; ++q) { 434 const PetscInt padj = tmpAdj[q]; 435 PetscInt ndof, ncdof, ngoff, nd; 436 437 if ((padj < pStart) || (padj >= pEnd)) continue; 438 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 439 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 440 PetscCall(PetscSectionGetOffset(sectionGlobal, padj, &ngoff)); 441 for (nd = 0; nd < ndof - ncdof; ++nd) { 442 rootAdj[aoff + i] = ngoff < 0 ? -(ngoff + 1) + nd : ngoff + nd; 443 --i; 444 } 445 } 446 } 447 } 448 /* Debugging */ 449 if (debug) { 450 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots before compression:\n")); 451 PetscCall(PetscSectionArrayView(rootSectionAdj, rootAdj, PETSC_INT, NULL)); 452 } 453 /* Compress indices */ 454 PetscCall(PetscSectionSetUp(rootSectionAdj)); 455 for (p = pStart; p < pEnd; ++p) { 456 PetscInt dof, cdof, off, d; 457 PetscInt adof, aoff; 458 459 PetscCall(PetscSectionGetDof(section, p, &dof)); 460 PetscCall(PetscSectionGetConstraintDof(section, p, &cdof)); 461 PetscCall(PetscSectionGetOffset(section, p, &off)); 462 if (!dof) continue; 463 PetscCall(PetscSectionGetDof(rootSectionAdj, off, &adof)); 464 if (adof <= 0) continue; 465 for (d = off; d < off + dof - cdof; ++d) { 466 PetscCall(PetscSectionGetDof(rootSectionAdj, d, &adof)); 467 PetscCall(PetscSectionGetOffset(rootSectionAdj, d, &aoff)); 468 PetscCall(PetscSortRemoveDupsInt(&adof, &rootAdj[aoff])); 469 PetscCall(PetscSectionSetDof(rootSectionAdj, d, adof)); 470 } 471 } 472 /* Debugging */ 473 if (debug) { 474 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots after compression:\n")); 475 PetscCall(PetscSectionArrayView(rootSectionAdj, rootAdj, PETSC_INT, NULL)); 476 } 477 /* Build adjacency section: Maps global indices to sets of adjacent global indices */ 478 PetscCall(PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd)); 479 PetscCall(PetscSectionCreate(comm, §ionAdj)); 480 PetscCall(PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd)); 481 482 PetscInt *exclude_leaves, num_exclude_leaves, max_adjacency_size = 0; 483 PetscCall(DMPlexGetMaxAdjacencySize_Internal(dm, useAnchors, &max_adjacency_size)); 484 PetscCall(PetscMalloc1(max_adjacency_size, &exclude_leaves)); 485 486 for (p = pStart; p < pEnd; ++p) { 487 PetscInt numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof; 488 PetscBool found = PETSC_TRUE; 489 490 PetscCall(PetscSectionGetDof(section, p, &dof)); 491 PetscCall(PetscSectionGetConstraintDof(section, p, &cdof)); 492 PetscCall(PetscSectionGetOffset(section, p, &off)); 493 PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff)); 494 for (d = 0; d < dof - cdof; ++d) { 495 PetscInt ldof, rdof; 496 497 PetscCall(PetscSectionGetDof(leafSectionAdj, off + d, &ldof)); 498 PetscCall(PetscSectionGetDof(rootSectionAdj, off + d, &rdof)); 499 if (ldof > 0) { 500 /* We do not own this point */ 501 } else if (rdof > 0) { 502 PetscCall(PetscSectionSetDof(sectionAdj, goff + d, rdof)); 503 } else { 504 found = PETSC_FALSE; 505 } 506 } 507 if (found) continue; 508 PetscCall(PetscSectionGetDof(section, p, &dof)); 509 PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff)); 510 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 511 PetscCall(AdjancencyContainsLeafRootPair(numMyRankPair, leavesMyRankPair, rootsMyRankPair, numAdj, tmpAdj, &num_exclude_leaves, exclude_leaves)); 512 for (q = 0; q < numAdj; ++q) { 513 const PetscInt padj = tmpAdj[q]; 514 PetscInt ndof, ncdof, noff, count; 515 516 if ((padj < pStart) || (padj >= pEnd)) continue; 517 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 518 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 519 PetscCall(PetscSectionGetOffset(section, padj, &noff)); 520 // If leaf-root pair are both on this rank, only count root 521 PetscCall(PetscFindInt(padj, num_exclude_leaves, exclude_leaves, &count)); 522 if (count >= 0) continue; 523 for (d = goff; d < goff + dof - cdof; ++d) PetscCall(PetscSectionAddDof(sectionAdj, d, ndof - ncdof)); 524 } 525 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 526 if (anDof) { 527 for (d = goff; d < goff + dof - cdof; ++d) PetscCall(PetscSectionAddDof(sectionAdj, d, anDof)); 528 } 529 } 530 PetscCall(PetscSectionSetUp(sectionAdj)); 531 if (debug) { 532 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation:\n")); 533 PetscCall(PetscSectionView(sectionAdj, NULL)); 534 } 535 /* Get adjacent indices */ 536 PetscCall(PetscSectionGetStorageSize(sectionAdj, &numCols)); 537 PetscCall(PetscMalloc1(numCols, &cols)); 538 for (p = pStart; p < pEnd; ++p) { 539 PetscInt numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof, anOff; 540 PetscBool found = PETSC_TRUE; 541 542 PetscCall(PetscSectionGetDof(section, p, &dof)); 543 PetscCall(PetscSectionGetConstraintDof(section, p, &cdof)); 544 PetscCall(PetscSectionGetOffset(section, p, &off)); 545 PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff)); 546 for (d = 0; d < dof - cdof; ++d) { 547 PetscInt ldof, rdof; 548 549 PetscCall(PetscSectionGetDof(leafSectionAdj, off + d, &ldof)); 550 PetscCall(PetscSectionGetDof(rootSectionAdj, off + d, &rdof)); 551 if (ldof > 0) { 552 /* We do not own this point */ 553 } else if (rdof > 0) { 554 PetscInt aoff, roff; 555 556 PetscCall(PetscSectionGetOffset(sectionAdj, goff + d, &aoff)); 557 PetscCall(PetscSectionGetOffset(rootSectionAdj, off + d, &roff)); 558 PetscCall(PetscArraycpy(&cols[aoff], &rootAdj[roff], rdof)); 559 } else { 560 found = PETSC_FALSE; 561 } 562 } 563 if (found) continue; 564 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 565 PetscCall(AdjancencyContainsLeafRootPair(numMyRankPair, leavesMyRankPair, rootsMyRankPair, numAdj, tmpAdj, &num_exclude_leaves, exclude_leaves)); 566 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 567 PetscCall(PetscSectionGetOffset(anchorSectionAdj, p, &anOff)); 568 for (d = goff; d < goff + dof - cdof; ++d) { 569 PetscInt adof, aoff, i = 0; 570 571 PetscCall(PetscSectionGetDof(sectionAdj, d, &adof)); 572 PetscCall(PetscSectionGetOffset(sectionAdj, d, &aoff)); 573 for (q = 0; q < numAdj; ++q) { 574 const PetscInt padj = tmpAdj[q], *ncind; 575 PetscInt ndof, ncdof, ngoff, nd, count; 576 577 /* Adjacent points may not be in the section chart */ 578 if ((padj < pStart) || (padj >= pEnd)) continue; 579 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 580 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 581 PetscCall(PetscSectionGetConstraintIndices(section, padj, &ncind)); 582 PetscCall(PetscSectionGetOffset(sectionGlobal, padj, &ngoff)); 583 // If leaf-root pair are both on this rank, only count root 584 PetscCall(PetscFindInt(padj, num_exclude_leaves, exclude_leaves, &count)); 585 if (count >= 0) continue; 586 for (nd = 0; nd < ndof - ncdof; ++nd, ++i) cols[aoff + i] = ngoff < 0 ? -(ngoff + 1) + nd : ngoff + nd; 587 } 588 for (q = 0; q < anDof; q++, i++) cols[aoff + i] = anchorAdj[anOff + q]; 589 PetscCheck(i == adof, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of entries %" PetscInt_FMT " != %" PetscInt_FMT " for dof %" PetscInt_FMT " (point %" PetscInt_FMT ")", i, adof, d, p); 590 } 591 } 592 PetscCall(PetscSectionDestroy(&anchorSectionAdj)); 593 PetscCall(PetscSectionDestroy(&leafSectionAdj)); 594 PetscCall(PetscSectionDestroy(&rootSectionAdj)); 595 PetscCall(PetscFree(exclude_leaves)); 596 PetscCall(PetscFree2(rootsMyRankPair, leavesMyRankPair)); 597 PetscCall(PetscFree(anchorAdj)); 598 PetscCall(PetscFree(rootAdj)); 599 PetscCall(PetscFree(tmpAdj)); 600 /* Debugging */ 601 if (debug) { 602 PetscCall(PetscPrintf(comm, "Column indices\n")); 603 PetscCall(PetscSectionArrayView(sectionAdj, cols, PETSC_INT, NULL)); 604 } 605 606 *sA = sectionAdj; 607 *colIdx = cols; 608 PetscFunctionReturn(PETSC_SUCCESS); 609 } 610 611 static PetscErrorCode DMPlexUpdateAllocation_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[]) 612 { 613 PetscSection section; 614 PetscInt rStart, rEnd, r, pStart, pEnd, p; 615 616 PetscFunctionBegin; 617 /* This loop needs to change to a loop over points, then field dofs, which means we need to look both sections */ 618 PetscCall(PetscLayoutGetRange(rLayout, &rStart, &rEnd)); 619 PetscCheck((rStart % bs) == 0 && (rEnd % bs) == 0, PetscObjectComm((PetscObject)rLayout), PETSC_ERR_ARG_WRONG, "Invalid layout [%" PetscInt_FMT ", %" PetscInt_FMT ") for matrix, must be divisible by block size %" PetscInt_FMT, rStart, rEnd, bs); 620 if (f >= 0 && bs == 1) { 621 PetscCall(DMGetLocalSection(dm, §ion)); 622 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 623 for (p = pStart; p < pEnd; ++p) { 624 PetscInt rS, rE; 625 626 PetscCall(DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE)); 627 for (r = rS; r < rE; ++r) { 628 PetscInt numCols, cStart, c; 629 630 PetscCall(PetscSectionGetDof(sectionAdj, r, &numCols)); 631 PetscCall(PetscSectionGetOffset(sectionAdj, r, &cStart)); 632 for (c = cStart; c < cStart + numCols; ++c) { 633 if ((cols[c] >= rStart) && (cols[c] < rEnd)) { 634 ++dnz[r - rStart]; 635 if (cols[c] >= r) ++dnzu[r - rStart]; 636 } else { 637 ++onz[r - rStart]; 638 if (cols[c] >= r) ++onzu[r - rStart]; 639 } 640 } 641 } 642 } 643 } else { 644 /* Only loop over blocks of rows */ 645 for (r = rStart / bs; r < rEnd / bs; ++r) { 646 const PetscInt row = r * bs; 647 PetscInt numCols, cStart, c; 648 649 PetscCall(PetscSectionGetDof(sectionAdj, row, &numCols)); 650 PetscCall(PetscSectionGetOffset(sectionAdj, row, &cStart)); 651 for (c = cStart; c < cStart + numCols; ++c) { 652 if ((cols[c] >= rStart) && (cols[c] < rEnd)) { 653 ++dnz[r - rStart / bs]; 654 if (cols[c] >= row) ++dnzu[r - rStart / bs]; 655 } else { 656 ++onz[r - rStart / bs]; 657 if (cols[c] >= row) ++onzu[r - rStart / bs]; 658 } 659 } 660 } 661 for (r = 0; r < (rEnd - rStart) / bs; ++r) { 662 dnz[r] /= bs; 663 onz[r] /= bs; 664 dnzu[r] /= bs; 665 onzu[r] /= bs; 666 } 667 } 668 PetscFunctionReturn(PETSC_SUCCESS); 669 } 670 671 static PetscErrorCode DMPlexFillMatrix_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], Mat A) 672 { 673 PetscSection section; 674 PetscScalar *values; 675 PetscInt rStart, rEnd, r, pStart, pEnd, p, len, maxRowLen = 0; 676 677 PetscFunctionBegin; 678 PetscCall(PetscLayoutGetRange(rLayout, &rStart, &rEnd)); 679 for (r = rStart; r < rEnd; ++r) { 680 PetscCall(PetscSectionGetDof(sectionAdj, r, &len)); 681 maxRowLen = PetscMax(maxRowLen, len); 682 } 683 PetscCall(PetscCalloc1(maxRowLen, &values)); 684 if (f >= 0 && bs == 1) { 685 PetscCall(DMGetLocalSection(dm, §ion)); 686 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 687 for (p = pStart; p < pEnd; ++p) { 688 PetscInt rS, rE; 689 690 PetscCall(DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE)); 691 for (r = rS; r < rE; ++r) { 692 PetscInt numCols, cStart; 693 694 PetscCall(PetscSectionGetDof(sectionAdj, r, &numCols)); 695 PetscCall(PetscSectionGetOffset(sectionAdj, r, &cStart)); 696 PetscCall(MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES)); 697 } 698 } 699 } else { 700 for (r = rStart; r < rEnd; ++r) { 701 PetscInt numCols, cStart; 702 703 PetscCall(PetscSectionGetDof(sectionAdj, r, &numCols)); 704 PetscCall(PetscSectionGetOffset(sectionAdj, r, &cStart)); 705 PetscCall(MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES)); 706 } 707 } 708 PetscCall(PetscFree(values)); 709 PetscFunctionReturn(PETSC_SUCCESS); 710 } 711 712 /*@ 713 DMPlexPreallocateOperator - Calculate the matrix nonzero pattern based upon the information in the `DM`, 714 the `PetscDS` it contains, and the default `PetscSection`. 715 716 Collective 717 718 Input Parameters: 719 + dm - The `DMPLEX` 720 . bs - The matrix blocksize 721 . dnz - An array to hold the number of nonzeros in the diagonal block 722 . onz - An array to hold the number of nonzeros in the off-diagonal block 723 . dnzu - An array to hold the number of nonzeros in the upper triangle of the diagonal block 724 . onzu - An array to hold the number of nonzeros in the upper triangle of the off-diagonal block 725 - fillMatrix - If `PETSC_TRUE`, fill the matrix with zeros 726 727 Output Parameter: 728 . A - The preallocated matrix 729 730 Level: advanced 731 732 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreateMatrix()` 733 @*/ 734 PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix) 735 { 736 MPI_Comm comm; 737 MatType mtype; 738 PetscSF sf, sfDof; 739 PetscSection section; 740 PetscInt *remoteOffsets; 741 PetscSection sectionAdj[4] = {NULL, NULL, NULL, NULL}; 742 PetscInt *cols[4] = {NULL, NULL, NULL, NULL}; 743 PetscBool useCone, useClosure; 744 PetscInt Nf, f, idx, locRows; 745 PetscLayout rLayout; 746 PetscBool isSymBlock, isSymSeqBlock, isSymMPIBlock, debug = PETSC_FALSE; 747 748 PetscFunctionBegin; 749 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 750 PetscValidHeaderSpecific(A, MAT_CLASSID, 7); 751 if (dnz) PetscAssertPointer(dnz, 3); 752 if (onz) PetscAssertPointer(onz, 4); 753 if (dnzu) PetscAssertPointer(dnzu, 5); 754 if (onzu) PetscAssertPointer(onzu, 6); 755 PetscCall(DMGetIsoperiodicPointSF_Internal(dm, &sf)); 756 PetscCall(DMGetLocalSection(dm, §ion)); 757 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_view_preallocation", &debug, NULL)); 758 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 759 PetscCall(PetscLogEventBegin(DMPLEX_Preallocate, dm, 0, 0, 0)); 760 /* Create dof SF based on point SF */ 761 if (debug) { 762 PetscSection section, sectionGlobal; 763 PetscSF sf; 764 765 PetscCall(DMGetIsoperiodicPointSF_Internal(dm, &sf)); 766 PetscCall(DMGetLocalSection(dm, §ion)); 767 PetscCall(DMGetGlobalSection(dm, §ionGlobal)); 768 PetscCall(PetscPrintf(comm, "Input Section for Preallocation:\n")); 769 PetscCall(PetscSectionView(section, NULL)); 770 PetscCall(PetscPrintf(comm, "Input Global Section for Preallocation:\n")); 771 PetscCall(PetscSectionView(sectionGlobal, NULL)); 772 PetscCall(PetscPrintf(comm, "Input SF for Preallocation:\n")); 773 PetscCall(PetscSFView(sf, NULL)); 774 } 775 PetscCall(PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets)); 776 PetscCall(PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof)); 777 PetscCall(PetscFree(remoteOffsets)); 778 if (debug) { 779 PetscCall(PetscPrintf(comm, "Dof SF for Preallocation:\n")); 780 PetscCall(PetscSFView(sfDof, NULL)); 781 } 782 /* Create allocation vectors from adjacency graph */ 783 PetscCall(MatGetLocalSize(A, &locRows, NULL)); 784 PetscCall(PetscLayoutCreate(comm, &rLayout)); 785 PetscCall(PetscLayoutSetLocalSize(rLayout, locRows)); 786 PetscCall(PetscLayoutSetBlockSize(rLayout, 1)); 787 PetscCall(PetscLayoutSetUp(rLayout)); 788 /* There are 4 types of adjacency */ 789 PetscCall(PetscSectionGetNumFields(section, &Nf)); 790 if (Nf < 1 || bs > 1) { 791 PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure)); 792 idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0); 793 PetscCall(DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, §ionAdj[idx], &cols[idx])); 794 PetscCall(DMPlexUpdateAllocation_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu)); 795 } else { 796 for (f = 0; f < Nf; ++f) { 797 PetscCall(DMGetAdjacency(dm, f, &useCone, &useClosure)); 798 idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0); 799 if (!sectionAdj[idx]) PetscCall(DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, §ionAdj[idx], &cols[idx])); 800 PetscCall(DMPlexUpdateAllocation_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu)); 801 } 802 } 803 PetscCall(PetscSFDestroy(&sfDof)); 804 /* Set matrix pattern */ 805 PetscCall(MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu)); 806 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 807 /* Check for symmetric storage */ 808 PetscCall(MatGetType(A, &mtype)); 809 PetscCall(PetscStrcmp(mtype, MATSBAIJ, &isSymBlock)); 810 PetscCall(PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock)); 811 PetscCall(PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock)); 812 if (isSymBlock || isSymSeqBlock || isSymMPIBlock) PetscCall(MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE)); 813 /* Fill matrix with zeros */ 814 if (fillMatrix) { 815 if (Nf < 1 || bs > 1) { 816 PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure)); 817 idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0); 818 PetscCall(DMPlexFillMatrix_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], A)); 819 } else { 820 for (f = 0; f < Nf; ++f) { 821 PetscCall(DMGetAdjacency(dm, f, &useCone, &useClosure)); 822 idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0); 823 PetscCall(DMPlexFillMatrix_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], A)); 824 } 825 } 826 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 827 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 828 } 829 PetscCall(PetscLayoutDestroy(&rLayout)); 830 for (idx = 0; idx < 4; ++idx) { 831 PetscCall(PetscSectionDestroy(§ionAdj[idx])); 832 PetscCall(PetscFree(cols[idx])); 833 } 834 PetscCall(PetscLogEventEnd(DMPLEX_Preallocate, dm, 0, 0, 0)); 835 PetscFunctionReturn(PETSC_SUCCESS); 836 } 837 838 #if 0 839 PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix) 840 { 841 PetscInt *tmpClosure,*tmpAdj,*visits; 842 PetscInt c,cStart,cEnd,pStart,pEnd; 843 844 PetscFunctionBegin; 845 PetscCall(DMGetDimension(dm, &dim)); 846 PetscCall(DMPlexGetDepth(dm, &depth)); 847 PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize)); 848 849 maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1)); 850 851 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 852 npoints = pEnd - pStart; 853 854 PetscCall(PetscMalloc3(maxClosureSize,&tmpClosure,npoints,&lvisits,npoints,&visits)); 855 PetscCall(PetscArrayzero(lvisits,pEnd-pStart)); 856 PetscCall(PetscArrayzero(visits,pEnd-pStart)); 857 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 858 for (c=cStart; c<cEnd; c++) { 859 PetscInt *support = tmpClosure; 860 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support)); 861 for (p=0; p<supportSize; p++) lvisits[support[p]]++; 862 } 863 PetscCall(PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM)); 864 PetscCall(PetscSFReduceEnd (sf,MPIU_INT,lvisits,visits,MPI_SUM)); 865 PetscCall(PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits,MPI_REPLACE)); 866 PetscCall(PetscSFBcastEnd (sf,MPIU_INT,visits,lvisits)); 867 868 PetscCall(PetscSFGetRootRanks()); 869 870 PetscCall(PetscMalloc2(maxClosureSize*maxClosureSize,&cellmat,npoints,&owner)); 871 for (c=cStart; c<cEnd; c++) { 872 PetscCall(PetscArrayzero(cellmat,maxClosureSize*maxClosureSize)); 873 /* 874 Depth-first walk of transitive closure. 875 At each leaf frame f of transitive closure that we see, add 1/visits[f] to each pair (p,q) not marked as done in cellmat. 876 This contribution is added to dnz if owning ranks of p and q match, to onz otherwise. 877 */ 878 } 879 880 PetscCall(PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM)); 881 PetscCall(PetscSFReduceEnd (sf,MPIU_INT,lonz,onz,MPI_SUM)); 882 PetscFunctionReturn(PETSC_SUCCESS); 883 } 884 #endif 885