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