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: DMGetLocalSection(dm, §ion);
14: DMGetGlobalSection(dm, §ionGlobal);
15: PetscSectionCreate(PetscObjectComm((PetscObject)section), &adjSec);
16: PetscSectionGetChart(section, &pStart, &pEnd);
17: PetscSectionSetChart(adjSec, pStart, pEnd);
19: DMPlexGetAnchors(dm, &aSec, &aIS);
20: if (aSec) {
21: const PetscInt *anchors;
22: PetscInt p, q, a, aSize, *offsets, aStart, aEnd, *inverse, iSize, *adj, adjSize;
23: PetscInt *tmpAdjP = NULL, *tmpAdjQ = NULL;
24: PetscSection inverseSec;
26: /* invert the constraint-to-anchor map */
27: PetscSectionCreate(PetscObjectComm((PetscObject)aSec), &inverseSec);
28: PetscSectionSetChart(inverseSec, pStart, pEnd);
29: ISGetLocalSize(aIS, &aSize);
30: ISGetIndices(aIS, &anchors);
32: for (p = 0; p < aSize; p++) {
33: PetscInt a = anchors[p];
35: PetscSectionAddDof(inverseSec, a, 1);
36: }
37: PetscSectionSetUp(inverseSec);
38: PetscSectionGetStorageSize(inverseSec, &iSize);
39: PetscMalloc1(iSize, &inverse);
40: PetscCalloc1(pEnd - pStart, &offsets);
41: PetscSectionGetChart(aSec, &aStart, &aEnd);
42: for (p = aStart; p < aEnd; p++) {
43: PetscInt dof, off;
45: PetscSectionGetDof(aSec, p, &dof);
46: PetscSectionGetOffset(aSec, p, &off);
48: for (q = 0; q < dof; q++) {
49: PetscInt iOff;
51: a = anchors[off + q];
52: PetscSectionGetOffset(inverseSec, a, &iOff);
53: inverse[iOff + offsets[a - pStart]++] = p;
54: }
55: }
56: ISRestoreIndices(aIS, &anchors);
57: PetscFree(offsets);
59: /* construct anchorAdj and adjSec
60: *
61: * loop over anchors:
62: * construct anchor adjacency
63: * loop over constrained:
64: * construct constrained adjacency
65: * if not in anchor adjacency, add to dofs
66: * setup adjSec, allocate anchorAdj
67: * loop over anchors:
68: * construct anchor adjacency
69: * loop over constrained:
70: * construct constrained adjacency
71: * if not in anchor adjacency
72: * if not already in list, put in list
73: * sort, unique, reduce dof count
74: * optional: compactify
75: */
76: for (p = pStart; p < pEnd; p++) {
77: PetscInt iDof, iOff, i, r, s, numAdjP = PETSC_DETERMINE;
79: PetscSectionGetDof(inverseSec, p, &iDof);
80: if (!iDof) continue;
81: PetscSectionGetOffset(inverseSec, p, &iOff);
82: DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, PETSC_TRUE, &numAdjP, &tmpAdjP);
83: for (i = 0; i < iDof; i++) {
84: PetscInt iNew = 0, qAdj, qAdjDof, qAdjCDof, numAdjQ = PETSC_DETERMINE;
86: q = inverse[iOff + i];
87: DMPlexGetAdjacency_Internal(dm, q, useCone, useClosure, PETSC_TRUE, &numAdjQ, &tmpAdjQ);
88: for (r = 0; r < numAdjQ; r++) {
89: qAdj = tmpAdjQ[r];
90: if ((qAdj < pStart) || (qAdj >= pEnd)) continue;
91: for (s = 0; s < numAdjP; s++) {
92: if (qAdj == tmpAdjP[s]) break;
93: }
94: if (s < numAdjP) continue;
95: PetscSectionGetDof(section, qAdj, &qAdjDof);
96: PetscSectionGetConstraintDof(section, qAdj, &qAdjCDof);
97: iNew += qAdjDof - qAdjCDof;
98: }
99: PetscSectionAddDof(adjSec, p, iNew);
100: }
101: }
103: PetscSectionSetUp(adjSec);
104: PetscSectionGetStorageSize(adjSec, &adjSize);
105: PetscMalloc1(adjSize, &adj);
107: for (p = pStart; p < pEnd; p++) {
108: PetscInt iDof, iOff, i, r, s, aOff, aOffOrig, aDof, numAdjP = PETSC_DETERMINE;
110: PetscSectionGetDof(inverseSec, p, &iDof);
111: if (!iDof) continue;
112: PetscSectionGetOffset(inverseSec, p, &iOff);
113: DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, PETSC_TRUE, &numAdjP, &tmpAdjP);
114: PetscSectionGetDof(adjSec, p, &aDof);
115: PetscSectionGetOffset(adjSec, p, &aOff);
116: aOffOrig = aOff;
117: for (i = 0; i < iDof; i++) {
118: PetscInt qAdj, qAdjDof, qAdjCDof, qAdjOff, nd, numAdjQ = PETSC_DETERMINE;
120: q = inverse[iOff + i];
121: DMPlexGetAdjacency_Internal(dm, q, useCone, useClosure, PETSC_TRUE, &numAdjQ, &tmpAdjQ);
122: for (r = 0; r < numAdjQ; r++) {
123: qAdj = tmpAdjQ[r];
124: if ((qAdj < pStart) || (qAdj >= pEnd)) continue;
125: for (s = 0; s < numAdjP; s++) {
126: if (qAdj == tmpAdjP[s]) break;
127: }
128: if (s < numAdjP) continue;
129: PetscSectionGetDof(section, qAdj, &qAdjDof);
130: PetscSectionGetConstraintDof(section, qAdj, &qAdjCDof);
131: PetscSectionGetOffset(sectionGlobal, qAdj, &qAdjOff);
132: for (nd = 0; nd < qAdjDof - qAdjCDof; ++nd) adj[aOff++] = (qAdjOff < 0 ? -(qAdjOff + 1) : qAdjOff) + nd;
133: }
134: }
135: PetscSortRemoveDupsInt(&aDof, &adj[aOffOrig]);
136: PetscSectionSetDof(adjSec, p, aDof);
137: }
138: *anchorAdj = adj;
140: /* clean up */
141: PetscSectionDestroy(&inverseSec);
142: PetscFree(inverse);
143: PetscFree(tmpAdjP);
144: PetscFree(tmpAdjQ);
145: } else {
146: *anchorAdj = NULL;
147: PetscSectionSetUp(adjSec);
148: }
149: *anchorSectionAdj = adjSec;
150: return 0;
151: }
153: static PetscErrorCode DMPlexCreateAdjacencySection_Static(DM dm, PetscInt bs, PetscSF sfDof, PetscBool useCone, PetscBool useClosure, PetscBool useAnchors, PetscSection *sA, PetscInt **colIdx)
154: {
155: MPI_Comm comm;
156: PetscMPIInt size;
157: PetscBool doCommLocal, doComm, debug = PETSC_FALSE;
158: PetscSF sf, sfAdj;
159: PetscSection section, sectionGlobal, leafSectionAdj, rootSectionAdj, sectionAdj, anchorSectionAdj;
160: PetscInt nroots, nleaves, l, p, r;
161: const PetscInt *leaves;
162: const PetscSFNode *remotes;
163: PetscInt dim, pStart, pEnd, numDof, globalOffStart, globalOffEnd, numCols;
164: PetscInt *tmpAdj = NULL, *adj, *rootAdj, *anchorAdj = NULL, *cols, *remoteOffsets;
165: PetscInt adjSize;
167: PetscObjectGetComm((PetscObject)dm, &comm);
168: PetscOptionsGetBool(NULL, NULL, "-dm_view_preallocation", &debug, NULL);
169: MPI_Comm_size(comm, &size);
170: DMGetDimension(dm, &dim);
171: DMGetPointSF(dm, &sf);
172: DMGetLocalSection(dm, §ion);
173: DMGetGlobalSection(dm, §ionGlobal);
174: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
175: doCommLocal = (size > 1) && (nroots >= 0) ? PETSC_TRUE : PETSC_FALSE;
176: MPIU_Allreduce(&doCommLocal, &doComm, 1, MPIU_BOOL, MPI_LAND, comm);
177: /* Create section for dof adjacency (dof ==> # adj dof) */
178: PetscSectionGetChart(section, &pStart, &pEnd);
179: PetscSectionGetStorageSize(section, &numDof);
180: PetscSectionCreate(comm, &leafSectionAdj);
181: PetscSectionSetChart(leafSectionAdj, 0, numDof);
182: PetscSectionCreate(comm, &rootSectionAdj);
183: PetscSectionSetChart(rootSectionAdj, 0, numDof);
184: /* Fill in the ghost dofs on the interface */
185: PetscSFGetGraph(sf, NULL, &nleaves, &leaves, &remotes);
186: /*
187: section - maps points to (# dofs, local dofs)
188: sectionGlobal - maps points to (# dofs, global dofs)
189: leafSectionAdj - maps unowned local dofs to # adj dofs
190: rootSectionAdj - maps owned local dofs to # adj dofs
191: adj - adj global dofs indexed by leafSectionAdj
192: rootAdj - adj global dofs indexed by rootSectionAdj
193: sf - describes shared points across procs
194: sfDof - describes shared dofs across procs
195: sfAdj - describes shared adjacent dofs across procs
196: ** The bootstrapping process involves six rounds with similar structure of visiting neighbors of each point.
197: (0). If there are point-to-point constraints, add the adjacencies of constrained points to anchors in anchorAdj
198: (This is done in DMPlexComputeAnchorAdjacencies())
199: 1. Visit unowned points on interface, count adjacencies placing in leafSectionAdj
200: Reduce those counts to rootSectionAdj (now redundantly counting some interface points)
201: 2. Visit owned points on interface, count adjacencies placing in rootSectionAdj
202: Create sfAdj connecting rootSectionAdj and leafSectionAdj
203: 3. Visit unowned points on interface, write adjacencies to adj
204: Gather adj to rootAdj (note that there is redundancy in rootAdj when multiple procs find the same adjacencies)
205: 4. Visit owned points on interface, write adjacencies to rootAdj
206: Remove redundancy in rootAdj
207: ** The last two traversals use transitive closure
208: 5. Visit all owned points in the subdomain, count dofs for each point (sectionAdj)
209: Allocate memory addressed by sectionAdj (cols)
210: 6. Visit all owned points in the subdomain, insert dof adjacencies into cols
211: ** Knowing all the column adjacencies, check ownership and sum into dnz and onz
212: */
213: DMPlexComputeAnchorAdjacencies(dm, useCone, useClosure, &anchorSectionAdj, &anchorAdj);
214: for (l = 0; l < nleaves; ++l) {
215: PetscInt dof, off, d, q, anDof;
216: PetscInt p = leaves[l], numAdj = PETSC_DETERMINE;
218: if ((p < pStart) || (p >= pEnd)) continue;
219: PetscSectionGetDof(section, p, &dof);
220: PetscSectionGetOffset(section, p, &off);
221: DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
222: for (q = 0; q < numAdj; ++q) {
223: const PetscInt padj = tmpAdj[q];
224: PetscInt ndof, ncdof;
226: if ((padj < pStart) || (padj >= pEnd)) continue;
227: PetscSectionGetDof(section, padj, &ndof);
228: PetscSectionGetConstraintDof(section, padj, &ncdof);
229: for (d = off; d < off + dof; ++d) PetscSectionAddDof(leafSectionAdj, d, ndof - ncdof);
230: }
231: PetscSectionGetDof(anchorSectionAdj, p, &anDof);
232: if (anDof) {
233: for (d = off; d < off + dof; ++d) PetscSectionAddDof(leafSectionAdj, d, anDof);
234: }
235: }
236: PetscSectionSetUp(leafSectionAdj);
237: if (debug) {
238: PetscPrintf(comm, "Adjacency Section for Preallocation on Leaves:\n");
239: PetscSectionView(leafSectionAdj, NULL);
240: }
241: /* Get maximum remote adjacency sizes for owned dofs on interface (roots) */
242: if (doComm) {
243: PetscSFReduceBegin(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);
244: PetscSFReduceEnd(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);
245: PetscSectionInvalidateMaxDof_Internal(rootSectionAdj);
246: }
247: if (debug) {
248: PetscPrintf(comm, "Adjancency Section for Preallocation on Roots:\n");
249: PetscSectionView(rootSectionAdj, NULL);
250: }
251: /* Add in local adjacency sizes for owned dofs on interface (roots) */
252: for (p = pStart; p < pEnd; ++p) {
253: PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof;
255: PetscSectionGetDof(section, p, &dof);
256: PetscSectionGetOffset(section, p, &off);
257: if (!dof) continue;
258: PetscSectionGetDof(rootSectionAdj, off, &adof);
259: if (adof <= 0) continue;
260: DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
261: for (q = 0; q < numAdj; ++q) {
262: const PetscInt padj = tmpAdj[q];
263: PetscInt ndof, ncdof;
265: if ((padj < pStart) || (padj >= pEnd)) continue;
266: PetscSectionGetDof(section, padj, &ndof);
267: PetscSectionGetConstraintDof(section, padj, &ncdof);
268: for (d = off; d < off + dof; ++d) PetscSectionAddDof(rootSectionAdj, d, ndof - ncdof);
269: }
270: PetscSectionGetDof(anchorSectionAdj, p, &anDof);
271: if (anDof) {
272: for (d = off; d < off + dof; ++d) PetscSectionAddDof(rootSectionAdj, d, anDof);
273: }
274: }
275: PetscSectionSetUp(rootSectionAdj);
276: if (debug) {
277: PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after local additions:\n");
278: PetscSectionView(rootSectionAdj, NULL);
279: }
280: /* Create adj SF based on dof SF */
281: PetscSFCreateRemoteOffsets(sfDof, rootSectionAdj, leafSectionAdj, &remoteOffsets);
282: PetscSFCreateSectionSF(sfDof, rootSectionAdj, remoteOffsets, leafSectionAdj, &sfAdj);
283: PetscFree(remoteOffsets);
284: if (debug && size > 1) {
285: PetscPrintf(comm, "Adjacency SF for Preallocation:\n");
286: PetscSFView(sfAdj, NULL);
287: }
288: /* Create leaf adjacency */
289: PetscSectionSetUp(leafSectionAdj);
290: PetscSectionGetStorageSize(leafSectionAdj, &adjSize);
291: PetscCalloc1(adjSize, &adj);
292: for (l = 0; l < nleaves; ++l) {
293: PetscInt dof, off, d, q, anDof, anOff;
294: PetscInt p = leaves[l], numAdj = PETSC_DETERMINE;
296: if ((p < pStart) || (p >= pEnd)) continue;
297: PetscSectionGetDof(section, p, &dof);
298: PetscSectionGetOffset(section, p, &off);
299: DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
300: PetscSectionGetDof(anchorSectionAdj, p, &anDof);
301: PetscSectionGetOffset(anchorSectionAdj, p, &anOff);
302: for (d = off; d < off + dof; ++d) {
303: PetscInt aoff, i = 0;
305: PetscSectionGetOffset(leafSectionAdj, d, &aoff);
306: for (q = 0; q < numAdj; ++q) {
307: const PetscInt padj = tmpAdj[q];
308: PetscInt ndof, ncdof, ngoff, nd;
310: if ((padj < pStart) || (padj >= pEnd)) continue;
311: PetscSectionGetDof(section, padj, &ndof);
312: PetscSectionGetConstraintDof(section, padj, &ncdof);
313: PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
314: for (nd = 0; nd < ndof - ncdof; ++nd) {
315: adj[aoff + i] = (ngoff < 0 ? -(ngoff + 1) : ngoff) + nd;
316: ++i;
317: }
318: }
319: for (q = 0; q < anDof; q++) {
320: adj[aoff + i] = anchorAdj[anOff + q];
321: ++i;
322: }
323: }
324: }
325: /* Debugging */
326: if (debug) {
327: IS tmp;
328: PetscPrintf(comm, "Leaf adjacency indices\n");
329: ISCreateGeneral(comm, adjSize, adj, PETSC_USE_POINTER, &tmp);
330: ISView(tmp, NULL);
331: ISDestroy(&tmp);
332: }
333: /* Gather adjacent indices to root */
334: PetscSectionGetStorageSize(rootSectionAdj, &adjSize);
335: PetscMalloc1(adjSize, &rootAdj);
336: for (r = 0; r < adjSize; ++r) rootAdj[r] = -1;
337: if (doComm) {
338: const PetscInt *indegree;
339: PetscInt *remoteadj, radjsize = 0;
341: PetscSFComputeDegreeBegin(sfAdj, &indegree);
342: PetscSFComputeDegreeEnd(sfAdj, &indegree);
343: for (p = 0; p < adjSize; ++p) radjsize += indegree[p];
344: PetscMalloc1(radjsize, &remoteadj);
345: PetscSFGatherBegin(sfAdj, MPIU_INT, adj, remoteadj);
346: PetscSFGatherEnd(sfAdj, MPIU_INT, adj, remoteadj);
347: for (p = 0, l = 0, r = 0; p < adjSize; ++p, l = PetscMax(p, l + indegree[p - 1])) {
348: PetscInt s;
349: for (s = 0; s < indegree[p]; ++s, ++r) rootAdj[l + s] = remoteadj[r];
350: }
353: PetscFree(remoteadj);
354: }
355: PetscSFDestroy(&sfAdj);
356: PetscFree(adj);
357: /* Debugging */
358: if (debug) {
359: IS tmp;
360: PetscPrintf(comm, "Root adjacency indices after gather\n");
361: ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
362: ISView(tmp, NULL);
363: ISDestroy(&tmp);
364: }
365: /* Add in local adjacency indices for owned dofs on interface (roots) */
366: for (p = pStart; p < pEnd; ++p) {
367: PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof, anOff;
369: PetscSectionGetDof(section, p, &dof);
370: PetscSectionGetOffset(section, p, &off);
371: if (!dof) continue;
372: PetscSectionGetDof(rootSectionAdj, off, &adof);
373: if (adof <= 0) continue;
374: DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
375: PetscSectionGetDof(anchorSectionAdj, p, &anDof);
376: PetscSectionGetOffset(anchorSectionAdj, p, &anOff);
377: for (d = off; d < off + dof; ++d) {
378: PetscInt adof, aoff, i;
380: PetscSectionGetDof(rootSectionAdj, d, &adof);
381: PetscSectionGetOffset(rootSectionAdj, d, &aoff);
382: i = adof - 1;
383: for (q = 0; q < anDof; q++) {
384: rootAdj[aoff + i] = anchorAdj[anOff + q];
385: --i;
386: }
387: for (q = 0; q < numAdj; ++q) {
388: const PetscInt padj = tmpAdj[q];
389: PetscInt ndof, ncdof, ngoff, nd;
391: if ((padj < pStart) || (padj >= pEnd)) continue;
392: PetscSectionGetDof(section, padj, &ndof);
393: PetscSectionGetConstraintDof(section, padj, &ncdof);
394: PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
395: for (nd = 0; nd < ndof - ncdof; ++nd) {
396: rootAdj[aoff + i] = ngoff < 0 ? -(ngoff + 1) + nd : ngoff + nd;
397: --i;
398: }
399: }
400: }
401: }
402: /* Debugging */
403: if (debug) {
404: IS tmp;
405: PetscPrintf(comm, "Root adjacency indices\n");
406: ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
407: ISView(tmp, NULL);
408: ISDestroy(&tmp);
409: }
410: /* Compress indices */
411: PetscSectionSetUp(rootSectionAdj);
412: for (p = pStart; p < pEnd; ++p) {
413: PetscInt dof, cdof, off, d;
414: PetscInt adof, aoff;
416: PetscSectionGetDof(section, p, &dof);
417: PetscSectionGetConstraintDof(section, p, &cdof);
418: PetscSectionGetOffset(section, p, &off);
419: if (!dof) continue;
420: PetscSectionGetDof(rootSectionAdj, off, &adof);
421: if (adof <= 0) continue;
422: for (d = off; d < off + dof - cdof; ++d) {
423: PetscSectionGetDof(rootSectionAdj, d, &adof);
424: PetscSectionGetOffset(rootSectionAdj, d, &aoff);
425: PetscSortRemoveDupsInt(&adof, &rootAdj[aoff]);
426: PetscSectionSetDof(rootSectionAdj, d, adof);
427: }
428: }
429: /* Debugging */
430: if (debug) {
431: IS tmp;
432: PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after compression:\n");
433: PetscSectionView(rootSectionAdj, NULL);
434: PetscPrintf(comm, "Root adjacency indices after compression\n");
435: ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
436: ISView(tmp, NULL);
437: ISDestroy(&tmp);
438: }
439: /* Build adjacency section: Maps global indices to sets of adjacent global indices */
440: PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd);
441: PetscSectionCreate(comm, §ionAdj);
442: PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd);
443: for (p = pStart; p < pEnd; ++p) {
444: PetscInt numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof;
445: PetscBool found = PETSC_TRUE;
447: PetscSectionGetDof(section, p, &dof);
448: PetscSectionGetConstraintDof(section, p, &cdof);
449: PetscSectionGetOffset(section, p, &off);
450: PetscSectionGetOffset(sectionGlobal, p, &goff);
451: for (d = 0; d < dof - cdof; ++d) {
452: PetscInt ldof, rdof;
454: PetscSectionGetDof(leafSectionAdj, off + d, &ldof);
455: PetscSectionGetDof(rootSectionAdj, off + d, &rdof);
456: if (ldof > 0) {
457: /* We do not own this point */
458: } else if (rdof > 0) {
459: PetscSectionSetDof(sectionAdj, goff + d, rdof);
460: } else {
461: found = PETSC_FALSE;
462: }
463: }
464: if (found) continue;
465: PetscSectionGetDof(section, p, &dof);
466: PetscSectionGetOffset(sectionGlobal, p, &goff);
467: DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
468: for (q = 0; q < numAdj; ++q) {
469: const PetscInt padj = tmpAdj[q];
470: PetscInt ndof, ncdof, noff;
472: if ((padj < pStart) || (padj >= pEnd)) continue;
473: PetscSectionGetDof(section, padj, &ndof);
474: PetscSectionGetConstraintDof(section, padj, &ncdof);
475: PetscSectionGetOffset(section, padj, &noff);
476: for (d = goff; d < goff + dof - cdof; ++d) PetscSectionAddDof(sectionAdj, d, ndof - ncdof);
477: }
478: PetscSectionGetDof(anchorSectionAdj, p, &anDof);
479: if (anDof) {
480: for (d = goff; d < goff + dof - cdof; ++d) PetscSectionAddDof(sectionAdj, d, anDof);
481: }
482: }
483: PetscSectionSetUp(sectionAdj);
484: if (debug) {
485: PetscPrintf(comm, "Adjacency Section for Preallocation:\n");
486: PetscSectionView(sectionAdj, NULL);
487: }
488: /* Get adjacent indices */
489: PetscSectionGetStorageSize(sectionAdj, &numCols);
490: PetscMalloc1(numCols, &cols);
491: for (p = pStart; p < pEnd; ++p) {
492: PetscInt numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof, anOff;
493: PetscBool found = PETSC_TRUE;
495: PetscSectionGetDof(section, p, &dof);
496: PetscSectionGetConstraintDof(section, p, &cdof);
497: PetscSectionGetOffset(section, p, &off);
498: PetscSectionGetOffset(sectionGlobal, p, &goff);
499: for (d = 0; d < dof - cdof; ++d) {
500: PetscInt ldof, rdof;
502: PetscSectionGetDof(leafSectionAdj, off + d, &ldof);
503: PetscSectionGetDof(rootSectionAdj, off + d, &rdof);
504: if (ldof > 0) {
505: /* We do not own this point */
506: } else if (rdof > 0) {
507: PetscInt aoff, roff;
509: PetscSectionGetOffset(sectionAdj, goff + d, &aoff);
510: PetscSectionGetOffset(rootSectionAdj, off + d, &roff);
511: PetscArraycpy(&cols[aoff], &rootAdj[roff], rdof);
512: } else {
513: found = PETSC_FALSE;
514: }
515: }
516: if (found) continue;
517: DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
518: PetscSectionGetDof(anchorSectionAdj, p, &anDof);
519: PetscSectionGetOffset(anchorSectionAdj, p, &anOff);
520: for (d = goff; d < goff + dof - cdof; ++d) {
521: PetscInt adof, aoff, i = 0;
523: PetscSectionGetDof(sectionAdj, d, &adof);
524: PetscSectionGetOffset(sectionAdj, d, &aoff);
525: for (q = 0; q < numAdj; ++q) {
526: const PetscInt padj = tmpAdj[q];
527: PetscInt ndof, ncdof, ngoff, nd;
528: const PetscInt *ncind;
530: /* Adjacent points may not be in the section chart */
531: if ((padj < pStart) || (padj >= pEnd)) continue;
532: PetscSectionGetDof(section, padj, &ndof);
533: PetscSectionGetConstraintDof(section, padj, &ncdof);
534: PetscSectionGetConstraintIndices(section, padj, &ncind);
535: PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
536: for (nd = 0; nd < ndof - ncdof; ++nd, ++i) cols[aoff + i] = ngoff < 0 ? -(ngoff + 1) + nd : ngoff + nd;
537: }
538: for (q = 0; q < anDof; q++, i++) cols[aoff + i] = anchorAdj[anOff + q];
540: }
541: }
542: PetscSectionDestroy(&anchorSectionAdj);
543: PetscSectionDestroy(&leafSectionAdj);
544: PetscSectionDestroy(&rootSectionAdj);
545: PetscFree(anchorAdj);
546: PetscFree(rootAdj);
547: PetscFree(tmpAdj);
548: /* Debugging */
549: if (debug) {
550: IS tmp;
551: PetscPrintf(comm, "Column indices\n");
552: ISCreateGeneral(comm, numCols, cols, PETSC_USE_POINTER, &tmp);
553: ISView(tmp, NULL);
554: ISDestroy(&tmp);
555: }
557: *sA = sectionAdj;
558: *colIdx = cols;
559: return 0;
560: }
562: 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[])
563: {
564: PetscSection section;
565: PetscInt rStart, rEnd, r, pStart, pEnd, p;
567: /* This loop needs to change to a loop over points, then field dofs, which means we need to look both sections */
568: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
570: if (f >= 0 && bs == 1) {
571: DMGetLocalSection(dm, §ion);
572: PetscSectionGetChart(section, &pStart, &pEnd);
573: for (p = pStart; p < pEnd; ++p) {
574: PetscInt rS, rE;
576: DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE);
577: for (r = rS; r < rE; ++r) {
578: PetscInt numCols, cStart, c;
580: PetscSectionGetDof(sectionAdj, r, &numCols);
581: PetscSectionGetOffset(sectionAdj, r, &cStart);
582: for (c = cStart; c < cStart + numCols; ++c) {
583: if ((cols[c] >= rStart) && (cols[c] < rEnd)) {
584: ++dnz[r - rStart];
585: if (cols[c] >= r) ++dnzu[r - rStart];
586: } else {
587: ++onz[r - rStart];
588: if (cols[c] >= r) ++onzu[r - rStart];
589: }
590: }
591: }
592: }
593: } else {
594: /* Only loop over blocks of rows */
595: for (r = rStart / bs; r < rEnd / bs; ++r) {
596: const PetscInt row = r * bs;
597: PetscInt numCols, cStart, c;
599: PetscSectionGetDof(sectionAdj, row, &numCols);
600: PetscSectionGetOffset(sectionAdj, row, &cStart);
601: for (c = cStart; c < cStart + numCols; ++c) {
602: if ((cols[c] >= rStart) && (cols[c] < rEnd)) {
603: ++dnz[r - rStart / bs];
604: if (cols[c] >= row) ++dnzu[r - rStart / bs];
605: } else {
606: ++onz[r - rStart / bs];
607: if (cols[c] >= row) ++onzu[r - rStart / bs];
608: }
609: }
610: }
611: for (r = 0; r < (rEnd - rStart) / bs; ++r) {
612: dnz[r] /= bs;
613: onz[r] /= bs;
614: dnzu[r] /= bs;
615: onzu[r] /= bs;
616: }
617: }
618: return 0;
619: }
621: static PetscErrorCode DMPlexFillMatrix_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], Mat A)
622: {
623: PetscSection section;
624: PetscScalar *values;
625: PetscInt rStart, rEnd, r, pStart, pEnd, p, len, maxRowLen = 0;
627: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
628: for (r = rStart; r < rEnd; ++r) {
629: PetscSectionGetDof(sectionAdj, r, &len);
630: maxRowLen = PetscMax(maxRowLen, len);
631: }
632: PetscCalloc1(maxRowLen, &values);
633: if (f >= 0 && bs == 1) {
634: DMGetLocalSection(dm, §ion);
635: PetscSectionGetChart(section, &pStart, &pEnd);
636: for (p = pStart; p < pEnd; ++p) {
637: PetscInt rS, rE;
639: DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE);
640: for (r = rS; r < rE; ++r) {
641: PetscInt numCols, cStart;
643: PetscSectionGetDof(sectionAdj, r, &numCols);
644: PetscSectionGetOffset(sectionAdj, r, &cStart);
645: MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);
646: }
647: }
648: } else {
649: for (r = rStart; r < rEnd; ++r) {
650: PetscInt numCols, cStart;
652: PetscSectionGetDof(sectionAdj, r, &numCols);
653: PetscSectionGetOffset(sectionAdj, r, &cStart);
654: MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);
655: }
656: }
657: PetscFree(values);
658: return 0;
659: }
661: /*@C
662: DMPlexPreallocateOperator - Calculate the matrix nonzero pattern based upon the information in the `DM`,
663: the `PetscDS` it contains, and the default `PetscSection`.
665: Collective
667: Input Parameters:
668: + dm - The `DMPLEX`
669: . bs - The matrix blocksize
670: . dnz - An array to hold the number of nonzeros in the diagonal block
671: . onz - An array to hold the number of nonzeros in the off-diagonal block
672: . dnzu - An array to hold the number of nonzeros in the upper triangle of the diagonal block
673: . onzu - An array to hold the number of nonzeros in the upper triangle of the off-diagonal block
674: - fillMatrix - If `PETSC_TRUE`, fill the matrix with zeros
676: Output Parameter:
677: . A - The preallocated matrix
679: Level: advanced
681: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMCreateMatrix()`
682: @*/
683: PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
684: {
685: MPI_Comm comm;
686: PetscDS prob;
687: MatType mtype;
688: PetscSF sf, sfDof;
689: PetscSection section;
690: PetscInt *remoteOffsets;
691: PetscSection sectionAdj[4] = {NULL, NULL, NULL, NULL};
692: PetscInt *cols[4] = {NULL, NULL, NULL, NULL};
693: PetscBool useCone, useClosure;
694: PetscInt Nf, f, idx, locRows;
695: PetscLayout rLayout;
696: PetscBool isSymBlock, isSymSeqBlock, isSymMPIBlock, debug = PETSC_FALSE;
697: PetscMPIInt size;
705: DMGetDS(dm, &prob);
706: DMGetPointSF(dm, &sf);
707: DMGetLocalSection(dm, §ion);
708: PetscOptionsGetBool(NULL, NULL, "-dm_view_preallocation", &debug, NULL);
709: PetscObjectGetComm((PetscObject)dm, &comm);
710: MPI_Comm_size(comm, &size);
711: PetscLogEventBegin(DMPLEX_Preallocate, dm, 0, 0, 0);
712: /* Create dof SF based on point SF */
713: if (debug) {
714: PetscSection section, sectionGlobal;
715: PetscSF sf;
717: DMGetPointSF(dm, &sf);
718: DMGetLocalSection(dm, §ion);
719: DMGetGlobalSection(dm, §ionGlobal);
720: PetscPrintf(comm, "Input Section for Preallocation:\n");
721: PetscSectionView(section, NULL);
722: PetscPrintf(comm, "Input Global Section for Preallocation:\n");
723: PetscSectionView(sectionGlobal, NULL);
724: if (size > 1) {
725: PetscPrintf(comm, "Input SF for Preallocation:\n");
726: PetscSFView(sf, NULL);
727: }
728: }
729: PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets);
730: PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof);
731: PetscFree(remoteOffsets);
732: if (debug && size > 1) {
733: PetscPrintf(comm, "Dof SF for Preallocation:\n");
734: PetscSFView(sfDof, NULL);
735: }
736: /* Create allocation vectors from adjacency graph */
737: MatGetLocalSize(A, &locRows, NULL);
738: PetscLayoutCreate(comm, &rLayout);
739: PetscLayoutSetLocalSize(rLayout, locRows);
740: PetscLayoutSetBlockSize(rLayout, 1);
741: PetscLayoutSetUp(rLayout);
742: /* There are 4 types of adjacency */
743: PetscSectionGetNumFields(section, &Nf);
744: if (Nf < 1 || bs > 1) {
745: DMGetBasicAdjacency(dm, &useCone, &useClosure);
746: idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
747: DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, §ionAdj[idx], &cols[idx]);
748: DMPlexUpdateAllocation_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu);
749: } else {
750: for (f = 0; f < Nf; ++f) {
751: DMGetAdjacency(dm, f, &useCone, &useClosure);
752: idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
753: if (!sectionAdj[idx]) DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, §ionAdj[idx], &cols[idx]);
754: DMPlexUpdateAllocation_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu);
755: }
756: }
757: PetscSFDestroy(&sfDof);
758: /* Set matrix pattern */
759: MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu);
760: MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
761: /* Check for symmetric storage */
762: MatGetType(A, &mtype);
763: PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
764: PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
765: PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
766: if (isSymBlock || isSymSeqBlock || isSymMPIBlock) MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);
767: /* Fill matrix with zeros */
768: if (fillMatrix) {
769: if (Nf < 1 || bs > 1) {
770: DMGetBasicAdjacency(dm, &useCone, &useClosure);
771: idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
772: DMPlexFillMatrix_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], A);
773: } else {
774: for (f = 0; f < Nf; ++f) {
775: DMGetAdjacency(dm, f, &useCone, &useClosure);
776: idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
777: DMPlexFillMatrix_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], A);
778: }
779: }
780: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
781: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
782: }
783: PetscLayoutDestroy(&rLayout);
784: for (idx = 0; idx < 4; ++idx) {
785: PetscSectionDestroy(§ionAdj[idx]);
786: PetscFree(cols[idx]);
787: }
788: PetscLogEventEnd(DMPLEX_Preallocate, dm, 0, 0, 0);
789: return 0;
790: }
792: #if 0
793: PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
794: {
795: PetscInt *tmpClosure,*tmpAdj,*visits;
796: PetscInt c,cStart,cEnd,pStart,pEnd;
798: DMGetDimension(dm, &dim);
799: DMPlexGetDepth(dm, &depth);
800: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
802: maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1));
804: PetscSectionGetChart(section, &pStart, &pEnd);
805: npoints = pEnd - pStart;
807: PetscMalloc3(maxClosureSize,&tmpClosure,npoints,&lvisits,npoints,&visits);
808: PetscArrayzero(lvisits,pEnd-pStart);
809: PetscArrayzero(visits,pEnd-pStart);
810: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
811: for (c=cStart; c<cEnd; c++) {
812: PetscInt *support = tmpClosure;
813: DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support);
814: for (p=0; p<supportSize; p++) lvisits[support[p]]++;
815: }
816: PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM);
817: PetscSFReduceEnd (sf,MPIU_INT,lvisits,visits,MPI_SUM);
818: PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits,MPI_REPLACE);
819: PetscSFBcastEnd (sf,MPIU_INT,visits,lvisits);
821: PetscSFGetRootRanks();
823: PetscMalloc2(maxClosureSize*maxClosureSize,&cellmat,npoints,&owner);
824: for (c=cStart; c<cEnd; c++) {
825: PetscArrayzero(cellmat,maxClosureSize*maxClosureSize);
826: /*
827: Depth-first walk of transitive closure.
828: 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.
829: This contribution is added to dnz if owning ranks of p and q match, to onz otherwise.
830: */
831: }
833: PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM);
834: PetscSFReduceEnd (sf,MPIU_INT,lonz,onz,MPI_SUM);
835: return 0;
836: }
837: #endif