Actual source code: plexpreallocate.c

  1: #include <petsc/private/dmpleximpl.h>
  2: #include <petsc/private/isimpl.h>
  3: #include <petscsf.h>
  4: #include <petscds.h>

  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;

 13:   PetscFunctionBegin;
 14:   PetscCall(DMGetLocalSection(dm, &section));
 15:   PetscCall(DMGetGlobalSection(dm, &sectionGlobal));
 16:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section), &adjSec));
 17:   PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
 18:   PetscCall(PetscSectionSetChart(adjSec, pStart, pEnd));

 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;

 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));

 33:     for (p = 0; p < aSize; p++) {
 34:       PetscInt a = anchors[p];

 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;

 46:       PetscCall(PetscSectionGetDof(aSec, p, &dof));
 47:       PetscCall(PetscSectionGetOffset(aSec, p, &off));

 49:       for (q = 0; q < dof; q++) {
 50:         PetscInt iOff;

 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));

 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;

 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;

 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:     }

104:     PetscCall(PetscSectionSetUp(adjSec));
105:     PetscCall(PetscSectionGetStorageSize(adjSec, &adjSize));
106:     PetscCall(PetscMalloc1(adjSize, &adj));

108:     for (p = pStart; p < pEnd; p++) {
109:       PetscInt iDof, iOff, i, r, s, aOff, aOffOrig, aDof, numAdjP = PETSC_DETERMINE;

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;

121:         q = inverse[iOff + i];
122:         PetscCall(DMPlexGetAdjacency_Internal(dm, q, useCone, useClosure, PETSC_TRUE, &numAdjQ, &tmpAdjQ));
123:         for (r = 0; r < numAdjQ; r++) {
124:           qAdj = tmpAdjQ[r];
125:           if ((qAdj < pStart) || (qAdj >= pEnd)) continue;
126:           for (s = 0; s < numAdjP; s++) {
127:             if (qAdj == tmpAdjP[s]) break;
128:           }
129:           if (s < numAdjP) continue;
130:           PetscCall(PetscSectionGetDof(section, qAdj, &qAdjDof));
131:           PetscCall(PetscSectionGetConstraintDof(section, qAdj, &qAdjCDof));
132:           PetscCall(PetscSectionGetOffset(sectionGlobal, qAdj, &qAdjOff));
133:           for (nd = 0; nd < qAdjDof - qAdjCDof; ++nd) adj[aOff++] = (qAdjOff < 0 ? -(qAdjOff + 1) : qAdjOff) + nd;
134:         }
135:       }
136:       PetscCall(PetscSortRemoveDupsInt(&aDof, PetscSafePointerPlusOffset(adj, aOffOrig)));
137:       PetscCall(PetscSectionSetDof(adjSec, p, aDof));
138:     }
139:     *anchorAdj = adj;

141:     /* clean up */
142:     PetscCall(PetscSectionDestroy(&inverseSec));
143:     PetscCall(PetscFree(inverse));
144:     PetscCall(PetscFree(tmpAdjP));
145:     PetscCall(PetscFree(tmpAdjQ));
146:   } else {
147:     *anchorAdj = NULL;
148:     PetscCall(PetscSectionSetUp(adjSec));
149:   }
150:   *anchorSectionAdj = adjSec;
151:   PetscFunctionReturn(PETSC_SUCCESS);
152: }

