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: for (PetscInt s = 0; s < indegree[p]; ++s, ++r) rootAdj[l + s] = remoteadj[r];
486: }
487: PetscCheck(r == radjsize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %" PetscInt_FMT " != %" PetscInt_FMT, r, radjsize);
488: PetscCheck(l == adjSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %" PetscInt_FMT " != %" PetscInt_FMT, l, adjSize);
489: PetscCall(PetscFree(remoteadj));
490: }
491: PetscCall(PetscSFDestroy(&sfAdj));
492: PetscCall(PetscFree(adj));
493: /* Debugging */
494: if (debug) {
495: PetscCall(PetscPrintf(comm, "Root adjacency indices after gather\n"));
496: PetscCall(PetscSectionArrayView(rootSectionAdj, rootAdj, PETSC_INT, NULL));
497: }
498: /* Add in local adjacency indices for owned dofs on interface (roots) */
499: for (p = pStart; p < pEnd; ++p) {
500: PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof, anOff;
502: PetscCall(PetscSectionGetDof(section, p, &dof));
503: PetscCall(PetscSectionGetOffset(section, p, &off));
504: if (!dof) continue;
505: PetscCall(PetscSectionGetDof(rootSectionAdj, off, &adof));
506: if (adof <= 0) continue;
507: PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj));
508: PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof));
509: PetscCall(PetscSectionGetOffset(anchorSectionAdj, p, &anOff));
510: for (d = off; d < off + dof; ++d) {
511: PetscInt adof, aoff, i;
513: PetscCall(PetscSectionGetDof(rootSectionAdj, d, &adof));
514: PetscCall(PetscSectionGetOffset(rootSectionAdj, d, &aoff));
515: i = adof - 1;
516: for (q = 0; q < anDof; q++) {
517: rootAdj[aoff + i] = anchorAdj[anOff + q];
518: --i;
519: }
520: for (q = 0; q < numAdj; ++q) {
521: const PetscInt padj = tmpAdj[q];
522: PetscInt ndof, ncdof, ngoff, nd;
524: if ((padj < pStart) || (padj >= pEnd)) continue;
525: PetscCall(PetscSectionGetDof(section, padj, &ndof));
526: PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof));
527: PetscCall(PetscSectionGetOffset(sectionGlobal, padj, &ngoff));
528: for (nd = 0; nd < ndof - ncdof; ++nd) {
529: rootAdj[aoff + i] = ngoff < 0 ? -(ngoff + 1) + nd : ngoff + nd;
530: --i;
531: }
532: }
533: }
534: }
535: /* Debugging */
536: if (debug) {
537: PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots before compression:\n"));
538: PetscCall(PetscSectionArrayView(rootSectionAdj, rootAdj, PETSC_INT, NULL));
539: }
540: /* Compress indices */
541: PetscCall(PetscSectionSetUp(rootSectionAdj));
542: for (p = pStart; p < pEnd; ++p) {
543: PetscInt dof, cdof, off, d;
544: PetscInt adof, aoff;
546: PetscCall(PetscSectionGetDof(section, p, &dof));
547: PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
548: PetscCall(PetscSectionGetOffset(section, p, &off));
549: if (!dof) continue;
550: PetscCall(PetscSectionGetDof(rootSectionAdj, off, &adof));
551: if (adof <= 0) continue;
552: for (d = off; d < off + dof - cdof; ++d) {
553: PetscCall(PetscSectionGetDof(rootSectionAdj, d, &adof));
554: PetscCall(PetscSectionGetOffset(rootSectionAdj, d, &aoff));
555: PetscCall(PetscSortRemoveDupsInt(&adof, &rootAdj[aoff]));
556: PetscCall(PetscSectionSetDof(rootSectionAdj, d, adof));
557: }
558: }
559: /* Debugging */
560: if (debug) {
561: PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots after compression:\n"));
562: PetscCall(PetscSectionArrayView(rootSectionAdj, rootAdj, PETSC_INT, NULL));
563: }
564: /* Build adjacency section: Maps global indices to sets of adjacent global indices */
565: PetscCall(PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd));
566: PetscCall(PetscSectionCreate(comm, §ionAdj));
567: PetscCall(PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd));
569: PetscInt *exclude_leaves, num_exclude_leaves = 0, max_adjacency_size = 0;
570: PetscCall(DMPlexGetMaxAdjacencySize_Internal(dm, useAnchors, &max_adjacency_size));
571: PetscCall(PetscMalloc1(max_adjacency_size, &exclude_leaves));
573: for (p = pStart; p < pEnd; ++p) {
574: PetscInt numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof;
575: PetscBool found = PETSC_TRUE;
577: PetscCall(PetscSectionGetDof(section, p, &dof));
578: PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
579: PetscCall(PetscSectionGetOffset(section, p, &off));
580: PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff));
581: for (d = 0; d < dof - cdof; ++d) {
582: PetscInt ldof, rdof;
584: PetscCall(PetscSectionGetDof(leafSectionAdj, off + d, &ldof));
585: PetscCall(PetscSectionGetDof(rootSectionAdj, off + d, &rdof));
586: if (ldof > 0) {
587: /* We do not own this point */
588: } else if (rdof > 0) {
589: PetscCall(PetscSectionSetDof(sectionAdj, goff + d, rdof));
590: } else {
591: found = PETSC_FALSE;
592: }
593: }
594: if (found) continue;
595: PetscCall(PetscSectionGetDof(section, p, &dof));
596: PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff));
597: PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj));
598: PetscCall(AdjancencyContainsLeafRootPair(myRankPairSection, myRankPairLeaves, myRankPairRoots, numAdj, tmpAdj, &num_exclude_leaves, exclude_leaves));
599: for (q = 0; q < numAdj; ++q) {
600: const PetscInt padj = tmpAdj[q];
601: PetscInt ndof, ncdof, noff, count;
603: if ((padj < pStart) || (padj >= pEnd)) continue;
604: PetscCall(PetscSectionGetDof(section, padj, &ndof));
605: PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof));
606: PetscCall(PetscSectionGetOffset(section, padj, &noff));
607: // If leaf-root pair are both on this rank, only count root
608: PetscCall(PetscFindInt(padj, num_exclude_leaves, exclude_leaves, &count));
609: if (count >= 0) continue;
610: for (d = goff; d < goff + dof - cdof; ++d) PetscCall(PetscSectionAddDof(sectionAdj, d, ndof - ncdof));
611: }
612: PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof));
613: if (anDof) {
614: for (d = goff; d < goff + dof - cdof; ++d) PetscCall(PetscSectionAddDof(sectionAdj, d, anDof));
615: }
616: }
617: PetscCall(PetscSectionSetUp(sectionAdj));
618: if (debug) {
619: PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation:\n"));
620: PetscCall(PetscSectionView(sectionAdj, NULL));
621: }
622: /* Get adjacent indices */
623: PetscCall(PetscSectionGetStorageSize(sectionAdj, &numCols));
624: PetscCall(PetscMalloc1(numCols, &cols));
625: for (p = pStart; p < pEnd; ++p) {
626: PetscInt numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof, anOff;
627: PetscBool found = PETSC_TRUE;
629: PetscCall(PetscSectionGetDof(section, p, &dof));
630: PetscCall(PetscSectionGetConstraintDof(section, p, &cdof));
631: PetscCall(PetscSectionGetOffset(section, p, &off));
632: PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff));
633: for (d = 0; d < dof - cdof; ++d) {
634: PetscInt ldof, rdof;
636: PetscCall(PetscSectionGetDof(leafSectionAdj, off + d, &ldof));
637: PetscCall(PetscSectionGetDof(rootSectionAdj, off + d, &rdof));
638: if (ldof > 0) {
639: /* We do not own this point */
640: } else if (rdof > 0) {
641: PetscInt aoff, roff;
643: PetscCall(PetscSectionGetOffset(sectionAdj, goff + d, &aoff));
644: PetscCall(PetscSectionGetOffset(rootSectionAdj, off + d, &roff));
645: PetscCall(PetscArraycpy(&cols[aoff], &rootAdj[roff], rdof));
646: } else {
647: found = PETSC_FALSE;
648: }
649: }
650: if (found) continue;
651: PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj));
652: PetscCall(AdjancencyContainsLeafRootPair(myRankPairSection, myRankPairLeaves, myRankPairRoots, numAdj, tmpAdj, &num_exclude_leaves, exclude_leaves));
653: PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof));
654: PetscCall(PetscSectionGetOffset(anchorSectionAdj, p, &anOff));
655: for (d = goff; d < goff + dof - cdof; ++d) {
656: PetscInt adof, aoff, i = 0;
658: PetscCall(PetscSectionGetDof(sectionAdj, d, &adof));
659: PetscCall(PetscSectionGetOffset(sectionAdj, d, &aoff));
660: for (q = 0; q < numAdj; ++q) {
661: const PetscInt padj = tmpAdj[q], *ncind;
662: PetscInt ndof, ncdof, ngoff, nd, count;
664: /* Adjacent points may not be in the section chart */
665: if ((padj < pStart) || (padj >= pEnd)) continue;
666: PetscCall(PetscSectionGetDof(section, padj, &ndof));
667: PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof));
668: PetscCall(PetscSectionGetConstraintIndices(section, padj, &ncind));
669: PetscCall(PetscSectionGetOffset(sectionGlobal, padj, &ngoff));
670: // If leaf-root pair are both on this rank, only count root
671: PetscCall(PetscFindInt(padj, num_exclude_leaves, exclude_leaves, &count));
672: if (count >= 0) continue;
673: for (nd = 0; nd < ndof - ncdof; ++nd, ++i) cols[aoff + i] = ngoff < 0 ? -(ngoff + 1) + nd : ngoff + nd;
674: }
675: for (q = 0; q < anDof; q++, i++) cols[aoff + i] = anchorAdj[anOff + q];
676: 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);
677: }
678: }
679: PetscCall(PetscSectionDestroy(&anchorSectionAdj));
680: PetscCall(PetscSectionDestroy(&leafSectionAdj));
681: PetscCall(PetscSectionDestroy(&rootSectionAdj));
682: PetscCall(PetscSectionDestroy(&myRankPairSection));
683: PetscCall(PetscFree(exclude_leaves));
684: PetscCall(PetscFree(myRankPairLeaves));
685: PetscCall(PetscFree(myRankPairRoots));
686: PetscCall(PetscFree(anchorAdj));
687: PetscCall(PetscFree(rootAdj));
688: PetscCall(PetscFree(tmpAdj));
689: /* Debugging */
690: if (debug) {
691: PetscCall(PetscPrintf(comm, "Column indices\n"));
692: PetscCall(PetscSectionArrayView(sectionAdj, cols, PETSC_INT, NULL));
693: }
695: *sA = sectionAdj;
696: *colIdx = cols;
697: PetscFunctionReturn(PETSC_SUCCESS);
698: }
700: 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[])
701: {
702: PetscSection section;
703: PetscInt rStart, rEnd, r, pStart, pEnd, p;
705: PetscFunctionBegin;
706: /* This loop needs to change to a loop over points, then field dofs, which means we need to look both sections */
707: PetscCall(PetscLayoutGetRange(rLayout, &rStart, &rEnd));
708: 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);
709: if (f >= 0 && bs == 1) {
710: PetscCall(DMGetLocalSection(dm, §ion));
711: PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
712: for (p = pStart; p < pEnd; ++p) {
713: PetscInt rS, rE;
715: PetscCall(DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE));
716: for (r = rS; r < rE; ++r) {
717: PetscInt numCols, cStart, c;
719: PetscCall(PetscSectionGetDof(sectionAdj, r, &numCols));
720: PetscCall(PetscSectionGetOffset(sectionAdj, r, &cStart));
721: for (c = cStart; c < cStart + numCols; ++c) {
722: if ((cols[c] >= rStart) && (cols[c] < rEnd)) {
723: ++dnz[r - rStart];
724: if (cols[c] >= r) ++dnzu[r - rStart];
725: } else {
726: ++onz[r - rStart];
727: if (cols[c] >= r) ++onzu[r - rStart];
728: }
729: }
730: }
731: }
732: } else {
733: /* Only loop over blocks of rows */
734: for (r = rStart / bs; r < rEnd / bs; ++r) {
735: const PetscInt row = r * bs;
736: PetscInt numCols, cStart, c;
738: PetscCall(PetscSectionGetDof(sectionAdj, row, &numCols));
739: PetscCall(PetscSectionGetOffset(sectionAdj, row, &cStart));
740: for (c = cStart; c < cStart + numCols; ++c) {
741: if ((cols[c] >= rStart) && (cols[c] < rEnd)) {
742: ++dnz[r - rStart / bs];
743: if (cols[c] >= row) ++dnzu[r - rStart / bs];
744: } else {
745: ++onz[r - rStart / bs];
746: if (cols[c] >= row) ++onzu[r - rStart / bs];
747: }
748: }
749: }
750: for (r = 0; r < (rEnd - rStart) / bs; ++r) {
751: dnz[r] /= bs;
752: onz[r] /= bs;
753: dnzu[r] /= bs;
754: onzu[r] /= bs;
755: }
756: }
757: PetscFunctionReturn(PETSC_SUCCESS);
758: }
760: static PetscErrorCode DMPlexFillMatrix_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], Mat A)
761: {
762: PetscSection section;
763: PetscScalar *values;
764: PetscInt rStart, rEnd, r, pStart, pEnd, p, len, maxRowLen = 0;
766: PetscFunctionBegin;
767: PetscCall(PetscLayoutGetRange(rLayout, &rStart, &rEnd));
768: for (r = rStart; r < rEnd; ++r) {
769: PetscCall(PetscSectionGetDof(sectionAdj, r, &len));
770: maxRowLen = PetscMax(maxRowLen, len);
771: }
772: PetscCall(PetscCalloc1(maxRowLen, &values));
773: if (f >= 0 && bs == 1) {
774: PetscCall(DMGetLocalSection(dm, §ion));
775: PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
776: for (p = pStart; p < pEnd; ++p) {
777: PetscInt rS, rE;
779: PetscCall(DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE));
780: for (r = rS; r < rE; ++r) {
781: PetscInt numCols, cStart;
783: PetscCall(PetscSectionGetDof(sectionAdj, r, &numCols));
784: PetscCall(PetscSectionGetOffset(sectionAdj, r, &cStart));
785: PetscCall(MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES));
786: }
787: }
788: } else {
789: for (r = rStart; r < rEnd; ++r) {
790: PetscInt numCols, cStart;
792: PetscCall(PetscSectionGetDof(sectionAdj, r, &numCols));
793: PetscCall(PetscSectionGetOffset(sectionAdj, r, &cStart));
794: PetscCall(MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES));
795: }
796: }
797: PetscCall(PetscFree(values));
798: PetscFunctionReturn(PETSC_SUCCESS);
799: }
801: /*@
802: DMPlexPreallocateOperator - Calculate the matrix nonzero pattern based upon the information in the `DM`,
803: the `PetscDS` it contains, and the default `PetscSection`.
805: Collective
807: Input Parameters:
808: + dm - The `DMPLEX`
809: . bs - The matrix blocksize
810: . dnz - An array to hold the number of nonzeros in the diagonal block
811: . onz - An array to hold the number of nonzeros in the off-diagonal block
812: . dnzu - An array to hold the number of nonzeros in the upper triangle of the diagonal block
813: . onzu - An array to hold the number of nonzeros in the upper triangle of the off-diagonal block
814: - fillMatrix - If `PETSC_TRUE`, fill the matrix with zeros
816: Output Parameter:
817: . A - The preallocated matrix
819: Level: advanced
821: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreateMatrix()`
822: @*/
823: PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
824: {
825: MPI_Comm comm;
826: MatType mtype;
827: PetscSF sf, sfDof;
828: PetscSection section;
829: PetscInt *remoteOffsets;
830: PetscSection sectionAdj[4] = {NULL, NULL, NULL, NULL};
831: PetscInt *cols[4] = {NULL, NULL, NULL, NULL};
832: PetscBool useCone, useClosure;
833: PetscInt Nf, f, idx, locRows;
834: PetscLayout rLayout;
835: PetscBool isSymBlock, isSymSeqBlock, isSymMPIBlock, debug = PETSC_FALSE;
837: PetscFunctionBegin;
840: if (dnz) PetscAssertPointer(dnz, 3);
841: if (onz) PetscAssertPointer(onz, 4);
842: if (dnzu) PetscAssertPointer(dnzu, 5);
843: if (onzu) PetscAssertPointer(onzu, 6);
844: PetscCall(DMGetIsoperiodicPointSF_Internal(dm, &sf));
845: PetscCall(DMGetLocalSection(dm, §ion));
846: PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_view_preallocation", &debug, NULL));
847: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
848: PetscCall(PetscLogEventBegin(DMPLEX_Preallocate, dm, 0, 0, 0));
849: /* Create dof SF based on point SF */
850: if (debug) {
851: PetscSection section, sectionGlobal;
852: PetscSF sf;
854: PetscCall(DMGetIsoperiodicPointSF_Internal(dm, &sf));
855: PetscCall(DMGetLocalSection(dm, §ion));
856: PetscCall(DMGetGlobalSection(dm, §ionGlobal));
857: PetscCall(PetscPrintf(comm, "Input Section for Preallocation:\n"));
858: PetscCall(PetscSectionView(section, NULL));
859: PetscCall(PetscPrintf(comm, "Input Global Section for Preallocation:\n"));
860: PetscCall(PetscSectionView(sectionGlobal, NULL));
861: PetscCall(PetscPrintf(comm, "Input SF for Preallocation:\n"));
862: PetscCall(PetscSFView(sf, NULL));
863: }
864: PetscCall(PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets));
865: PetscCall(PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof));
866: PetscCall(PetscFree(remoteOffsets));
867: if (debug) {
868: PetscCall(PetscPrintf(comm, "Dof SF for Preallocation:\n"));
869: PetscCall(PetscSFView(sfDof, NULL));
870: }
871: /* Create allocation vectors from adjacency graph */
872: PetscCall(MatGetLocalSize(A, &locRows, NULL));
873: PetscCall(PetscLayoutCreate(comm, &rLayout));
874: PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
875: PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
876: PetscCall(PetscLayoutSetUp(rLayout));
877: /* There are 4 types of adjacency */
878: PetscCall(PetscSectionGetNumFields(section, &Nf));
879: if (Nf < 1 || bs > 1) {
880: PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure));
881: idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
882: PetscCall(DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, §ionAdj[idx], &cols[idx]));
883: PetscCall(DMPlexUpdateAllocation_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu));
884: } else {
885: for (f = 0; f < Nf; ++f) {
886: PetscCall(DMGetAdjacency(dm, f, &useCone, &useClosure));
887: idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
888: if (!sectionAdj[idx]) PetscCall(DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, §ionAdj[idx], &cols[idx]));
889: PetscCall(DMPlexUpdateAllocation_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu));
890: }
891: }
892: PetscCall(PetscSFDestroy(&sfDof));
893: /* Set matrix pattern */
894: PetscCall(MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu));
895: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
896: /* Check for symmetric storage */
897: PetscCall(MatGetType(A, &mtype));
898: PetscCall(PetscStrcmp(mtype, MATSBAIJ, &isSymBlock));
899: PetscCall(PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock));
900: PetscCall(PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock));
901: if (isSymBlock || isSymSeqBlock || isSymMPIBlock) PetscCall(MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE));
902: /* Fill matrix with zeros */
903: if (fillMatrix) {
904: if (Nf < 1 || bs > 1) {
905: PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure));
906: idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
907: PetscCall(DMPlexFillMatrix_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], A));
908: } else {
909: for (f = 0; f < Nf; ++f) {
910: PetscCall(DMGetAdjacency(dm, f, &useCone, &useClosure));
911: idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
912: PetscCall(DMPlexFillMatrix_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], A));
913: }
914: }
915: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
916: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
917: }
918: PetscCall(PetscLayoutDestroy(&rLayout));
919: for (idx = 0; idx < 4; ++idx) {
920: PetscCall(PetscSectionDestroy(§ionAdj[idx]));
921: PetscCall(PetscFree(cols[idx]));
922: }
923: PetscCall(PetscLogEventEnd(DMPLEX_Preallocate, dm, 0, 0, 0));
924: PetscFunctionReturn(PETSC_SUCCESS);
925: }
927: #if 0
928: PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
929: {
930: PetscInt *tmpClosure,*tmpAdj,*visits;
931: PetscInt c,cStart,cEnd,pStart,pEnd;
933: PetscFunctionBegin;
934: PetscCall(DMGetDimension(dm, &dim));
935: PetscCall(DMPlexGetDepth(dm, &depth));
936: PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize));
938: maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1));
940: PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
941: npoints = pEnd - pStart;
943: PetscCall(PetscMalloc3(maxClosureSize,&tmpClosure,npoints,&lvisits,npoints,&visits));
944: PetscCall(PetscArrayzero(lvisits,pEnd-pStart));
945: PetscCall(PetscArrayzero(visits,pEnd-pStart));
946: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
947: for (c=cStart; c<cEnd; c++) {
948: PetscInt *support = tmpClosure;
949: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support));
950: for (p=0; p<supportSize; p++) lvisits[support[p]]++;
951: }
952: PetscCall(PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM));
953: PetscCall(PetscSFReduceEnd (sf,MPIU_INT,lvisits,visits,MPI_SUM));
954: PetscCall(PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits,MPI_REPLACE));
955: PetscCall(PetscSFBcastEnd (sf,MPIU_INT,visits,lvisits));
957: PetscCall(PetscSFGetRootRanks());
959: PetscCall(PetscMalloc2(maxClosureSize*maxClosureSize,&cellmat,npoints,&owner));
960: for (c=cStart; c<cEnd; c++) {
961: PetscCall(PetscArrayzero(cellmat,maxClosureSize*maxClosureSize));
962: /*
963: Depth-first walk of transitive closure.
964: 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.
965: This contribution is added to dnz if owning ranks of p and q match, to onz otherwise.
966: */
967: }
969: PetscCall(PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM));
970: PetscCall(PetscSFReduceEnd (sf,MPIU_INT,lonz,onz,MPI_SUM));
971: PetscFunctionReturn(PETSC_SUCCESS);
972: }
973: #endif