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