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