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