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