Actual source code: plexpreallocate.c
petsc-3.6.1 2015-08-06
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, "-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: MPI_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: if (debug) {
302: PetscPrintf(comm, "Adjacency SF for Preallocation:\n");
303: PetscSFView(sfAdj, NULL);
304: }
305: /* Create leaf adjacency */
306: PetscSectionSetUp(leafSectionAdj);
307: PetscSectionGetStorageSize(leafSectionAdj, &adjSize);
308: PetscCalloc1(adjSize, &adj);
309: for (l = 0; l < nleaves; ++l) {
310: PetscInt dof, off, d, q, anDof, anOff;
311: PetscInt p = leaves[l], numAdj = PETSC_DETERMINE;
313: if ((p < pStart) || (p >= pEnd)) continue;
314: PetscSectionGetDof(section, p, &dof);
315: PetscSectionGetOffset(section, p, &off);
316: DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
317: PetscSectionGetDof(anchorSectionAdj, p, &anDof);
318: PetscSectionGetOffset(anchorSectionAdj, p, &anOff);
319: for (d = off; d < off+dof; ++d) {
320: PetscInt aoff, i = 0;
322: PetscSectionGetOffset(leafSectionAdj, d, &aoff);
323: for (q = 0; q < numAdj; ++q) {
324: const PetscInt padj = tmpAdj[q];
325: PetscInt ndof, ncdof, ngoff, nd;
327: if ((padj < pStart) || (padj >= pEnd)) continue;
328: PetscSectionGetDof(section, padj, &ndof);
329: PetscSectionGetConstraintDof(section, padj, &ncdof);
330: PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
331: for (nd = 0; nd < ndof-ncdof; ++nd) {
332: adj[aoff+i] = (ngoff < 0 ? -(ngoff+1) : ngoff) + nd;
333: ++i;
334: }
335: }
336: for (q = 0; q < anDof; q++) {
337: adj[aoff+i] = anchorAdj[anOff+q];
338: ++i;
339: }
340: }
341: }
342: /* Debugging */
343: if (debug) {
344: IS tmp;
345: PetscPrintf(comm, "Leaf adjacency indices\n");
346: ISCreateGeneral(comm, adjSize, adj, PETSC_USE_POINTER, &tmp);
347: ISView(tmp, NULL);
348: ISDestroy(&tmp);
349: }
350: /* Gather adjacent indices to root */
351: PetscSectionGetStorageSize(rootSectionAdj, &adjSize);
352: PetscMalloc1(adjSize, &rootAdj);
353: for (r = 0; r < adjSize; ++r) rootAdj[r] = -1;
354: if (doComm) {
355: const PetscInt *indegree;
356: PetscInt *remoteadj, radjsize = 0;
358: PetscSFComputeDegreeBegin(sfAdj, &indegree);
359: PetscSFComputeDegreeEnd(sfAdj, &indegree);
360: for (p = 0; p < adjSize; ++p) radjsize += indegree[p];
361: PetscMalloc1(radjsize, &remoteadj);
362: PetscSFGatherBegin(sfAdj, MPIU_INT, adj, remoteadj);
363: PetscSFGatherEnd(sfAdj, MPIU_INT, adj, remoteadj);
364: for (p = 0, l = 0, r = 0; p < adjSize; ++p, l = PetscMax(p, l + indegree[p-1])) {
365: PetscInt s;
366: for (s = 0; s < indegree[p]; ++s, ++r) rootAdj[l+s] = remoteadj[r];
367: }
368: if (r != radjsize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %d != %d", r, radjsize);
369: if (l != adjSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %d != %d", l, adjSize);
370: PetscFree(remoteadj);
371: }
372: PetscSFDestroy(&sfAdj);
373: PetscFree(adj);
374: /* Debugging */
375: if (debug) {
376: IS tmp;
377: PetscPrintf(comm, "Root adjacency indices after gather\n");
378: ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
379: ISView(tmp, NULL);
380: ISDestroy(&tmp);
381: }
382: /* Add in local adjacency indices for owned dofs on interface (roots) */
383: for (p = pStart; p < pEnd; ++p) {
384: PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof, anOff;
386: PetscSectionGetDof(section, p, &dof);
387: PetscSectionGetOffset(section, p, &off);
388: if (!dof) continue;
389: PetscSectionGetDof(rootSectionAdj, off, &adof);
390: if (adof <= 0) continue;
391: DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
392: PetscSectionGetDof(anchorSectionAdj, p, &anDof);
393: PetscSectionGetOffset(anchorSectionAdj, p, &anOff);
394: for (d = off; d < off+dof; ++d) {
395: PetscInt adof, aoff, i;
397: PetscSectionGetDof(rootSectionAdj, d, &adof);
398: PetscSectionGetOffset(rootSectionAdj, d, &aoff);
399: i = adof-1;
400: for (q = 0; q < anDof; q++) {
401: rootAdj[aoff+i] = anchorAdj[anOff+q];
402: --i;
403: }
404: for (q = 0; q < numAdj; ++q) {
405: const PetscInt padj = tmpAdj[q];
406: PetscInt ndof, ncdof, ngoff, nd;
408: if ((padj < pStart) || (padj >= pEnd)) continue;
409: PetscSectionGetDof(section, padj, &ndof);
410: PetscSectionGetConstraintDof(section, padj, &ncdof);
411: PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
412: for (nd = 0; nd < ndof-ncdof; ++nd) {
413: rootAdj[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
414: --i;
415: }
416: }
417: }
418: }
419: /* Debugging */
420: if (debug) {
421: IS tmp;
422: PetscPrintf(comm, "Root adjacency indices\n");
423: ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
424: ISView(tmp, NULL);
425: ISDestroy(&tmp);
426: }
427: /* Compress indices */
428: PetscSectionSetUp(rootSectionAdj);
429: for (p = pStart; p < pEnd; ++p) {
430: PetscInt dof, cdof, off, d;
431: PetscInt adof, aoff;
433: PetscSectionGetDof(section, p, &dof);
434: PetscSectionGetConstraintDof(section, p, &cdof);
435: PetscSectionGetOffset(section, p, &off);
436: if (!dof) continue;
437: PetscSectionGetDof(rootSectionAdj, off, &adof);
438: if (adof <= 0) continue;
439: for (d = off; d < off+dof-cdof; ++d) {
440: PetscSectionGetDof(rootSectionAdj, d, &adof);
441: PetscSectionGetOffset(rootSectionAdj, d, &aoff);
442: PetscSortRemoveDupsInt(&adof, &rootAdj[aoff]);
443: PetscSectionSetDof(rootSectionAdj, d, adof);
444: }
445: }
446: /* Debugging */
447: if (debug) {
448: IS tmp;
449: PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after compression:\n");
450: PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);
451: PetscPrintf(comm, "Root adjacency indices after compression\n");
452: ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
453: ISView(tmp, NULL);
454: ISDestroy(&tmp);
455: }
456: /* Build adjacency section: Maps global indices to sets of adjacent global indices */
457: PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd);
458: PetscSectionCreate(comm, §ionAdj);
459: PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd);
460: for (p = pStart; p < pEnd; ++p) {
461: PetscInt numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof;
462: PetscBool found = PETSC_TRUE;
464: PetscSectionGetDof(section, p, &dof);
465: PetscSectionGetConstraintDof(section, p, &cdof);
466: PetscSectionGetOffset(section, p, &off);
467: PetscSectionGetOffset(sectionGlobal, p, &goff);
468: for (d = 0; d < dof-cdof; ++d) {
469: PetscInt ldof, rdof;
471: PetscSectionGetDof(leafSectionAdj, off+d, &ldof);
472: PetscSectionGetDof(rootSectionAdj, off+d, &rdof);
473: if (ldof > 0) {
474: /* We do not own this point */
475: } else if (rdof > 0) {
476: PetscSectionSetDof(sectionAdj, goff+d, rdof);
477: } else {
478: found = PETSC_FALSE;
479: }
480: }
481: if (found) continue;
482: PetscSectionGetDof(section, p, &dof);
483: PetscSectionGetOffset(sectionGlobal, p, &goff);
484: DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
485: for (q = 0; q < numAdj; ++q) {
486: const PetscInt padj = tmpAdj[q];
487: PetscInt ndof, ncdof, noff;
489: if ((padj < pStart) || (padj >= pEnd)) continue;
490: PetscSectionGetDof(section, padj, &ndof);
491: PetscSectionGetConstraintDof(section, padj, &ncdof);
492: PetscSectionGetOffset(section, padj, &noff);
493: for (d = goff; d < goff+dof-cdof; ++d) {
494: PetscSectionAddDof(sectionAdj, d, ndof-ncdof);
495: }
496: }
497: PetscSectionGetDof(anchorSectionAdj, p, &anDof);
498: if (anDof) {
499: for (d = goff; d < goff+dof-cdof; ++d) {
500: PetscSectionAddDof(sectionAdj, d, anDof);
501: }
502: }
503: }
504: PetscSectionSetUp(sectionAdj);
505: if (debug) {
506: PetscPrintf(comm, "Adjacency Section for Preallocation:\n");
507: PetscSectionView(sectionAdj, PETSC_VIEWER_STDOUT_WORLD);
508: }
509: /* Get adjacent indices */
510: PetscSectionGetStorageSize(sectionAdj, &numCols);
511: PetscMalloc1(numCols, &cols);
512: for (p = pStart; p < pEnd; ++p) {
513: PetscInt numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof, anOff;
514: PetscBool found = PETSC_TRUE;
516: PetscSectionGetDof(section, p, &dof);
517: PetscSectionGetConstraintDof(section, p, &cdof);
518: PetscSectionGetOffset(section, p, &off);
519: PetscSectionGetOffset(sectionGlobal, p, &goff);
520: for (d = 0; d < dof-cdof; ++d) {
521: PetscInt ldof, rdof;
523: PetscSectionGetDof(leafSectionAdj, off+d, &ldof);
524: PetscSectionGetDof(rootSectionAdj, off+d, &rdof);
525: if (ldof > 0) {
526: /* We do not own this point */
527: } else if (rdof > 0) {
528: PetscInt aoff, roff;
530: PetscSectionGetOffset(sectionAdj, goff+d, &aoff);
531: PetscSectionGetOffset(rootSectionAdj, off+d, &roff);
532: PetscMemcpy(&cols[aoff], &rootAdj[roff], rdof * sizeof(PetscInt));
533: } else {
534: found = PETSC_FALSE;
535: }
536: }
537: if (found) continue;
538: DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
539: PetscSectionGetDof(anchorSectionAdj, p, &anDof);
540: PetscSectionGetOffset(anchorSectionAdj, p, &anOff);
541: for (d = goff; d < goff+dof-cdof; ++d) {
542: PetscInt adof, aoff, i = 0;
544: PetscSectionGetDof(sectionAdj, d, &adof);
545: PetscSectionGetOffset(sectionAdj, d, &aoff);
546: for (q = 0; q < numAdj; ++q) {
547: const PetscInt padj = tmpAdj[q];
548: PetscInt ndof, ncdof, ngoff, nd;
549: const PetscInt *ncind;
551: /* Adjacent points may not be in the section chart */
552: if ((padj < pStart) || (padj >= pEnd)) continue;
553: PetscSectionGetDof(section, padj, &ndof);
554: PetscSectionGetConstraintDof(section, padj, &ncdof);
555: PetscSectionGetConstraintIndices(section, padj, &ncind);
556: PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
557: for (nd = 0; nd < ndof-ncdof; ++nd, ++i) {
558: cols[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
559: }
560: }
561: for (q = 0; q < anDof; q++, i++) {
562: cols[aoff+i] = anchorAdj[anOff + q];
563: }
564: 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);
565: }
566: }
567: PetscSectionDestroy(&anchorSectionAdj);
568: PetscSectionDestroy(&leafSectionAdj);
569: PetscSectionDestroy(&rootSectionAdj);
570: PetscFree(anchorAdj);
571: PetscFree(rootAdj);
572: PetscFree(tmpAdj);
573: /* Debugging */
574: if (debug) {
575: IS tmp;
576: PetscPrintf(comm, "Column indices\n");
577: ISCreateGeneral(comm, numCols, cols, PETSC_USE_POINTER, &tmp);
578: ISView(tmp, NULL);
579: ISDestroy(&tmp);
580: }
582: *sA = sectionAdj;
583: *colIdx = cols;
584: return(0);
585: }
589: PetscErrorCode DMPlexUpdateAllocation_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[])
590: {
591: PetscSection section;
592: PetscInt rStart, rEnd, r, pStart, pEnd, p;
596: /* This loop needs to change to a loop over points, then field dofs, which means we need to look both sections */
597: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
598: 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);
599: if (f >= 0 && bs == 1) {
600: DMGetDefaultSection(dm, §ion);
601: PetscSectionGetChart(section, &pStart, &pEnd);
602: for (p = pStart; p < pEnd; ++p) {
603: PetscInt rS, rE;
605: DMPlexGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE);
606: for (r = rS; r < rE; ++r) {
607: PetscInt numCols, cStart, c;
609: PetscSectionGetDof(sectionAdj, r, &numCols);
610: PetscSectionGetOffset(sectionAdj, r, &cStart);
611: for (c = cStart; c < cStart+numCols; ++c) {
612: if ((cols[c] >= rStart) && (cols[c] < rEnd)) {
613: ++dnz[r-rStart];
614: if (cols[c] >= r) ++dnzu[r-rStart];
615: } else {
616: ++onz[r-rStart];
617: if (cols[c] >= r) ++onzu[r-rStart];
618: }
619: }
620: }
621: }
622: } else {
623: /* Only loop over blocks of rows */
624: for (r = rStart/bs; r < rEnd/bs; ++r) {
625: const PetscInt row = r*bs;
626: PetscInt numCols, cStart, c;
628: PetscSectionGetDof(sectionAdj, row, &numCols);
629: PetscSectionGetOffset(sectionAdj, row, &cStart);
630: for (c = cStart; c < cStart+numCols; ++c) {
631: if ((cols[c] >= rStart*bs) && (cols[c] < rEnd*bs)) {
632: ++dnz[r-rStart];
633: if (cols[c] >= row) ++dnzu[r-rStart];
634: } else {
635: ++onz[r-rStart];
636: if (cols[c] >= row) ++onzu[r-rStart];
637: }
638: }
639: }
640: for (r = 0; r < (rEnd - rStart)/bs; ++r) {
641: dnz[r] /= bs;
642: onz[r] /= bs;
643: dnzu[r] /= bs;
644: onzu[r] /= bs;
645: }
646: }
647: return(0);
648: }
652: PetscErrorCode DMPlexFillMatrix_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], Mat A)
653: {
654: PetscSection section;
655: PetscScalar *values;
656: PetscInt rStart, rEnd, r, pStart, pEnd, p, len, maxRowLen = 0;
660: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
661: for (r = rStart; r < rEnd; ++r) {
662: PetscSectionGetDof(sectionAdj, r, &len);
663: maxRowLen = PetscMax(maxRowLen, len);
664: }
665: PetscCalloc1(maxRowLen, &values);
666: if (f >=0 && bs == 1) {
667: DMGetDefaultSection(dm, §ion);
668: PetscSectionGetChart(section, &pStart, &pEnd);
669: for (p = pStart; p < pEnd; ++p) {
670: PetscInt rS, rE;
672: DMPlexGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE);
673: for (r = rS; r < rE; ++r) {
674: PetscInt numCols, cStart;
676: PetscSectionGetDof(sectionAdj, r, &numCols);
677: PetscSectionGetOffset(sectionAdj, r, &cStart);
678: MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);
679: }
680: }
681: } else {
682: for (r = rStart; r < rEnd; ++r) {
683: PetscInt numCols, cStart;
685: PetscSectionGetDof(sectionAdj, r, &numCols);
686: PetscSectionGetOffset(sectionAdj, r, &cStart);
687: MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);
688: }
689: }
690: PetscFree(values);
691: return(0);
692: }
696: PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
697: {
698: MPI_Comm comm;
699: PetscDS prob;
700: MatType mtype;
701: PetscSF sf, sfDof;
702: PetscSection section;
703: PetscInt *remoteOffsets;
704: PetscSection sectionAdj[4] = {NULL, NULL, NULL, NULL};
705: PetscInt *cols[4] = {NULL, NULL, NULL, NULL};
706: PetscBool useCone, useClosure;
707: PetscInt Nf, f, idx, locRows;
708: PetscLayout rLayout;
709: PetscBool isSymBlock, isSymSeqBlock, isSymMPIBlock, debug = PETSC_FALSE;
717: DMGetDS(dm, &prob);
718: DMGetPointSF(dm, &sf);
719: DMGetDefaultSection(dm, §ion);
720: PetscOptionsGetBool(NULL, "-dm_view_preallocation", &debug, NULL);
721: PetscObjectGetComm((PetscObject) dm, &comm);
722: PetscLogEventBegin(DMPLEX_Preallocate,dm,0,0,0);
723: /* Create dof SF based on point SF */
724: if (debug) {
725: PetscSection section, sectionGlobal;
726: PetscSF sf;
728: DMGetPointSF(dm, &sf);
729: DMGetDefaultSection(dm, §ion);
730: DMGetDefaultGlobalSection(dm, §ionGlobal);
731: PetscPrintf(comm, "Input Section for Preallocation:\n");
732: PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD);
733: PetscPrintf(comm, "Input Global Section for Preallocation:\n");
734: PetscSectionView(sectionGlobal, PETSC_VIEWER_STDOUT_WORLD);
735: PetscPrintf(comm, "Input SF for Preallocation:\n");
736: PetscSFView(sf, NULL);
737: }
738: PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets);
739: PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof);
740: if (debug) {
741: PetscPrintf(comm, "Dof SF for Preallocation:\n");
742: PetscSFView(sfDof, NULL);
743: }
744: /* Create allocation vectors from adjacency graph */
745: MatGetLocalSize(A, &locRows, NULL);
746: PetscLayoutCreate(comm, &rLayout);
747: PetscLayoutSetLocalSize(rLayout, locRows);
748: PetscLayoutSetBlockSize(rLayout, 1);
749: PetscLayoutSetUp(rLayout);
750: /* There are 4 types of adjacency */
751: PetscSectionGetNumFields(section, &Nf);
752: if (Nf < 1 || bs > 1) {
753: DMPlexGetAdjacencyUseCone(dm, &useCone);
754: DMPlexGetAdjacencyUseClosure(dm, &useClosure);
755: idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
756: DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, §ionAdj[idx], &cols[idx]);
757: DMPlexUpdateAllocation_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu);
758: } else {
759: for (f = 0; f < Nf; ++f) {
760: PetscDSGetAdjacency(prob, f, &useCone, &useClosure);
761: idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
762: if (!sectionAdj[idx]) {DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, §ionAdj[idx], &cols[idx]);}
763: DMPlexUpdateAllocation_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu);
764: }
765: }
766: PetscSFDestroy(&sfDof);
767: /* Set matrix pattern */
768: MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu);
769: MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
770: /* Check for symmetric storage */
771: MatGetType(A, &mtype);
772: PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
773: PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
774: PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
775: if (isSymBlock || isSymSeqBlock || isSymMPIBlock) {MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);}
776: /* Fill matrix with zeros */
777: if (fillMatrix) {
778: if (Nf < 1 || bs > 1) {
779: DMPlexGetAdjacencyUseCone(dm, &useCone);
780: DMPlexGetAdjacencyUseClosure(dm, &useClosure);
781: idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
782: DMPlexFillMatrix_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], A);
783: } else {
784: for (f = 0; f < Nf; ++f) {
785: PetscDSGetAdjacency(prob, f, &useCone, &useClosure);
786: idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
787: DMPlexFillMatrix_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], A);
788: }
789: }
790: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
791: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
792: }
793: PetscLayoutDestroy(&rLayout);
794: for (idx = 0; idx < 4; ++idx) {PetscSectionDestroy(§ionAdj[idx]); PetscFree(cols[idx]);}
795: PetscLogEventEnd(DMPLEX_Preallocate,dm,0,0,0);
796: return(0);
797: }
799: #if 0
802: PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
803: {
804: PetscInt *tmpClosure,*tmpAdj,*visits;
805: PetscInt c,cStart,cEnd,pStart,pEnd;
806: PetscErrorCode ierr;
809: DMGetDimension(dm, &dim);
810: DMPlexGetDepth(dm, &depth);
811: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
813: maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1));
815: PetscSectionGetChart(section, &pStart, &pEnd);
816: npoints = pEnd - pStart;
818: PetscMalloc3(maxClosureSize,&tmpClosure,npoints,&lvisits,npoints,&visits);
819: PetscMemzero(lvisits,(pEnd-pStart)*sizeof(PetscInt));
820: PetscMemzero(visits,(pEnd-pStart)*sizeof(PetscInt));
821: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
822: for (c=cStart; c<cEnd; c++) {
823: PetscInt *support = tmpClosure;
824: DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support);
825: for (p=0; p<supportSize; p++) lvisits[support[p]]++;
826: }
827: PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM);
828: PetscSFReduceEnd (sf,MPIU_INT,lvisits,visits,MPI_SUM);
829: PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits);
830: PetscSFBcastEnd (sf,MPIU_INT,visits,lvisits);
832: PetscSFGetRanks();
835: PetscMalloc2(maxClosureSize*maxClosureSize,&cellmat,npoints,&owner);
836: for (c=cStart; c<cEnd; c++) {
837: PetscMemzero(cellmat,maxClosureSize*maxClosureSize*sizeof(PetscInt));
838: /*
839: Depth-first walk of transitive closure.
840: 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.
841: This contribution is added to dnz if owning ranks of p and q match, to onz otherwise.
842: */
843: }
845: PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM);
846: PetscSFReduceEnd (sf,MPIU_INT,lonz,onz,MPI_SUM);
847: return(0);
848: }
849: #endif