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