154: static PetscErrorCode DMPlexCreateAdjacencySection_Static(DM dm, PetscInt bs, PetscSF sfDof, PetscBool useCone, PetscBool useClosure, PetscBool useAnchors, PetscSection *sA, PetscInt **colIdx)
155: {
156:   MPI_Comm           comm;
157:   PetscMPIInt        size;
158:   PetscBool          doCommLocal, doComm, debug = PETSC_FALSE;
159:   PetscSF            sf, sfAdj;
160:   PetscSection       section, sectionGlobal, leafSectionAdj, rootSectionAdj, sectionAdj, anchorSectionAdj;
161:   PetscInt           nroots, nleaves, l, p, r;
162:   const PetscInt    *leaves;
163:   const PetscSFNode *remotes;
164:   PetscInt           dim, pStart, pEnd, numDof, globalOffStart, globalOffEnd, numCols;
165:   PetscInt          *tmpAdj = NULL, *adj, *rootAdj, *anchorAdj = NULL, *cols, *remoteOffsets;
166:   PetscInt           adjSize;

168:   PetscFunctionBegin;
169:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
170:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_view_preallocation", &debug, NULL));
171:   PetscCallMPI(MPI_Comm_size(comm, &size));
172:   PetscCall(DMGetDimension(dm, &dim));
173:   PetscCall(DMGetPointSF(dm, &sf));
174:   PetscCall(DMGetLocalSection(dm, &section));
175:   PetscCall(DMGetGlobalSection(dm, &sectionGlobal));
176:   PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
177:   doCommLocal = (size > 1) && (nroots >= 0) ? PETSC_TRUE : PETSC_FALSE;
178:   PetscCallMPI(MPIU_Allreduce(&doCommLocal, &doComm, 1, MPIU_BOOL, MPI_LAND, comm));
179:   /* Create section for dof adjacency (dof ==> # adj dof) */
180:   PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
181:   PetscCall(PetscSectionGetStorageSize(section, &numDof));
182:   PetscCall(PetscSectionCreate(comm, &leafSectionAdj));
183:   PetscCall(PetscSectionSetChart(leafSectionAdj, 0, numDof));
184:   PetscCall(PetscSectionCreate(comm, &rootSectionAdj));
185:   PetscCall(PetscSectionSetChart(rootSectionAdj, 0, numDof));
186:   /*   Fill in the ghost dofs on the interface */
187:   PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, &remotes));
188:   /*
189:    section        - maps points to (# dofs, local dofs)
190:    sectionGlobal  - maps points to (# dofs, global dofs)
191:    leafSectionAdj - maps unowned local dofs to # adj dofs
192:    rootSectionAdj - maps   owned local dofs to # adj dofs
193:    adj            - adj global dofs indexed by leafSectionAdj
194:    rootAdj        - adj global dofs indexed by rootSectionAdj
195:    sf    - describes shared points across procs
196:    sfDof - describes shared dofs across procs
197:    sfAdj - describes shared adjacent dofs across procs
198:    ** The bootstrapping process involves six rounds with similar structure of visiting neighbors of each point.
199:   (0). If there are point-to-point constraints, add the adjacencies of constrained points to anchors in anchorAdj
200:        (This is done in DMPlexComputeAnchorAdjacencies())
201:     1. Visit unowned points on interface, count adjacencies placing in leafSectionAdj
202:        Reduce those counts to rootSectionAdj (now redundantly counting some interface points)
203:     2. Visit owned points on interface, count adjacencies placing in rootSectionAdj
204:        Create sfAdj connecting rootSectionAdj and leafSectionAdj
205:     3. Visit unowned points on interface, write adjacencies to adj
206:        Gather adj to rootAdj (note that there is redundancy in rootAdj when multiple procs find the same adjacencies)
207:     4. Visit owned points on interface, write adjacencies to rootAdj
208:        Remove redundancy in rootAdj
209:    ** The last two traversals use transitive closure
210:     5. Visit all owned points in the subdomain, count dofs for each point (sectionAdj)
211:        Allocate memory addressed by sectionAdj (cols)
212:     6. Visit all owned points in the subdomain, insert dof adjacencies into cols
213:    ** Knowing all the column adjacencies, check ownership and sum into dnz and onz
214:   */
215:   PetscCall(DMPlexComputeAnchorAdjacencies(dm, useCone, useClosure, &anchorSectionAdj, &anchorAdj));
216:   for (l = 0; l < nleaves; ++l) {
217:     PetscInt dof, off, d, q, anDof;
218:     PetscInt p = leaves[l], numAdj = PETSC_DETERMINE;

220:     if ((p < pStart) || (p >= pEnd)) continue;
221:     PetscCall(PetscSectionGetDof(section, p, &dof));
222:     PetscCall(PetscSectionGetOffset(section, p, &off));
223:     PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj));
224:     for (q = 0; q < numAdj; ++q) {
225:       const PetscInt padj = tmpAdj[q];
226:       PetscInt       ndof, ncdof;

228:       if ((padj < pStart) || (padj >= pEnd)) continue;
229:       PetscCall(PetscSectionGetDof(section, padj, &ndof));
230:       PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof));
231:       for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(leafSectionAdj, d, ndof - ncdof));
232:     }
233:     PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof));
234:     if (anDof) {
235:       for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(leafSectionAdj, d, anDof));
236:     }
237:   }
238:   PetscCall(PetscSectionSetUp(leafSectionAdj));
239:   if (debug) {
240:     PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Leaves:\n"));
241:     PetscCall(PetscSectionView(leafSectionAdj, NULL));
242:   }
243:   /* Get maximum remote adjacency sizes for owned dofs on interface (roots) */
244:   if (doComm) {
245:     PetscCall(PetscSFReduceBegin(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM));
246:     PetscCall(PetscSFReduceEnd(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM));
247:     PetscCall(PetscSectionInvalidateMaxDof_Internal(rootSectionAdj));
248:   }
249:   if (debug) {
250:     PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots:\n"));
251:     PetscCall(PetscSectionView(rootSectionAdj, NULL));
252:   }
253:   /* Add in local adjacency sizes for owned dofs on interface (roots) */
254:   for (p = pStart; p < pEnd; ++p) {
255:     PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof;

257:     PetscCall(PetscSectionGetDof(section, p, &dof));
258:     PetscCall(PetscSectionGetOffset(section, p, &off));
259:     if (!dof) continue;
260:     PetscCall(PetscSectionGetDof(rootSectionAdj, off, &adof));
261:     if (adof <= 0) continue;
262:     PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj));
263:     for (q = 0; q < numAdj; ++q) {
264:       const PetscInt padj = tmpAdj[q];
265:       PetscInt       ndof, ncdof;

267:       if ((padj < pStart) || (padj >= pEnd)) continue;
268:       PetscCall(PetscSectionGetDof(section, padj, &ndof));
269:       PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof));
270:       for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(rootSectionAdj, d, ndof - ncdof));
271:     }
272:     PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof));
273:     if (anDof) {
274:       for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(rootSectionAdj, d, anDof));
275:     }
276:   }
277:   PetscCall(PetscSectionSetUp(rootSectionAdj));
278:   if (debug) {
279:     PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots after local additions:\n"));
280:     PetscCall(PetscSectionView(rootSectionAdj, NULL));
281:   }
282:   /* Create adj SF based on dof SF */
283:   PetscCall(PetscSFCreateRemoteOffsets(sfDof, rootSectionAdj, leafSectionAdj, &remoteOffsets));
284:   PetscCall(PetscSFCreateSectionSF(sfDof, rootSectionAdj, remoteOffsets, leafSectionAdj, &sfAdj));
285:   PetscCall(PetscFree(remoteOffsets));
286:   if (debug && size > 1) {
287:     PetscCall(PetscPrintf(comm, "Adjacency SF for Preallocation:\n"));
288:     PetscCall(PetscSFView(sfAdj, NULL));
289:   }
290:   /* Create leaf adjacency */
291:   PetscCall(PetscSectionSetUp(leafSectionAdj));
292:   PetscCall(PetscSectionGetStorageSize(leafSectionAdj, &adjSize));
293:   PetscCall(PetscCalloc1(adjSize, &adj));
294:   for (l = 0; l < nleaves; ++l) {
295:     PetscInt dof, off, d, q, anDof, anOff;
296:     PetscInt p = leaves[l], numAdj = PETSC_DETERMINE;

298:     if ((p < pStart) || (p >= pEnd)) continue;
299:     PetscCall(PetscSectionGetDof(section, p, &dof));
300:     PetscCall(PetscSectionGetOffset(section, p, &off));
301:     PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj));
302:     PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof));
303:     PetscCall(PetscSectionGetOffset(anchorSectionAdj, p, &anOff));
304:     for (d = off; d < off + dof; ++d) {
305:       PetscInt aoff, i = 0;

307:       PetscCall(PetscSectionGetOffset(leafSectionAdj, d, &aoff));
308:       for (q = 0; q < numAdj; ++q) {
309:         const PetscInt padj = tmpAdj[q];
310:         PetscInt       ndof, ncdof, ngoff, nd;

312:         if ((padj < pStart) || (padj >= pEnd)) continue;
313:         PetscCall(PetscSectionGetDof(section, padj, &ndof));
314:         PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof));
315:         PetscCall(PetscSectionGetOffset(sectionGlobal, padj, &ngoff));
316:         for (nd = 0; nd < ndof - ncdof; ++nd) {
317:           adj[aoff + i] = (ngoff < 0 ? -(ngoff + 1) : ngoff) + nd;
318:           ++i;
319:         }
320:       }
321:       for (q = 0; q < anDof; q++) {
322:         adj[aoff + i] = anchorAdj[anOff + q];
323:         ++i;
324:       }
325:     }
326:   }
327:   /* Debugging */
328:   if (debug) {
329:     IS tmp;
330:     PetscCall(PetscPrintf(comm, "Leaf adjacency indices\n"));
331:     PetscCall(ISCreateGeneral(comm, adjSize, adj, PETSC_USE_POINTER, &tmp));
332:     PetscCall(ISView(tmp, NULL));
333:     PetscCall(ISDestroy(&tmp));
334:   }
335:   /* Gather adjacent indices to root */
336:   PetscCall(PetscSectionGetStorageSize(rootSectionAdj, &adjSize));
337:   PetscCall(PetscMalloc1(adjSize, &rootAdj));
338:   for (r = 0; r < adjSize; ++r) rootAdj[r] = -1;
339:   if (doComm) {
340:     const PetscInt *indegree;
341:     PetscInt       *remoteadj, radjsize = 0;

343:     PetscCall(PetscSFComputeDegreeBegin(sfAdj, &indegree));
344:     PetscCall(PetscSFComputeDegreeEnd(sfAdj, &indegree));
345:     for (p = 0; p < adjSize; ++p) radjsize += indegree[p];
346:     PetscCall(PetscMalloc1(radjsize, &remoteadj));
347:     PetscCall(PetscSFGatherBegin(sfAdj, MPIU_INT, adj, remoteadj));
348:     PetscCall(PetscSFGatherEnd(sfAdj, MPIU_INT, adj, remoteadj));
349:     for (p = 0, l = 0, r = 0; p < adjSize; ++p, l = PetscMax(p, l + indegree[p - 1])) {
350:       PetscInt s;
351:       for (s = 0; s < indegree[p]; ++s, ++r) rootAdj[l + s] = remoteadj[r];
352:     }
353:     PetscCheck(r == radjsize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %" PetscInt_FMT " != %" PetscInt_FMT, r, radjsize);
354:     PetscCheck(l == adjSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %" PetscInt_FMT " != %" PetscInt_FMT, l, adjSize);
355:     PetscCall(PetscFree(remoteadj));
356:   }
357:   PetscCall(PetscSFDestroy(&sfAdj));
358:   PetscCall(PetscFree(adj));
359:   /* Debugging */
360:   if (debug) {
361:     IS tmp;
362:     PetscCall(PetscPrintf(comm, "Root adjacency indices after gather\n"));
363:     PetscCall(ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp));
364:     PetscCall(ISView(tmp, NULL));
365:     PetscCall(ISDestroy(&tmp));
366:   }
367:   /* Add in local adjacency indices for owned dofs on interface (roots) */
368:   for (p = pStart; p < pEnd; ++p) {
369:     PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof, anOff;

371:     PetscCall(PetscSectionGetDof(section, p, &dof));
372:     PetscCall(PetscSectionGetOffset(section, p, &off));
373:     if (!dof) continue;
374:     PetscCall(PetscSectionGetDof(rootSectionAdj, off, &adof));
375:     if (adof <= 0) continue;
376:     PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj));
377:     PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof));
378:     PetscCall(PetscSectionGetOffset(anchorSectionAdj, p, &anOff));
379:     for (d = off; d < off + dof; ++d) {
380:       PetscInt adof, aoff, i;

382:       PetscCall(PetscSectionGetDof(rootSectionAdj, d, &adof));
383:       PetscCall(PetscSectionGetOffset(rootSectionAdj, d, &aoff));
384:       i = adof - 1;
385:       for (q = 0; q < anDof; q++) {
386:         rootAdj[aoff + i] = anchorAdj[anOff + q];
387:         --i;
388:       }
389:       for (q = 0; q < numAdj; ++q) {
390:         const PetscInt padj = tmpAdj[q];
391:         PetscInt       ndof, ncdof, ngoff, nd;

393:         if ((padj < pStart) || (padj >= pEnd)) continue;
394:         PetscCall(PetscSectionGetDof(section, padj, &ndof));
395:         PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof));
396:         PetscCall(PetscSectionGetOffset(sectionGlobal, padj, &ngoff));
397:         for (nd = 0; nd < ndof - ncdof; ++nd) {
398:           rootAdj[aoff + i] = ngoff < 0 ? -(ngoff + 1) + nd : ngoff + nd;
399:           --i;
400:         }
401:       }
402:     }
403:   }
404:   /* Debugging */
405:   if (debug) {
406:     IS tmp;
407:     PetscCall(PetscPrintf(comm, "Root adjacency indices\n"));
408:     PetscCall(ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp));
409:     PetscCall(ISView(tmp, NULL));
410:     PetscCall(ISDestroy(&tmp));
411:   }
412:   /* Compress indices */
413:   PetscCall(PetscSectionSetUp(rootSectionAdj));
414:   for (p = pStart; p < pEnd; ++p) {
415:     PetscInt dof, cdof, off, d;
416:     PetscInt adof, aoff;

418:     PetscCall(PetscSectionGetDof(section, p, &dof));
419:     PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
420:     PetscCall(PetscSectionGetOffset(section, p, &off));
421:     if (!dof) continue;
422:     PetscCall(PetscSectionGetDof(rootSectionAdj, off, &adof));
423:     if (adof <= 0) continue;
424:     for (d = off; d < off + dof - cdof; ++d) {
425:       PetscCall(PetscSectionGetDof(rootSectionAdj, d, &adof));
426:       PetscCall(PetscSectionGetOffset(rootSectionAdj, d, &aoff));
427:       PetscCall(PetscSortRemoveDupsInt(&adof, &rootAdj[aoff]));
428:       PetscCall(PetscSectionSetDof(rootSectionAdj, d, adof));
429:     }
430:   }
431:   /* Debugging */
432:   if (debug) {
433:     IS tmp;
434:     PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots after compression:\n"));
435:     PetscCall(PetscSectionView(rootSectionAdj, NULL));
436:     PetscCall(PetscPrintf(comm, "Root adjacency indices after compression\n"));
437:     PetscCall(ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp));
438:     PetscCall(ISView(tmp, NULL));
439:     PetscCall(ISDestroy(&tmp));
440:   }
441:   /* Build adjacency section: Maps global indices to sets of adjacent global indices */
442:   PetscCall(PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd));
443:   PetscCall(PetscSectionCreate(comm, &sectionAdj));
444:   PetscCall(PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd));
445:   for (p = pStart; p < pEnd; ++p) {
446:     PetscInt  numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof;
447:     PetscBool found  = PETSC_TRUE;

449:     PetscCall(PetscSectionGetDof(section, p, &dof));
450:     PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
451:     PetscCall(PetscSectionGetOffset(section, p, &off));
452:     PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff));
453:     for (d = 0; d < dof - cdof; ++d) {
454:       PetscInt ldof, rdof;

456:       PetscCall(PetscSectionGetDof(leafSectionAdj, off + d, &ldof));
457:       PetscCall(PetscSectionGetDof(rootSectionAdj, off + d, &rdof));
458:       if (ldof > 0) {
459:         /* We do not own this point */
460:       } else if (rdof > 0) {
461:         PetscCall(PetscSectionSetDof(sectionAdj, goff + d, rdof));
462:       } else {
463:         found = PETSC_FALSE;
464:       }
465:     }
466:     if (found) continue;
467:     PetscCall(PetscSectionGetDof(section, p, &dof));
468:     PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff));
469:     PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj));
470:     for (q = 0; q < numAdj; ++q) {
471:       const PetscInt padj = tmpAdj[q];
472:       PetscInt       ndof, ncdof, noff;

474:       if ((padj < pStart) || (padj >= pEnd)) continue;
475:       PetscCall(PetscSectionGetDof(section, padj, &ndof));
476:       PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof));
477:       PetscCall(PetscSectionGetOffset(section, padj, &noff));
478:       for (d = goff; d < goff + dof - cdof; ++d) PetscCall(PetscSectionAddDof(sectionAdj, d, ndof - ncdof));
479:     }
480:     PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof));
481:     if (anDof) {
482:       for (d = goff; d < goff + dof - cdof; ++d) PetscCall(PetscSectionAddDof(sectionAdj, d, anDof));
483:     }
484:   }
485:   PetscCall(PetscSectionSetUp(sectionAdj));
486:   if (debug) {
487:     PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation:\n"));
488:     PetscCall(PetscSectionView(sectionAdj, NULL));
489:   }
490:   /* Get adjacent indices */
491:   PetscCall(PetscSectionGetStorageSize(sectionAdj, &numCols));
492:   PetscCall(PetscMalloc1(numCols, &cols));
493:   for (p = pStart; p < pEnd; ++p) {
494:     PetscInt  numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof, anOff;
495:     PetscBool found  = PETSC_TRUE;

497:     PetscCall(PetscSectionGetDof(section, p, &dof));
498:     PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
499:     PetscCall(PetscSectionGetOffset(section, p, &off));
500:     PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff));
501:     for (d = 0; d < dof - cdof; ++d) {
502:       PetscInt ldof, rdof;

504:       PetscCall(PetscSectionGetDof(leafSectionAdj, off + d, &ldof));
505:       PetscCall(PetscSectionGetDof(rootSectionAdj, off + d, &rdof));
506:       if (ldof > 0) {
507:         /* We do not own this point */
508:       } else if (rdof > 0) {
509:         PetscInt aoff, roff;

511:         PetscCall(PetscSectionGetOffset(sectionAdj, goff + d, &aoff));
512:         PetscCall(PetscSectionGetOffset(rootSectionAdj, off + d, &roff));
513:         PetscCall(PetscArraycpy(&cols[aoff], &rootAdj[roff], rdof));
514:       } else {
515:         found = PETSC_FALSE;
516:       }
517:     }
518:     if (found) continue;
519:     PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj));
520:     PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof));
521:     PetscCall(PetscSectionGetOffset(anchorSectionAdj, p, &anOff));
522:     for (d = goff; d < goff + dof - cdof; ++d) {
523:       PetscInt adof, aoff, i = 0;

525:       PetscCall(PetscSectionGetDof(sectionAdj, d, &adof));
526:       PetscCall(PetscSectionGetOffset(sectionAdj, d, &aoff));
527:       for (q = 0; q < numAdj; ++q) {
528:         const PetscInt  padj = tmpAdj[q];
529:         PetscInt        ndof, ncdof, ngoff, nd;
530:         const PetscInt *ncind;

532:         /* Adjacent points may not be in the section chart */
533:         if ((padj < pStart) || (padj >= pEnd)) continue;
534:         PetscCall(PetscSectionGetDof(section, padj, &ndof));
535:         PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof));
536:         PetscCall(PetscSectionGetConstraintIndices(section, padj, &ncind));
537:         PetscCall(PetscSectionGetOffset(sectionGlobal, padj, &ngoff));
538:         for (nd = 0; nd < ndof - ncdof; ++nd, ++i) cols[aoff + i] = ngoff < 0 ? -(ngoff + 1) + nd : ngoff + nd;
539:       }
540:       for (q = 0; q < anDof; q++, i++) cols[aoff + i] = anchorAdj[anOff + q];
541:       PetscCheck(i == adof, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of entries %" PetscInt_FMT " != %" PetscInt_FMT " for dof %" PetscInt_FMT " (point %" PetscInt_FMT ")", i, adof, d, p);
542:     }
543:   }
544:   PetscCall(PetscSectionDestroy(&anchorSectionAdj));
545:   PetscCall(PetscSectionDestroy(&leafSectionAdj));
546:   PetscCall(PetscSectionDestroy(&rootSectionAdj));
547:   PetscCall(PetscFree(anchorAdj));
548:   PetscCall(PetscFree(rootAdj));
549:   PetscCall(PetscFree(tmpAdj));
550:   /* Debugging */
551:   if (debug) {
552:     IS tmp;
553:     PetscCall(PetscPrintf(comm, "Column indices\n"));
554:     PetscCall(ISCreateGeneral(comm, numCols, cols, PETSC_USE_POINTER, &tmp));
555:     PetscCall(ISView(tmp, NULL));
556:     PetscCall(ISDestroy(&tmp));
557:   }

