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