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