559:   *sA     = sectionAdj;
560:   *colIdx = cols;
561:   PetscFunctionReturn(PETSC_SUCCESS);
562: }

564: static PetscErrorCode DMPlexUpdateAllocation_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[])
565: {
566:   PetscSection section;
567:   PetscInt     rStart, rEnd, r, pStart, pEnd, p;

569:   PetscFunctionBegin;
570:   /* This loop needs to change to a loop over points, then field dofs, which means we need to look both sections */
571:   PetscCall(PetscLayoutGetRange(rLayout, &rStart, &rEnd));
572:   PetscCheck((rStart % bs) == 0 && (rEnd % bs) == 0, PetscObjectComm((PetscObject)rLayout), PETSC_ERR_ARG_WRONG, "Invalid layout [%" PetscInt_FMT ", %" PetscInt_FMT ") for matrix, must be divisible by block size %" PetscInt_FMT, rStart, rEnd, bs);
573:   if (f >= 0 && bs == 1) {
574:     PetscCall(DMGetLocalSection(dm, &section));
575:     PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
576:     for (p = pStart; p < pEnd; ++p) {
577:       PetscInt rS, rE;

579:       PetscCall(DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE));
580:       for (r = rS; r < rE; ++r) {
581:         PetscInt numCols, cStart, c;

583:         PetscCall(PetscSectionGetDof(sectionAdj, r, &numCols));
584:         PetscCall(PetscSectionGetOffset(sectionAdj, r, &cStart));
585:         for (c = cStart; c < cStart + numCols; ++c) {
586:           if ((cols[c] >= rStart) && (cols[c] < rEnd)) {
587:             ++dnz[r - rStart];
588:             if (cols[c] >= r) ++dnzu[r - rStart];
589:           } else {
590:             ++onz[r - rStart];
591:             if (cols[c] >= r) ++onzu[r - rStart];
592:           }
593:         }
594:       }
595:     }
596:   } else {
597:     /* Only loop over blocks of rows */
598:     for (r = rStart / bs; r < rEnd / bs; ++r) {
599:       const PetscInt row = r * bs;
600:       PetscInt       numCols, cStart, c;

602:       PetscCall(PetscSectionGetDof(sectionAdj, row, &numCols));
603:       PetscCall(PetscSectionGetOffset(sectionAdj, row, &cStart));
604:       for (c = cStart; c < cStart + numCols; ++c) {
605:         if ((cols[c] >= rStart) && (cols[c] < rEnd)) {
606:           ++dnz[r - rStart / bs];
607:           if (cols[c] >= row) ++dnzu[r - rStart / bs];
608:         } else {
609:           ++onz[r - rStart / bs];
610:           if (cols[c] >= row) ++onzu[r - rStart / bs];
611:         }
612:       }
613:     }
614:     for (r = 0; r < (rEnd - rStart) / bs; ++r) {
615:       dnz[r] /= bs;
616:       onz[r] /= bs;
617:       dnzu[r] /= bs;
618:       onzu[r] /= bs;
619:     }
620:   }
621:   PetscFunctionReturn(PETSC_SUCCESS);
622: }

