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