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 static PetscErrorCode DMPlexCreateAdjacencySection_Static(DM dm, PetscInt bs, PetscSF sfDof, PetscBool useCone, PetscBool useClosure, PetscBool useAnchors, PetscSection *sA, PetscInt **colIdx) 155 { 156 MPI_Comm comm; 157 PetscMPIInt size; 158 PetscBool doCommLocal, doComm, debug = PETSC_FALSE; 159 PetscSF sf, sfAdj; 160 PetscSection section, sectionGlobal, leafSectionAdj, rootSectionAdj, sectionAdj, anchorSectionAdj; 161 PetscInt nroots, nleaves, l, p, r; 162 const PetscInt *leaves; 163 const PetscSFNode *remotes; 164 PetscInt dim, pStart, pEnd, numDof, globalOffStart, globalOffEnd, numCols; 165 PetscInt *tmpAdj = NULL, *adj, *rootAdj, *anchorAdj = NULL, *cols, *remoteOffsets; 166 PetscInt adjSize; 167 168 PetscFunctionBegin; 169 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 170 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_view_preallocation", &debug, NULL)); 171 PetscCallMPI(MPI_Comm_size(comm, &size)); 172 PetscCall(DMGetDimension(dm, &dim)); 173 PetscCall(DMGetPointSF(dm, &sf)); 174 PetscCall(DMGetLocalSection(dm, §ion)); 175 PetscCall(DMGetGlobalSection(dm, §ionGlobal)); 176 PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL)); 177 doCommLocal = (size > 1) && (nroots >= 0) ? PETSC_TRUE : PETSC_FALSE; 178 PetscCall(MPIU_Allreduce(&doCommLocal, &doComm, 1, MPIU_BOOL, MPI_LAND, comm)); 179 /* Create section for dof adjacency (dof ==> # adj dof) */ 180 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 181 PetscCall(PetscSectionGetStorageSize(section, &numDof)); 182 PetscCall(PetscSectionCreate(comm, &leafSectionAdj)); 183 PetscCall(PetscSectionSetChart(leafSectionAdj, 0, numDof)); 184 PetscCall(PetscSectionCreate(comm, &rootSectionAdj)); 185 PetscCall(PetscSectionSetChart(rootSectionAdj, 0, numDof)); 186 /* Fill in the ghost dofs on the interface */ 187 PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, &remotes)); 188 /* 189 section - maps points to (# dofs, local dofs) 190 sectionGlobal - maps points to (# dofs, global dofs) 191 leafSectionAdj - maps unowned local dofs to # adj dofs 192 rootSectionAdj - maps owned local dofs to # adj dofs 193 adj - adj global dofs indexed by leafSectionAdj 194 rootAdj - adj global dofs indexed by rootSectionAdj 195 sf - describes shared points across procs 196 sfDof - describes shared dofs across procs 197 sfAdj - describes shared adjacent dofs across procs 198 ** The bootstrapping process involves six rounds with similar structure of visiting neighbors of each point. 199 (0). If there are point-to-point constraints, add the adjacencies of constrained points to anchors in anchorAdj 200 (This is done in DMPlexComputeAnchorAdjacencies()) 201 1. Visit unowned points on interface, count adjacencies placing in leafSectionAdj 202 Reduce those counts to rootSectionAdj (now redundantly counting some interface points) 203 2. Visit owned points on interface, count adjacencies placing in rootSectionAdj 204 Create sfAdj connecting rootSectionAdj and leafSectionAdj 205 3. Visit unowned points on interface, write adjacencies to adj 206 Gather adj to rootAdj (note that there is redundancy in rootAdj when multiple procs find the same adjacencies) 207 4. Visit owned points on interface, write adjacencies to rootAdj 208 Remove redundancy in rootAdj 209 ** The last two traversals use transitive closure 210 5. Visit all owned points in the subdomain, count dofs for each point (sectionAdj) 211 Allocate memory addressed by sectionAdj (cols) 212 6. Visit all owned points in the subdomain, insert dof adjacencies into cols 213 ** Knowing all the column adjacencies, check ownership and sum into dnz and onz 214 */ 215 PetscCall(DMPlexComputeAnchorAdjacencies(dm, useCone, useClosure, &anchorSectionAdj, &anchorAdj)); 216 for (l = 0; l < nleaves; ++l) { 217 PetscInt dof, off, d, q, anDof; 218 PetscInt p = leaves[l], numAdj = PETSC_DETERMINE; 219 220 if ((p < pStart) || (p >= pEnd)) continue; 221 PetscCall(PetscSectionGetDof(section, p, &dof)); 222 PetscCall(PetscSectionGetOffset(section, p, &off)); 223 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 224 for (q = 0; q < numAdj; ++q) { 225 const PetscInt padj = tmpAdj[q]; 226 PetscInt ndof, ncdof; 227 228 if ((padj < pStart) || (padj >= pEnd)) continue; 229 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 230 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 231 for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(leafSectionAdj, d, ndof - ncdof)); 232 } 233 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 234 if (anDof) { 235 for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(leafSectionAdj, d, anDof)); 236 } 237 } 238 PetscCall(PetscSectionSetUp(leafSectionAdj)); 239 if (debug) { 240 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Leaves:\n")); 241 PetscCall(PetscSectionView(leafSectionAdj, NULL)); 242 } 243 /* Get maximum remote adjacency sizes for owned dofs on interface (roots) */ 244 if (doComm) { 245 PetscCall(PetscSFReduceBegin(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM)); 246 PetscCall(PetscSFReduceEnd(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM)); 247 PetscCall(PetscSectionInvalidateMaxDof_Internal(rootSectionAdj)); 248 } 249 if (debug) { 250 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots:\n")); 251 PetscCall(PetscSectionView(rootSectionAdj, NULL)); 252 } 253 /* Add in local adjacency sizes for owned dofs on interface (roots) */ 254 for (p = pStart; p < pEnd; ++p) { 255 PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof; 256 257 PetscCall(PetscSectionGetDof(section, p, &dof)); 258 PetscCall(PetscSectionGetOffset(section, p, &off)); 259 if (!dof) continue; 260 PetscCall(PetscSectionGetDof(rootSectionAdj, off, &adof)); 261 if (adof <= 0) continue; 262 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 263 for (q = 0; q < numAdj; ++q) { 264 const PetscInt padj = tmpAdj[q]; 265 PetscInt ndof, ncdof; 266 267 if ((padj < pStart) || (padj >= pEnd)) continue; 268 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 269 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 270 for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(rootSectionAdj, d, ndof - ncdof)); 271 } 272 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 273 if (anDof) { 274 for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(rootSectionAdj, d, anDof)); 275 } 276 } 277 PetscCall(PetscSectionSetUp(rootSectionAdj)); 278 if (debug) { 279 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots after local additions:\n")); 280 PetscCall(PetscSectionView(rootSectionAdj, NULL)); 281 } 282 /* Create adj SF based on dof SF */ 283 PetscCall(PetscSFCreateRemoteOffsets(sfDof, rootSectionAdj, leafSectionAdj, &remoteOffsets)); 284 PetscCall(PetscSFCreateSectionSF(sfDof, rootSectionAdj, remoteOffsets, leafSectionAdj, &sfAdj)); 285 PetscCall(PetscFree(remoteOffsets)); 286 if (debug && size > 1) { 287 PetscCall(PetscPrintf(comm, "Adjacency SF for Preallocation:\n")); 288 PetscCall(PetscSFView(sfAdj, NULL)); 289 } 290 /* Create leaf adjacency */ 291 PetscCall(PetscSectionSetUp(leafSectionAdj)); 292 PetscCall(PetscSectionGetStorageSize(leafSectionAdj, &adjSize)); 293 PetscCall(PetscCalloc1(adjSize, &adj)); 294 for (l = 0; l < nleaves; ++l) { 295 PetscInt dof, off, d, q, anDof, anOff; 296 PetscInt p = leaves[l], numAdj = PETSC_DETERMINE; 297 298 if ((p < pStart) || (p >= pEnd)) continue; 299 PetscCall(PetscSectionGetDof(section, p, &dof)); 300 PetscCall(PetscSectionGetOffset(section, p, &off)); 301 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 302 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 303 PetscCall(PetscSectionGetOffset(anchorSectionAdj, p, &anOff)); 304 for (d = off; d < off + dof; ++d) { 305 PetscInt aoff, i = 0; 306 307 PetscCall(PetscSectionGetOffset(leafSectionAdj, d, &aoff)); 308 for (q = 0; q < numAdj; ++q) { 309 const PetscInt padj = tmpAdj[q]; 310 PetscInt ndof, ncdof, ngoff, nd; 311 312 if ((padj < pStart) || (padj >= pEnd)) continue; 313 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 314 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 315 PetscCall(PetscSectionGetOffset(sectionGlobal, padj, &ngoff)); 316 for (nd = 0; nd < ndof - ncdof; ++nd) { 317 adj[aoff + i] = (ngoff < 0 ? -(ngoff + 1) : ngoff) + nd; 318 ++i; 319 } 320 } 321 for (q = 0; q < anDof; q++) { 322 adj[aoff + i] = anchorAdj[anOff + q]; 323 ++i; 324 } 325 } 326 } 327 /* Debugging */ 328 if (debug) { 329 IS tmp; 330 PetscCall(PetscPrintf(comm, "Leaf adjacency indices\n")); 331 PetscCall(ISCreateGeneral(comm, adjSize, adj, PETSC_USE_POINTER, &tmp)); 332 PetscCall(ISView(tmp, NULL)); 333 PetscCall(ISDestroy(&tmp)); 334 } 335 /* Gather adjacent indices to root */ 336 PetscCall(PetscSectionGetStorageSize(rootSectionAdj, &adjSize)); 337 PetscCall(PetscMalloc1(adjSize, &rootAdj)); 338 for (r = 0; r < adjSize; ++r) rootAdj[r] = -1; 339 if (doComm) { 340 const PetscInt *indegree; 341 PetscInt *remoteadj, radjsize = 0; 342 343 PetscCall(PetscSFComputeDegreeBegin(sfAdj, &indegree)); 344 PetscCall(PetscSFComputeDegreeEnd(sfAdj, &indegree)); 345 for (p = 0; p < adjSize; ++p) radjsize += indegree[p]; 346 PetscCall(PetscMalloc1(radjsize, &remoteadj)); 347 PetscCall(PetscSFGatherBegin(sfAdj, MPIU_INT, adj, remoteadj)); 348 PetscCall(PetscSFGatherEnd(sfAdj, MPIU_INT, adj, remoteadj)); 349 for (p = 0, l = 0, r = 0; p < adjSize; ++p, l = PetscMax(p, l + indegree[p - 1])) { 350 PetscInt s; 351 for (s = 0; s < indegree[p]; ++s, ++r) rootAdj[l + s] = remoteadj[r]; 352 } 353 PetscCheck(r == radjsize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %" PetscInt_FMT " != %" PetscInt_FMT, r, radjsize); 354 PetscCheck(l == adjSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %" PetscInt_FMT " != %" PetscInt_FMT, l, adjSize); 355 PetscCall(PetscFree(remoteadj)); 356 } 357 PetscCall(PetscSFDestroy(&sfAdj)); 358 PetscCall(PetscFree(adj)); 359 /* Debugging */ 360 if (debug) { 361 IS tmp; 362 PetscCall(PetscPrintf(comm, "Root adjacency indices after gather\n")); 363 PetscCall(ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp)); 364 PetscCall(ISView(tmp, NULL)); 365 PetscCall(ISDestroy(&tmp)); 366 } 367 /* Add in local adjacency indices for owned dofs on interface (roots) */ 368 for (p = pStart; p < pEnd; ++p) { 369 PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof, anOff; 370 371 PetscCall(PetscSectionGetDof(section, p, &dof)); 372 PetscCall(PetscSectionGetOffset(section, p, &off)); 373 if (!dof) continue; 374 PetscCall(PetscSectionGetDof(rootSectionAdj, off, &adof)); 375 if (adof <= 0) continue; 376 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 377 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 378 PetscCall(PetscSectionGetOffset(anchorSectionAdj, p, &anOff)); 379 for (d = off; d < off + dof; ++d) { 380 PetscInt adof, aoff, i; 381 382 PetscCall(PetscSectionGetDof(rootSectionAdj, d, &adof)); 383 PetscCall(PetscSectionGetOffset(rootSectionAdj, d, &aoff)); 384 i = adof - 1; 385 for (q = 0; q < anDof; q++) { 386 rootAdj[aoff + i] = anchorAdj[anOff + q]; 387 --i; 388 } 389 for (q = 0; q < numAdj; ++q) { 390 const PetscInt padj = tmpAdj[q]; 391 PetscInt ndof, ncdof, ngoff, nd; 392 393 if ((padj < pStart) || (padj >= pEnd)) continue; 394 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 395 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 396 PetscCall(PetscSectionGetOffset(sectionGlobal, padj, &ngoff)); 397 for (nd = 0; nd < ndof - ncdof; ++nd) { 398 rootAdj[aoff + i] = ngoff < 0 ? -(ngoff + 1) + nd : ngoff + nd; 399 --i; 400 } 401 } 402 } 403 } 404 /* Debugging */ 405 if (debug) { 406 IS tmp; 407 PetscCall(PetscPrintf(comm, "Root adjacency indices\n")); 408 PetscCall(ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp)); 409 PetscCall(ISView(tmp, NULL)); 410 PetscCall(ISDestroy(&tmp)); 411 } 412 /* Compress indices */ 413 PetscCall(PetscSectionSetUp(rootSectionAdj)); 414 for (p = pStart; p < pEnd; ++p) { 415 PetscInt dof, cdof, off, d; 416 PetscInt adof, aoff; 417 418 PetscCall(PetscSectionGetDof(section, p, &dof)); 419 PetscCall(PetscSectionGetConstraintDof(section, p, &cdof)); 420 PetscCall(PetscSectionGetOffset(section, p, &off)); 421 if (!dof) continue; 422 PetscCall(PetscSectionGetDof(rootSectionAdj, off, &adof)); 423 if (adof <= 0) continue; 424 for (d = off; d < off + dof - cdof; ++d) { 425 PetscCall(PetscSectionGetDof(rootSectionAdj, d, &adof)); 426 PetscCall(PetscSectionGetOffset(rootSectionAdj, d, &aoff)); 427 PetscCall(PetscSortRemoveDupsInt(&adof, &rootAdj[aoff])); 428 PetscCall(PetscSectionSetDof(rootSectionAdj, d, adof)); 429 } 430 } 431 /* Debugging */ 432 if (debug) { 433 IS tmp; 434 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots after compression:\n")); 435 PetscCall(PetscSectionView(rootSectionAdj, NULL)); 436 PetscCall(PetscPrintf(comm, "Root adjacency indices after compression\n")); 437 PetscCall(ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp)); 438 PetscCall(ISView(tmp, NULL)); 439 PetscCall(ISDestroy(&tmp)); 440 } 441 /* Build adjacency section: Maps global indices to sets of adjacent global indices */ 442 PetscCall(PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd)); 443 PetscCall(PetscSectionCreate(comm, §ionAdj)); 444 PetscCall(PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd)); 445 for (p = pStart; p < pEnd; ++p) { 446 PetscInt numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof; 447 PetscBool found = PETSC_TRUE; 448 449 PetscCall(PetscSectionGetDof(section, p, &dof)); 450 PetscCall(PetscSectionGetConstraintDof(section, p, &cdof)); 451 PetscCall(PetscSectionGetOffset(section, p, &off)); 452 PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff)); 453 for (d = 0; d < dof - cdof; ++d) { 454 PetscInt ldof, rdof; 455 456 PetscCall(PetscSectionGetDof(leafSectionAdj, off + d, &ldof)); 457 PetscCall(PetscSectionGetDof(rootSectionAdj, off + d, &rdof)); 458 if (ldof > 0) { 459 /* We do not own this point */ 460 } else if (rdof > 0) { 461 PetscCall(PetscSectionSetDof(sectionAdj, goff + d, rdof)); 462 } else { 463 found = PETSC_FALSE; 464 } 465 } 466 if (found) continue; 467 PetscCall(PetscSectionGetDof(section, p, &dof)); 468 PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff)); 469 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 470 for (q = 0; q < numAdj; ++q) { 471 const PetscInt padj = tmpAdj[q]; 472 PetscInt ndof, ncdof, noff; 473 474 if ((padj < pStart) || (padj >= pEnd)) continue; 475 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 476 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 477 PetscCall(PetscSectionGetOffset(section, padj, &noff)); 478 for (d = goff; d < goff + dof - cdof; ++d) PetscCall(PetscSectionAddDof(sectionAdj, d, ndof - ncdof)); 479 } 480 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 481 if (anDof) { 482 for (d = goff; d < goff + dof - cdof; ++d) PetscCall(PetscSectionAddDof(sectionAdj, d, anDof)); 483 } 484 } 485 PetscCall(PetscSectionSetUp(sectionAdj)); 486 if (debug) { 487 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation:\n")); 488 PetscCall(PetscSectionView(sectionAdj, NULL)); 489 } 490 /* Get adjacent indices */ 491 PetscCall(PetscSectionGetStorageSize(sectionAdj, &numCols)); 492 PetscCall(PetscMalloc1(numCols, &cols)); 493 for (p = pStart; p < pEnd; ++p) { 494 PetscInt numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof, anOff; 495 PetscBool found = PETSC_TRUE; 496 497 PetscCall(PetscSectionGetDof(section, p, &dof)); 498 PetscCall(PetscSectionGetConstraintDof(section, p, &cdof)); 499 PetscCall(PetscSectionGetOffset(section, p, &off)); 500 PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff)); 501 for (d = 0; d < dof - cdof; ++d) { 502 PetscInt ldof, rdof; 503 504 PetscCall(PetscSectionGetDof(leafSectionAdj, off + d, &ldof)); 505 PetscCall(PetscSectionGetDof(rootSectionAdj, off + d, &rdof)); 506 if (ldof > 0) { 507 /* We do not own this point */ 508 } else if (rdof > 0) { 509 PetscInt aoff, roff; 510 511 PetscCall(PetscSectionGetOffset(sectionAdj, goff + d, &aoff)); 512 PetscCall(PetscSectionGetOffset(rootSectionAdj, off + d, &roff)); 513 PetscCall(PetscArraycpy(&cols[aoff], &rootAdj[roff], rdof)); 514 } else { 515 found = PETSC_FALSE; 516 } 517 } 518 if (found) continue; 519 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 520 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 521 PetscCall(PetscSectionGetOffset(anchorSectionAdj, p, &anOff)); 522 for (d = goff; d < goff + dof - cdof; ++d) { 523 PetscInt adof, aoff, i = 0; 524 525 PetscCall(PetscSectionGetDof(sectionAdj, d, &adof)); 526 PetscCall(PetscSectionGetOffset(sectionAdj, d, &aoff)); 527 for (q = 0; q < numAdj; ++q) { 528 const PetscInt padj = tmpAdj[q]; 529 PetscInt ndof, ncdof, ngoff, nd; 530 const PetscInt *ncind; 531 532 /* Adjacent points may not be in the section chart */ 533 if ((padj < pStart) || (padj >= pEnd)) continue; 534 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 535 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 536 PetscCall(PetscSectionGetConstraintIndices(section, padj, &ncind)); 537 PetscCall(PetscSectionGetOffset(sectionGlobal, padj, &ngoff)); 538 for (nd = 0; nd < ndof - ncdof; ++nd, ++i) cols[aoff + i] = ngoff < 0 ? -(ngoff + 1) + nd : ngoff + nd; 539 } 540 for (q = 0; q < anDof; q++, i++) cols[aoff + i] = anchorAdj[anOff + q]; 541 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); 542 } 543 } 544 PetscCall(PetscSectionDestroy(&anchorSectionAdj)); 545 PetscCall(PetscSectionDestroy(&leafSectionAdj)); 546 PetscCall(PetscSectionDestroy(&rootSectionAdj)); 547 PetscCall(PetscFree(anchorAdj)); 548 PetscCall(PetscFree(rootAdj)); 549 PetscCall(PetscFree(tmpAdj)); 550 /* Debugging */ 551 if (debug) { 552 IS tmp; 553 PetscCall(PetscPrintf(comm, "Column indices\n")); 554 PetscCall(ISCreateGeneral(comm, numCols, cols, PETSC_USE_POINTER, &tmp)); 555 PetscCall(ISView(tmp, NULL)); 556 PetscCall(ISDestroy(&tmp)); 557 } 558 559 *sA = sectionAdj; 560 *colIdx = cols; 561 PetscFunctionReturn(PETSC_SUCCESS); 562 } 563 564 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[]) 565 { 566 PetscSection section; 567 PetscInt rStart, rEnd, r, pStart, pEnd, p; 568 569 PetscFunctionBegin; 570 /* This loop needs to change to a loop over points, then field dofs, which means we need to look both sections */ 571 PetscCall(PetscLayoutGetRange(rLayout, &rStart, &rEnd)); 572 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); 573 if (f >= 0 && bs == 1) { 574 PetscCall(DMGetLocalSection(dm, §ion)); 575 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 576 for (p = pStart; p < pEnd; ++p) { 577 PetscInt rS, rE; 578 579 PetscCall(DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE)); 580 for (r = rS; r < rE; ++r) { 581 PetscInt numCols, cStart, c; 582 583 PetscCall(PetscSectionGetDof(sectionAdj, r, &numCols)); 584 PetscCall(PetscSectionGetOffset(sectionAdj, r, &cStart)); 585 for (c = cStart; c < cStart + numCols; ++c) { 586 if ((cols[c] >= rStart) && (cols[c] < rEnd)) { 587 ++dnz[r - rStart]; 588 if (cols[c] >= r) ++dnzu[r - rStart]; 589 } else { 590 ++onz[r - rStart]; 591 if (cols[c] >= r) ++onzu[r - rStart]; 592 } 593 } 594 } 595 } 596 } else { 597 /* Only loop over blocks of rows */ 598 for (r = rStart / bs; r < rEnd / bs; ++r) { 599 const PetscInt row = r * bs; 600 PetscInt numCols, cStart, c; 601 602 PetscCall(PetscSectionGetDof(sectionAdj, row, &numCols)); 603 PetscCall(PetscSectionGetOffset(sectionAdj, row, &cStart)); 604 for (c = cStart; c < cStart + numCols; ++c) { 605 if ((cols[c] >= rStart) && (cols[c] < rEnd)) { 606 ++dnz[r - rStart / bs]; 607 if (cols[c] >= row) ++dnzu[r - rStart / bs]; 608 } else { 609 ++onz[r - rStart / bs]; 610 if (cols[c] >= row) ++onzu[r - rStart / bs]; 611 } 612 } 613 } 614 for (r = 0; r < (rEnd - rStart) / bs; ++r) { 615 dnz[r] /= bs; 616 onz[r] /= bs; 617 dnzu[r] /= bs; 618 onzu[r] /= bs; 619 } 620 } 621 PetscFunctionReturn(PETSC_SUCCESS); 622 } 623 624 static PetscErrorCode DMPlexFillMatrix_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], Mat A) 625 { 626 PetscSection section; 627 PetscScalar *values; 628 PetscInt rStart, rEnd, r, pStart, pEnd, p, len, maxRowLen = 0; 629 630 PetscFunctionBegin; 631 PetscCall(PetscLayoutGetRange(rLayout, &rStart, &rEnd)); 632 for (r = rStart; r < rEnd; ++r) { 633 PetscCall(PetscSectionGetDof(sectionAdj, r, &len)); 634 maxRowLen = PetscMax(maxRowLen, len); 635 } 636 PetscCall(PetscCalloc1(maxRowLen, &values)); 637 if (f >= 0 && bs == 1) { 638 PetscCall(DMGetLocalSection(dm, §ion)); 639 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 640 for (p = pStart; p < pEnd; ++p) { 641 PetscInt rS, rE; 642 643 PetscCall(DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE)); 644 for (r = rS; r < rE; ++r) { 645 PetscInt numCols, cStart; 646 647 PetscCall(PetscSectionGetDof(sectionAdj, r, &numCols)); 648 PetscCall(PetscSectionGetOffset(sectionAdj, r, &cStart)); 649 PetscCall(MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES)); 650 } 651 } 652 } else { 653 for (r = rStart; r < rEnd; ++r) { 654 PetscInt numCols, cStart; 655 656 PetscCall(PetscSectionGetDof(sectionAdj, r, &numCols)); 657 PetscCall(PetscSectionGetOffset(sectionAdj, r, &cStart)); 658 PetscCall(MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES)); 659 } 660 } 661 PetscCall(PetscFree(values)); 662 PetscFunctionReturn(PETSC_SUCCESS); 663 } 664 665 /*@C 666 DMPlexPreallocateOperator - Calculate the matrix nonzero pattern based upon the information in the `DM`, 667 the `PetscDS` it contains, and the default `PetscSection`. 668 669 Collective 670 671 Input Parameters: 672 + dm - The `DMPLEX` 673 . bs - The matrix blocksize 674 . dnz - An array to hold the number of nonzeros in the diagonal block 675 . onz - An array to hold the number of nonzeros in the off-diagonal block 676 . dnzu - An array to hold the number of nonzeros in the upper triangle of the diagonal block 677 . onzu - An array to hold the number of nonzeros in the upper triangle of the off-diagonal block 678 - fillMatrix - If `PETSC_TRUE`, fill the matrix with zeros 679 680 Output Parameter: 681 . A - The preallocated matrix 682 683 Level: advanced 684 685 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreateMatrix()` 686 @*/ 687 PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix) 688 { 689 MPI_Comm comm; 690 MatType mtype; 691 PetscSF sf, sfDof; 692 PetscSection section; 693 PetscInt *remoteOffsets; 694 PetscSection sectionAdj[4] = {NULL, NULL, NULL, NULL}; 695 PetscInt *cols[4] = {NULL, NULL, NULL, NULL}; 696 PetscBool useCone, useClosure; 697 PetscInt Nf, f, idx, locRows; 698 PetscLayout rLayout; 699 PetscBool isSymBlock, isSymSeqBlock, isSymMPIBlock, debug = PETSC_FALSE; 700 PetscMPIInt size; 701 702 PetscFunctionBegin; 703 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 704 PetscValidHeaderSpecific(A, MAT_CLASSID, 7); 705 if (dnz) PetscAssertPointer(dnz, 3); 706 if (onz) PetscAssertPointer(onz, 4); 707 if (dnzu) PetscAssertPointer(dnzu, 5); 708 if (onzu) PetscAssertPointer(onzu, 6); 709 PetscCall(DMGetPointSF(dm, &sf)); 710 PetscCall(DMGetLocalSection(dm, §ion)); 711 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_view_preallocation", &debug, NULL)); 712 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 713 PetscCallMPI(MPI_Comm_size(comm, &size)); 714 PetscCall(PetscLogEventBegin(DMPLEX_Preallocate, dm, 0, 0, 0)); 715 /* Create dof SF based on point SF */ 716 if (debug) { 717 PetscSection section, sectionGlobal; 718 PetscSF sf; 719 720 PetscCall(DMGetPointSF(dm, &sf)); 721 PetscCall(DMGetLocalSection(dm, §ion)); 722 PetscCall(DMGetGlobalSection(dm, §ionGlobal)); 723 PetscCall(PetscPrintf(comm, "Input Section for Preallocation:\n")); 724 PetscCall(PetscSectionView(section, NULL)); 725 PetscCall(PetscPrintf(comm, "Input Global Section for Preallocation:\n")); 726 PetscCall(PetscSectionView(sectionGlobal, NULL)); 727 if (size > 1) { 728 PetscCall(PetscPrintf(comm, "Input SF for Preallocation:\n")); 729 PetscCall(PetscSFView(sf, NULL)); 730 } 731 } 732 PetscCall(PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets)); 733 PetscCall(PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof)); 734 PetscCall(PetscFree(remoteOffsets)); 735 if (debug && size > 1) { 736 PetscCall(PetscPrintf(comm, "Dof SF for Preallocation:\n")); 737 PetscCall(PetscSFView(sfDof, NULL)); 738 } 739 /* Create allocation vectors from adjacency graph */ 740 PetscCall(MatGetLocalSize(A, &locRows, NULL)); 741 PetscCall(PetscLayoutCreate(comm, &rLayout)); 742 PetscCall(PetscLayoutSetLocalSize(rLayout, locRows)); 743 PetscCall(PetscLayoutSetBlockSize(rLayout, 1)); 744 PetscCall(PetscLayoutSetUp(rLayout)); 745 /* There are 4 types of adjacency */ 746 PetscCall(PetscSectionGetNumFields(section, &Nf)); 747 if (Nf < 1 || bs > 1) { 748 PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure)); 749 idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0); 750 PetscCall(DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, §ionAdj[idx], &cols[idx])); 751 PetscCall(DMPlexUpdateAllocation_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu)); 752 } else { 753 for (f = 0; f < Nf; ++f) { 754 PetscCall(DMGetAdjacency(dm, f, &useCone, &useClosure)); 755 idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0); 756 if (!sectionAdj[idx]) PetscCall(DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, §ionAdj[idx], &cols[idx])); 757 PetscCall(DMPlexUpdateAllocation_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu)); 758 } 759 } 760 PetscCall(PetscSFDestroy(&sfDof)); 761 /* Set matrix pattern */ 762 PetscCall(MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu)); 763 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 764 /* Check for symmetric storage */ 765 PetscCall(MatGetType(A, &mtype)); 766 PetscCall(PetscStrcmp(mtype, MATSBAIJ, &isSymBlock)); 767 PetscCall(PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock)); 768 PetscCall(PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock)); 769 if (isSymBlock || isSymSeqBlock || isSymMPIBlock) PetscCall(MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE)); 770 /* Fill matrix with zeros */ 771 if (fillMatrix) { 772 if (Nf < 1 || bs > 1) { 773 PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure)); 774 idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0); 775 PetscCall(DMPlexFillMatrix_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], A)); 776 } else { 777 for (f = 0; f < Nf; ++f) { 778 PetscCall(DMGetAdjacency(dm, f, &useCone, &useClosure)); 779 idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0); 780 PetscCall(DMPlexFillMatrix_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], A)); 781 } 782 } 783 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 784 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 785 } 786 PetscCall(PetscLayoutDestroy(&rLayout)); 787 for (idx = 0; idx < 4; ++idx) { 788 PetscCall(PetscSectionDestroy(§ionAdj[idx])); 789 PetscCall(PetscFree(cols[idx])); 790 } 791 PetscCall(PetscLogEventEnd(DMPLEX_Preallocate, dm, 0, 0, 0)); 792 PetscFunctionReturn(PETSC_SUCCESS); 793 } 794 795 #if 0 796 PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix) 797 { 798 PetscInt *tmpClosure,*tmpAdj,*visits; 799 PetscInt c,cStart,cEnd,pStart,pEnd; 800 801 PetscFunctionBegin; 802 PetscCall(DMGetDimension(dm, &dim)); 803 PetscCall(DMPlexGetDepth(dm, &depth)); 804 PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize)); 805 806 maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1)); 807 808 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 809 npoints = pEnd - pStart; 810 811 PetscCall(PetscMalloc3(maxClosureSize,&tmpClosure,npoints,&lvisits,npoints,&visits)); 812 PetscCall(PetscArrayzero(lvisits,pEnd-pStart)); 813 PetscCall(PetscArrayzero(visits,pEnd-pStart)); 814 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 815 for (c=cStart; c<cEnd; c++) { 816 PetscInt *support = tmpClosure; 817 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support)); 818 for (p=0; p<supportSize; p++) lvisits[support[p]]++; 819 } 820 PetscCall(PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM)); 821 PetscCall(PetscSFReduceEnd (sf,MPIU_INT,lvisits,visits,MPI_SUM)); 822 PetscCall(PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits,MPI_REPLACE)); 823 PetscCall(PetscSFBcastEnd (sf,MPIU_INT,visits,lvisits)); 824 825 PetscCall(PetscSFGetRootRanks()); 826 827 PetscCall(PetscMalloc2(maxClosureSize*maxClosureSize,&cellmat,npoints,&owner)); 828 for (c=cStart; c<cEnd; c++) { 829 PetscCall(PetscArrayzero(cellmat,maxClosureSize*maxClosureSize)); 830 /* 831 Depth-first walk of transitive closure. 832 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. 833 This contribution is added to dnz if owning ranks of p and q match, to onz otherwise. 834 */ 835 } 836 837 PetscCall(PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM)); 838 PetscCall(PetscSFReduceEnd (sf,MPIU_INT,lonz,onz,MPI_SUM)); 839 PetscFunctionReturn(PETSC_SUCCESS); 840 } 841 #endif 842