624: static PetscErrorCode DMPlexFillMatrix_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], Mat A)
625: {
626:   PetscSection section;
627:   PetscScalar *values;
628:   PetscInt     rStart, rEnd, r, pStart, pEnd, p, len, maxRowLen = 0;

630:   PetscFunctionBegin;
631:   PetscCall(PetscLayoutGetRange(rLayout, &rStart, &rEnd));
632:   for (r = rStart; r < rEnd; ++r) {
633:     PetscCall(PetscSectionGetDof(sectionAdj, r, &len));
634:     maxRowLen = PetscMax(maxRowLen, len);
635:   }
636:   PetscCall(PetscCalloc1(maxRowLen, &values));
637:   if (f >= 0 && bs == 1) {
638:     PetscCall(DMGetLocalSection(dm, &section));
639:     PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
640:     for (p = pStart; p < pEnd; ++p) {
641:       PetscInt rS, rE;

643:       PetscCall(DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE));
644:       for (r = rS; r < rE; ++r) {
645:         PetscInt numCols, cStart;

647:         PetscCall(PetscSectionGetDof(sectionAdj, r, &numCols));
648:         PetscCall(PetscSectionGetOffset(sectionAdj, r, &cStart));
649:         PetscCall(MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES));
650:       }
651:     }
652:   } else {
653:     for (r = rStart; r < rEnd; ++r) {
654:       PetscInt numCols, cStart;

656:       PetscCall(PetscSectionGetDof(sectionAdj, r, &numCols));
657:       PetscCall(PetscSectionGetOffset(sectionAdj, r, &cStart));
658:       PetscCall(MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES));
659:     }
660:   }
661:   PetscCall(PetscFree(values));
662:   PetscFunctionReturn(PETSC_SUCCESS);
663: }

665: /*@
666:   DMPlexPreallocateOperator - Calculate the matrix nonzero pattern based upon the information in the `DM`,
667:   the `PetscDS` it contains, and the default `PetscSection`.

669:   Collective

671:   Input Parameters:
672: + dm         - The `DMPLEX`
673: . bs         - The matrix blocksize
674: . dnz        - An array to hold the number of nonzeros in the diagonal block
675: . onz        - An array to hold the number of nonzeros in the off-diagonal block
676: . dnzu       - An array to hold the number of nonzeros in the upper triangle of the diagonal block
677: . onzu       - An array to hold the number of nonzeros in the upper triangle of the off-diagonal block
678: - fillMatrix - If `PETSC_TRUE`, fill the matrix with zeros

680:   Output Parameter:
681: . A - The preallocated matrix

683:   Level: advanced

685: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreateMatrix()`
686: @*/
687: PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
688: {
689:   MPI_Comm     comm;
690:   MatType      mtype;
691:   PetscSF      sf, sfDof;
692:   PetscSection section;
693:   PetscInt    *remoteOffsets;
694:   PetscSection sectionAdj[4] = {NULL, NULL, NULL, NULL};
695:   PetscInt    *cols[4]       = {NULL, NULL, NULL, NULL};
696:   PetscBool    useCone, useClosure;
697:   PetscInt     Nf, f, idx, locRows;
698:   PetscLayout  rLayout;
699:   PetscBool    isSymBlock, isSymSeqBlock, isSymMPIBlock, debug = PETSC_FALSE;
700:   PetscMPIInt  size;

702:   PetscFunctionBegin;
705:   if (dnz) PetscAssertPointer(dnz, 3);
706:   if (onz) PetscAssertPointer(onz, 4);
707:   if (dnzu) PetscAssertPointer(dnzu, 5);
708:   if (onzu) PetscAssertPointer(onzu, 6);
709:   PetscCall(DMGetPointSF(dm, &sf));
710:   PetscCall(DMGetLocalSection(dm, &section));
711:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_view_preallocation", &debug, NULL));
712:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
713:   PetscCallMPI(MPI_Comm_size(comm, &size));
714:   PetscCall(PetscLogEventBegin(DMPLEX_Preallocate, dm, 0, 0, 0));
715:   /* Create dof SF based on point SF */
716:   if (debug) {
717:     PetscSection section, sectionGlobal;
718:     PetscSF      sf;

720:     PetscCall(DMGetPointSF(dm, &sf));
721:     PetscCall(DMGetLocalSection(dm, &section));
722:     PetscCall(DMGetGlobalSection(dm, &sectionGlobal));
723:     PetscCall(PetscPrintf(comm, "Input Section for Preallocation:\n"));
724:     PetscCall(PetscSectionView(section, NULL));
725:     PetscCall(PetscPrintf(comm, "Input Global Section for Preallocation:\n"));
726:     PetscCall(PetscSectionView(sectionGlobal, NULL));
727:     if (size > 1) {
728:       PetscCall(PetscPrintf(comm, "Input SF for Preallocation:\n"));
729:       PetscCall(PetscSFView(sf, NULL));
730:     }
731:   }
732:   PetscCall(PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets));
733:   PetscCall(PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof));
734:   PetscCall(PetscFree(remoteOffsets));
735:   if (debug && size > 1) {
736:     PetscCall(PetscPrintf(comm, "Dof SF for Preallocation:\n"));
737:     PetscCall(PetscSFView(sfDof, NULL));
738:   }
739:   /* Create allocation vectors from adjacency graph */
740:   PetscCall(MatGetLocalSize(A, &locRows, NULL));
741:   PetscCall(PetscLayoutCreate(comm, &rLayout));
742:   PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
743:   PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
744:   PetscCall(PetscLayoutSetUp(rLayout));
745:   /* There are 4 types of adjacency */
746:   PetscCall(PetscSectionGetNumFields(section, &Nf));
747:   if (Nf < 1 || bs > 1) {
748:     PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure));
749:     idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
750:     PetscCall(DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, &sectionAdj[idx], &cols[idx]));
751:     PetscCall(DMPlexUpdateAllocation_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu));
752:   } else {
753:     for (f = 0; f < Nf; ++f) {
754:       PetscCall(DMGetAdjacency(dm, f, &useCone, &useClosure));
755:       idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
756:       if (!sectionAdj[idx]) PetscCall(DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, &sectionAdj[idx], &cols[idx]));
757:       PetscCall(DMPlexUpdateAllocation_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu));
758:     }
759:   }
760:   PetscCall(PetscSFDestroy(&sfDof));
761:   /* Set matrix pattern */
762:   PetscCall(MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu));
763:   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
764:   /* Check for symmetric storage */
765:   PetscCall(MatGetType(A, &mtype));
766:   PetscCall(PetscStrcmp(mtype, MATSBAIJ, &isSymBlock));
767:   PetscCall(PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock));
768:   PetscCall(PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock));
769:   if (isSymBlock || isSymSeqBlock || isSymMPIBlock) PetscCall(MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE));
770:   /* Fill matrix with zeros */
771:   if (fillMatrix) {
772:     if (Nf < 1 || bs > 1) {
773:       PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure));
774:       idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
775:       PetscCall(DMPlexFillMatrix_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], A));
776:     } else {
777:       for (f = 0; f < Nf; ++f) {
778:         PetscCall(DMGetAdjacency(dm, f, &useCone, &useClosure));
779:         idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
780:         PetscCall(DMPlexFillMatrix_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], A));
781:       }
782:     }
783:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
784:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
785:   }
786:   PetscCall(PetscLayoutDestroy(&rLayout));
787:   for (idx = 0; idx < 4; ++idx) {
788:     PetscCall(PetscSectionDestroy(&sectionAdj[idx]));
789:     PetscCall(PetscFree(cols[idx]));
790:   }
791:   PetscCall(PetscLogEventEnd(DMPLEX_Preallocate, dm, 0, 0, 0));
792:   PetscFunctionReturn(PETSC_SUCCESS);
793: }

795: #if 0
796: PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
797: {
798:   PetscInt       *tmpClosure,*tmpAdj,*visits;
799:   PetscInt        c,cStart,cEnd,pStart,pEnd;

801:   PetscFunctionBegin;
802:   PetscCall(DMGetDimension(dm, &dim));
803:   PetscCall(DMPlexGetDepth(dm, &depth));
804:   PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize));

806:   maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1));

808:   PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
809:   npoints = pEnd - pStart;

811:   PetscCall(PetscMalloc3(maxClosureSize,&tmpClosure,npoints,&lvisits,npoints,&visits));
812:   PetscCall(PetscArrayzero(lvisits,pEnd-pStart));
813:   PetscCall(PetscArrayzero(visits,pEnd-pStart));
814:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
815:   for (c=cStart; c<cEnd; c++) {
816:     PetscInt *support = tmpClosure;
817:     PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support));
818:     for (p=0; p<supportSize; p++) lvisits[support[p]]++;
819:   }
820:   PetscCall(PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM));
821:   PetscCall(PetscSFReduceEnd  (sf,MPIU_INT,lvisits,visits,MPI_SUM));
822:   PetscCall(PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits,MPI_REPLACE));
823:   PetscCall(PetscSFBcastEnd  (sf,MPIU_INT,visits,lvisits));

825:   PetscCall(PetscSFGetRootRanks());

827:   PetscCall(PetscMalloc2(maxClosureSize*maxClosureSize,&cellmat,npoints,&owner));
828:   for (c=cStart; c<cEnd; c++) {
829:     PetscCall(PetscArrayzero(cellmat,maxClosureSize*maxClosureSize));
830:     /*
831:      Depth-first walk of transitive closure.
832:      At each leaf frame f of transitive closure that we see, add 1/visits[f] to each pair (p,q) not marked as done in cellmat.
833:      This contribution is added to dnz if owning ranks of p and q match, to onz otherwise.
834:      */
835:   }

837:   PetscCall(PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM));
838:   PetscCall(PetscSFReduceEnd  (sf,MPIU_INT,lonz,onz,MPI_SUM));
839:   PetscFunctionReturn(PETSC_SUCCESS);
840: }
841: #endif