xref: /petsc/src/dm/impls/plex/plexpreallocate.c (revision 7c8652dd9fb051dfaf30896d504f41ba028df3ea) ! 
1 #include <petsc-private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
2 #include <petsc-private/isimpl.h>
3 #include <petscsf.h>
4 
5 #undef __FUNCT__
6 #define __FUNCT__ "DMPlexGetAdjacency_Internal"
7 static PetscErrorCode DMPlexGetAdjacency_Internal(DM dm, PetscInt p, PetscBool useClosure, const PetscInt *tmpClosure, PetscInt *adjSize, PetscInt adj[])
8 {
9   const PetscInt *star  = tmpClosure;
10   PetscInt        numAdj = 0, maxAdjSize = *adjSize, starSize, s;
11   PetscErrorCode  ierr;
12 
13   PetscFunctionBegin;
14   ierr = DMPlexGetTransitiveClosure(dm, p, useClosure, &starSize, (PetscInt**) &star);CHKERRQ(ierr);
15   for (s = 2; s < starSize*2; s += 2) {
16     const PetscInt *closure = NULL;
17     PetscInt        closureSize, c, q;
18 
19     ierr = DMPlexGetTransitiveClosure(dm, star[s], (PetscBool)!useClosure, &closureSize, (PetscInt**) &closure);CHKERRQ(ierr);
20     for (c = 0; c < closureSize*2; c += 2) {
21       for (q = 0; q < numAdj || (adj[numAdj++] = closure[c],0); ++q) {
22         if (closure[c] == adj[q]) break;
23       }
24       if (numAdj > maxAdjSize) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid mesh exceeded adjacency allocation (%D)", maxAdjSize);
25     }
26     ierr = DMPlexRestoreTransitiveClosure(dm, star[s], (PetscBool)!useClosure, &closureSize, (PetscInt**) &closure);CHKERRQ(ierr);
27   }
28   *adjSize = numAdj;
29   PetscFunctionReturn(0);
30 }
31 
32 #undef __FUNCT__
33 #define __FUNCT__ "DMPlexSetPreallocationCenterDimension"
34 /*@
35   DMPlexSetPreallocationCenterDimension - Determine the topology used to determine adjacency
36 
37   Input Parameters:
38 + dm - The DM object
39 - preallocCenterDim - The dimension of points which connect adjacent entries
40 
41   Level: developer
42 
43   Notes:
44 $     FEM:   Two points p and q are adjacent if q \in closure(star(p)), preallocCenterDim = dim
45 $     FVM:   Two points p and q are adjacent if q \in star(cone(p)),    preallocCenterDim = dim-1
46 $     FVM++: Two points p and q are adjacent if q \in star(closure(p)), preallocCenterDim = 0
47 
48 .seealso: DMCreateMatrix(), DMPlexPreallocateOperator(), DMPlexGetPreallocationCenterDimension()
49 @*/
50 PetscErrorCode DMPlexSetPreallocationCenterDimension(DM dm, PetscInt preallocCenterDim)
51 {
52   DM_Plex *mesh = (DM_Plex*) dm->data;
53 
54   PetscFunctionBegin;
55   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
56   mesh->preallocCenterDim = preallocCenterDim;
57   PetscFunctionReturn(0);
58 }
59 
60 #undef __FUNCT__
61 #define __FUNCT__ "DMPlexGetPreallocationCenterDimension"
62 /*@
63   DMPlexGetPreallocationCenterDimension - Return the topology used to determine adjacency
64 
65   Input Parameter:
66 . dm - The DM object
67 
68   Output Parameter:
69 . preallocCenterDim - The dimension of points which connect adjacent entries
70 
71   Level: developer
72 
73   Notes:
74 $     FEM:   Two points p and q are adjacent if q \in closure(star(p)), preallocCenterDim = dim
75 $     FVM:   Two points p and q are adjacent if q \in star(cone(p)),    preallocCenterDim = dim-1
76 $     FVM++: Two points p and q are adjacent if q \in star(closure(p)), preallocCenterDim = 0
77 
78 .seealso: DMCreateMatrix(), DMPlexPreallocateOperator(), DMPlexSetPreallocationCenterDimension()
79 @*/
80 PetscErrorCode DMPlexGetPreallocationCenterDimension(DM dm, PetscInt *preallocCenterDim)
81 {
82   DM_Plex *mesh = (DM_Plex*) dm->data;
83 
84   PetscFunctionBegin;
85   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
86   PetscValidIntPointer(preallocCenterDim, 2);
87   *preallocCenterDim = mesh->preallocCenterDim;
88   PetscFunctionReturn(0);
89 }
90 
91 #undef __FUNCT__
92 #define __FUNCT__ "DMPlexPreallocateOperator"
93 PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
94 {
95   MPI_Comm           comm;
96   MatType            mtype;
97   PetscSF            sf, sfDof, sfAdj;
98   PetscSection       leafSectionAdj, rootSectionAdj, sectionAdj;
99   PetscInt           nroots, nleaves, l, p;
100   const PetscInt    *leaves;
101   const PetscSFNode *remotes;
102   PetscInt           dim, pStart, pEnd, numDof, globalOffStart, globalOffEnd, numCols;
103   PetscInt          *tmpClosure, *tmpAdj, *adj, *rootAdj, *cols, *remoteOffsets;
104   PetscInt           depth, centerDim, maxConeSize, maxSupportSize, maxClosureSize, maxAdjSize, adjSize;
105   PetscLayout        rLayout;
106   PetscInt           locRows, rStart, rEnd, r;
107   PetscMPIInt        size;
108   PetscBool          useClosure, doCommLocal, doComm, debug = PETSC_FALSE, isSymBlock, isSymSeqBlock, isSymMPIBlock;
109   PetscErrorCode     ierr;
110 
111   PetscFunctionBegin;
112   ierr = PetscLogEventBegin(DMPLEX_Preallocate,dm,0,0,0);CHKERRQ(ierr);
113   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
114   ierr = PetscOptionsGetBool(NULL, "-dm_view_preallocation", &debug, NULL);CHKERRQ(ierr);
115   ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
116   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
117   ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
118   ierr = PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);CHKERRQ(ierr);
119   doCommLocal = (size > 1) && (nroots >= 0) ? PETSC_TRUE : PETSC_FALSE;
120   ierr = MPI_Allreduce(&doCommLocal, &doComm, 1, MPIU_BOOL, MPI_LAND, comm);CHKERRQ(ierr);
121   /* Create dof SF based on point SF */
122   if (debug) {
123     ierr = PetscPrintf(comm, "Input Section for Preallocation:\n");CHKERRQ(ierr);
124     ierr = PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
125     ierr = PetscPrintf(comm, "Input Global Section for Preallocation:\n");CHKERRQ(ierr);
126     ierr = PetscSectionView(sectionGlobal, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
127     ierr = PetscPrintf(comm, "Input SF for Preallocation:\n");CHKERRQ(ierr);
128     ierr = PetscSFView(sf, NULL);CHKERRQ(ierr);
129   }
130   ierr = PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets);CHKERRQ(ierr);
131   ierr = PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof);CHKERRQ(ierr);
132   if (debug) {
133     ierr = PetscPrintf(comm, "Dof SF for Preallocation:\n");CHKERRQ(ierr);
134     ierr = PetscSFView(sfDof, NULL);CHKERRQ(ierr);
135   }
136   /* Create section for dof adjacency (dof ==> # adj dof) */
137   /*   FEM:   Two points p and q are adjacent if q \in closure(star(p)), preallocCenterDim = dim   */
138   /*   FVM:   Two points p and q are adjacent if q \in star(cone(p)),    preallocCenterDim = dim-1 */
139   /*   FVM++: Two points p and q are adjacent if q \in star(closure(p)), preallocCenterDim = 0     */
140   ierr = DMPlexGetPreallocationCenterDimension(dm, &centerDim);CHKERRQ(ierr);
141   if (centerDim == dim) {
142     useClosure = PETSC_FALSE;
143   } else if (centerDim == 0) {
144     useClosure = PETSC_TRUE;
145   } else SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Do not support preallocation with center points of dimension %d", centerDim);
146 
147   ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
148   ierr = PetscSectionGetStorageSize(section, &numDof);CHKERRQ(ierr);
149   ierr = PetscSectionCreate(comm, &leafSectionAdj);CHKERRQ(ierr);
150   ierr = PetscSectionSetChart(leafSectionAdj, 0, numDof);CHKERRQ(ierr);
151   ierr = PetscSectionCreate(comm, &rootSectionAdj);CHKERRQ(ierr);
152   ierr = PetscSectionSetChart(rootSectionAdj, 0, numDof);CHKERRQ(ierr);
153   /*   Fill in the ghost dofs on the interface */
154   ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, &remotes);CHKERRQ(ierr);
155   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
156   ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
157 
158   maxClosureSize = 2*PetscMax(PetscPowInt(maxConeSize,depth+1),PetscPowInt(maxSupportSize,depth+1));
159   maxAdjSize     = PetscPowInt(maxConeSize,depth+1) * PetscPowInt(maxSupportSize,depth+1) + 1;
160 
161   ierr = PetscMalloc2(maxClosureSize,&tmpClosure,maxAdjSize,&tmpAdj);CHKERRQ(ierr);
162 
163   /*
164    ** The bootstrapping process involves six rounds with similar structure of visiting neighbors of each point.
165     1. Visit unowned points on interface, count adjacencies placing in leafSectionAdj
166        Reduce those counts to rootSectionAdj (now redundantly counting some interface points)
167     2. Visit owned points on interface, count adjacencies placing in rootSectionAdj
168        Create sfAdj connecting rootSectionAdj and leafSectionAdj
169     3. Visit unowned points on interface, write adjacencies to adj
170        Gather adj to rootAdj (note that there is redundancy in rootAdj when multiple procs find the same adjacencies)
171     4. Visit owned points on interface, write adjacencies to rootAdj
172        Remove redundancy in rootAdj
173    ** The last two traversals use transitive closure
174     5. Visit all owned points in the subdomain, count dofs for each point (sectionAdj)
175        Allocate memory addressed by sectionAdj (cols)
176     6. Visit all owned points in the subdomain, insert dof adjacencies into cols
177    ** Knowing all the column adjacencies, check ownership and sum into dnz and onz
178   */
179 
180   for (l = 0; l < nleaves; ++l) {
181     PetscInt dof, off, d, q;
182     PetscInt p = leaves[l], numAdj = maxAdjSize;
183 
184     if ((p < pStart) || (p >= pEnd)) continue;
185     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
186     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
187     ierr = DMPlexGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);CHKERRQ(ierr);
188     for (q = 0; q < numAdj; ++q) {
189       const PetscInt padj = tmpAdj[q];
190       PetscInt ndof, ncdof;
191 
192       if ((padj < pStart) || (padj >= pEnd)) continue;
193       ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
194       ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
195       for (d = off; d < off+dof; ++d) {
196         ierr = PetscSectionAddDof(leafSectionAdj, d, ndof-ncdof);CHKERRQ(ierr);
197       }
198     }
199   }
200   ierr = PetscSectionSetUp(leafSectionAdj);CHKERRQ(ierr);
201   if (debug) {
202     ierr = PetscPrintf(comm, "Adjacency Section for Preallocation on Leaves:\n");CHKERRQ(ierr);
203     ierr = PetscSectionView(leafSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
204   }
205   /* Get maximum remote adjacency sizes for owned dofs on interface (roots) */
206   if (doComm) {
207     ierr = PetscSFReduceBegin(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);CHKERRQ(ierr);
208     ierr = PetscSFReduceEnd(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);CHKERRQ(ierr);
209   }
210   if (debug) {
211     ierr = PetscPrintf(comm, "Adjancency Section for Preallocation on Roots:\n");CHKERRQ(ierr);
212     ierr = PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
213   }
214   /* Add in local adjacency sizes for owned dofs on interface (roots) */
215   for (p = pStart; p < pEnd; ++p) {
216     PetscInt numAdj = maxAdjSize, adof, dof, off, d, q;
217 
218     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
219     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
220     if (!dof) continue;
221     ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
222     if (adof <= 0) continue;
223     ierr = DMPlexGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);CHKERRQ(ierr);
224     for (q = 0; q < numAdj; ++q) {
225       const PetscInt padj = tmpAdj[q];
226       PetscInt ndof, ncdof;
227 
228       if ((padj < pStart) || (padj >= pEnd)) continue;
229       ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
230       ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
231       for (d = off; d < off+dof; ++d) {
232         ierr = PetscSectionAddDof(rootSectionAdj, d, ndof-ncdof);CHKERRQ(ierr);
233       }
234     }
235   }
236   ierr = PetscSectionSetUp(rootSectionAdj);CHKERRQ(ierr);
237   if (debug) {
238     ierr = PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after local additions:\n");CHKERRQ(ierr);
239     ierr = PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
240   }
241   /* Create adj SF based on dof SF */
242   ierr = PetscSFCreateRemoteOffsets(sfDof, rootSectionAdj, leafSectionAdj, &remoteOffsets);CHKERRQ(ierr);
243   ierr = PetscSFCreateSectionSF(sfDof, rootSectionAdj, remoteOffsets, leafSectionAdj, &sfAdj);CHKERRQ(ierr);
244   if (debug) {
245     ierr = PetscPrintf(comm, "Adjacency SF for Preallocation:\n");CHKERRQ(ierr);
246     ierr = PetscSFView(sfAdj, NULL);CHKERRQ(ierr);
247   }
248   ierr = PetscSFDestroy(&sfDof);CHKERRQ(ierr);
249   /* Create leaf adjacency */
250   ierr = PetscSectionSetUp(leafSectionAdj);CHKERRQ(ierr);
251   ierr = PetscSectionGetStorageSize(leafSectionAdj, &adjSize);CHKERRQ(ierr);
252   ierr = PetscCalloc1(adjSize, &adj);CHKERRQ(ierr);
253   for (l = 0; l < nleaves; ++l) {
254     PetscInt dof, off, d, q;
255     PetscInt p = leaves[l], numAdj = maxAdjSize;
256 
257     if ((p < pStart) || (p >= pEnd)) continue;
258     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
259     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
260     ierr = DMPlexGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);CHKERRQ(ierr);
261     for (d = off; d < off+dof; ++d) {
262       PetscInt aoff, i = 0;
263 
264       ierr = PetscSectionGetOffset(leafSectionAdj, d, &aoff);CHKERRQ(ierr);
265       for (q = 0; q < numAdj; ++q) {
266         const PetscInt padj = tmpAdj[q];
267         PetscInt ndof, ncdof, ngoff, nd;
268 
269         if ((padj < pStart) || (padj >= pEnd)) continue;
270         ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
271         ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
272         ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
273         for (nd = 0; nd < ndof-ncdof; ++nd) {
274           adj[aoff+i] = (ngoff < 0 ? -(ngoff+1) : ngoff) + nd;
275           ++i;
276         }
277       }
278     }
279   }
280   /* Debugging */
281   if (debug) {
282     IS tmp;
283     ierr = PetscPrintf(comm, "Leaf adjacency indices\n");CHKERRQ(ierr);
284     ierr = ISCreateGeneral(comm, adjSize, adj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
285     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
286     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
287   }
288   /* Gather adjacenct indices to root */
289   ierr = PetscSectionGetStorageSize(rootSectionAdj, &adjSize);CHKERRQ(ierr);
290   ierr = PetscMalloc1(adjSize, &rootAdj);CHKERRQ(ierr);
291   for (r = 0; r < adjSize; ++r) rootAdj[r] = -1;
292   if (doComm) {
293     ierr = PetscSFGatherBegin(sfAdj, MPIU_INT, adj, rootAdj);CHKERRQ(ierr);
294     ierr = PetscSFGatherEnd(sfAdj, MPIU_INT, adj, rootAdj);CHKERRQ(ierr);
295   }
296   ierr = PetscSFDestroy(&sfAdj);CHKERRQ(ierr);
297   ierr = PetscFree(adj);CHKERRQ(ierr);
298   /* Debugging */
299   if (debug) {
300     IS tmp;
301     ierr = PetscPrintf(comm, "Root adjacency indices after gather\n");CHKERRQ(ierr);
302     ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
303     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
304     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
305   }
306   /* Add in local adjacency indices for owned dofs on interface (roots) */
307   for (p = pStart; p < pEnd; ++p) {
308     PetscInt numAdj = maxAdjSize, adof, dof, off, d, q;
309 
310     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
311     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
312     if (!dof) continue;
313     ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
314     if (adof <= 0) continue;
315     ierr = DMPlexGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);CHKERRQ(ierr);
316     for (d = off; d < off+dof; ++d) {
317       PetscInt adof, aoff, i;
318 
319       ierr = PetscSectionGetDof(rootSectionAdj, d, &adof);CHKERRQ(ierr);
320       ierr = PetscSectionGetOffset(rootSectionAdj, d, &aoff);CHKERRQ(ierr);
321       i    = adof-1;
322       for (q = 0; q < numAdj; ++q) {
323         const PetscInt padj = tmpAdj[q];
324         PetscInt ndof, ncdof, ngoff, nd;
325 
326         if ((padj < pStart) || (padj >= pEnd)) continue;
327         ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
328         ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
329         ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
330         for (nd = 0; nd < ndof-ncdof; ++nd) {
331           rootAdj[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
332           --i;
333         }
334       }
335     }
336   }
337   /* Debugging */
338   if (debug) {
339     IS tmp;
340     ierr = PetscPrintf(comm, "Root adjacency indices\n");CHKERRQ(ierr);
341     ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
342     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
343     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
344   }
345   /* Compress indices */
346   ierr = PetscSectionSetUp(rootSectionAdj);CHKERRQ(ierr);
347   for (p = pStart; p < pEnd; ++p) {
348     PetscInt dof, cdof, off, d;
349     PetscInt adof, aoff;
350 
351     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
352     ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
353     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
354     if (!dof) continue;
355     ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
356     if (adof <= 0) continue;
357     for (d = off; d < off+dof-cdof; ++d) {
358       ierr = PetscSectionGetDof(rootSectionAdj, d, &adof);CHKERRQ(ierr);
359       ierr = PetscSectionGetOffset(rootSectionAdj, d, &aoff);CHKERRQ(ierr);
360       ierr = PetscSortRemoveDupsInt(&adof, &rootAdj[aoff]);CHKERRQ(ierr);
361       ierr = PetscSectionSetDof(rootSectionAdj, d, adof);CHKERRQ(ierr);
362     }
363   }
364   /* Debugging */
365   if (debug) {
366     IS tmp;
367     ierr = PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after compression:\n");CHKERRQ(ierr);
368     ierr = PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
369     ierr = PetscPrintf(comm, "Root adjacency indices after compression\n");CHKERRQ(ierr);
370     ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
371     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
372     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
373   }
374   /* Build adjacency section: Maps global indices to sets of adjacent global indices */
375   ierr = PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd);CHKERRQ(ierr);
376   ierr = PetscSectionCreate(comm, &sectionAdj);CHKERRQ(ierr);
377   ierr = PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd);CHKERRQ(ierr);
378   for (p = pStart; p < pEnd; ++p) {
379     PetscInt  numAdj = maxAdjSize, dof, cdof, off, goff, d, q;
380     PetscBool found  = PETSC_TRUE;
381 
382     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
383     ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
384     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
385     ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
386     for (d = 0; d < dof-cdof; ++d) {
387       PetscInt ldof, rdof;
388 
389       ierr = PetscSectionGetDof(leafSectionAdj, off+d, &ldof);CHKERRQ(ierr);
390       ierr = PetscSectionGetDof(rootSectionAdj, off+d, &rdof);CHKERRQ(ierr);
391       if (ldof > 0) {
392         /* We do not own this point */
393       } else if (rdof > 0) {
394         ierr = PetscSectionSetDof(sectionAdj, goff+d, rdof);CHKERRQ(ierr);
395       } else {
396         found = PETSC_FALSE;
397       }
398     }
399     if (found) continue;
400     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
401     ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
402     ierr = DMPlexGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);CHKERRQ(ierr);
403     for (q = 0; q < numAdj; ++q) {
404       const PetscInt padj = tmpAdj[q];
405       PetscInt ndof, ncdof, noff;
406 
407       if ((padj < pStart) || (padj >= pEnd)) continue;
408       ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
409       ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
410       ierr = PetscSectionGetOffset(section, padj, &noff);CHKERRQ(ierr);
411       for (d = goff; d < goff+dof-cdof; ++d) {
412         ierr = PetscSectionAddDof(sectionAdj, d, ndof-ncdof);CHKERRQ(ierr);
413       }
414     }
415   }
416   ierr = PetscSectionSetUp(sectionAdj);CHKERRQ(ierr);
417   if (debug) {
418     ierr = PetscPrintf(comm, "Adjacency Section for Preallocation:\n");CHKERRQ(ierr);
419     ierr = PetscSectionView(sectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
420   }
421   /* Get adjacent indices */
422   ierr = PetscSectionGetStorageSize(sectionAdj, &numCols);CHKERRQ(ierr);
423   ierr = PetscMalloc1(numCols, &cols);CHKERRQ(ierr);
424   for (p = pStart; p < pEnd; ++p) {
425     PetscInt  numAdj = maxAdjSize, dof, cdof, off, goff, d, q;
426     PetscBool found  = PETSC_TRUE;
427 
428     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
429     ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
430     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
431     ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
432     for (d = 0; d < dof-cdof; ++d) {
433       PetscInt ldof, rdof;
434 
435       ierr = PetscSectionGetDof(leafSectionAdj, off+d, &ldof);CHKERRQ(ierr);
436       ierr = PetscSectionGetDof(rootSectionAdj, off+d, &rdof);CHKERRQ(ierr);
437       if (ldof > 0) {
438         /* We do not own this point */
439       } else if (rdof > 0) {
440         PetscInt aoff, roff;
441 
442         ierr = PetscSectionGetOffset(sectionAdj, goff+d, &aoff);CHKERRQ(ierr);
443         ierr = PetscSectionGetOffset(rootSectionAdj, off+d, &roff);CHKERRQ(ierr);
444         ierr = PetscMemcpy(&cols[aoff], &rootAdj[roff], rdof * sizeof(PetscInt));CHKERRQ(ierr);
445       } else {
446         found = PETSC_FALSE;
447       }
448     }
449     if (found) continue;
450     ierr = DMPlexGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);CHKERRQ(ierr);
451     for (d = goff; d < goff+dof-cdof; ++d) {
452       PetscInt adof, aoff, i = 0;
453 
454       ierr = PetscSectionGetDof(sectionAdj, d, &adof);CHKERRQ(ierr);
455       ierr = PetscSectionGetOffset(sectionAdj, d, &aoff);CHKERRQ(ierr);
456       for (q = 0; q < numAdj; ++q) {
457         const PetscInt  padj = tmpAdj[q];
458         PetscInt        ndof, ncdof, ngoff, nd;
459         const PetscInt *ncind;
460 
461         /* Adjacent points may not be in the section chart */
462         if ((padj < pStart) || (padj >= pEnd)) continue;
463         ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
464         ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
465         ierr = PetscSectionGetConstraintIndices(section, padj, &ncind);CHKERRQ(ierr);
466         ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
467         for (nd = 0; nd < ndof-ncdof; ++nd, ++i) {
468           cols[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
469         }
470       }
471       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);
472     }
473   }
474   ierr = PetscSectionDestroy(&leafSectionAdj);CHKERRQ(ierr);
475   ierr = PetscSectionDestroy(&rootSectionAdj);CHKERRQ(ierr);
476   ierr = PetscFree(rootAdj);CHKERRQ(ierr);
477   ierr = PetscFree2(tmpClosure, tmpAdj);CHKERRQ(ierr);
478   /* Debugging */
479   if (debug) {
480     IS tmp;
481     ierr = PetscPrintf(comm, "Column indices\n");CHKERRQ(ierr);
482     ierr = ISCreateGeneral(comm, numCols, cols, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
483     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
484     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
485   }
486   /* Create allocation vectors from adjacency graph */
487   ierr = MatGetLocalSize(A, &locRows, NULL);CHKERRQ(ierr);
488   ierr = PetscLayoutCreate(PetscObjectComm((PetscObject)A), &rLayout);CHKERRQ(ierr);
489   ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr);
490   ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr);
491   ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr);
492   ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr);
493   ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr);
494   /* Only loop over blocks of rows */
495   if (rStart%bs || rEnd%bs) SETERRQ3(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid layout [%d, %d) for matrix, must be divisible by block size %d", rStart, rEnd, bs);
496   for (r = rStart/bs; r < rEnd/bs; ++r) {
497     const PetscInt row = r*bs;
498     PetscInt       numCols, cStart, c;
499 
500     ierr = PetscSectionGetDof(sectionAdj, row, &numCols);CHKERRQ(ierr);
501     ierr = PetscSectionGetOffset(sectionAdj, row, &cStart);CHKERRQ(ierr);
502     for (c = cStart; c < cStart+numCols; ++c) {
503       if ((cols[c] >= rStart*bs) && (cols[c] < rEnd*bs)) {
504         ++dnz[r-rStart];
505         if (cols[c] >= row) ++dnzu[r-rStart];
506       } else {
507         ++onz[r-rStart];
508         if (cols[c] >= row) ++onzu[r-rStart];
509       }
510     }
511   }
512   if (bs > 1) {
513     for (r = 0; r < locRows/bs; ++r) {
514       dnz[r]  /= bs;
515       onz[r]  /= bs;
516       dnzu[r] /= bs;
517       onzu[r] /= bs;
518     }
519   }
520   /* Set matrix pattern */
521   ierr = MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu);CHKERRQ(ierr);
522   ierr = MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
523   /* Check for symmetric storage */
524   ierr = MatGetType(A, &mtype);CHKERRQ(ierr);
525   ierr = PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);CHKERRQ(ierr);
526   ierr = PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);CHKERRQ(ierr);
527   ierr = PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);CHKERRQ(ierr);
528   if (isSymBlock || isSymSeqBlock || isSymMPIBlock) {ierr = MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);CHKERRQ(ierr);}
529   /* Fill matrix with zeros */
530   if (fillMatrix) {
531     PetscScalar *values;
532     PetscInt     maxRowLen = 0;
533 
534     for (r = rStart; r < rEnd; ++r) {
535       PetscInt len;
536 
537       ierr      = PetscSectionGetDof(sectionAdj, r, &len);CHKERRQ(ierr);
538       maxRowLen = PetscMax(maxRowLen, len);
539     }
540     ierr = PetscCalloc1(maxRowLen, &values);CHKERRQ(ierr);
541     for (r = rStart; r < rEnd; ++r) {
542       PetscInt numCols, cStart;
543 
544       ierr = PetscSectionGetDof(sectionAdj, r, &numCols);CHKERRQ(ierr);
545       ierr = PetscSectionGetOffset(sectionAdj, r, &cStart);CHKERRQ(ierr);
546       ierr = MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);CHKERRQ(ierr);
547     }
548     ierr = PetscFree(values);CHKERRQ(ierr);
549     ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
550     ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
551   }
552   ierr = PetscSectionDestroy(&sectionAdj);CHKERRQ(ierr);
553   ierr = PetscFree(cols);CHKERRQ(ierr);
554   ierr = PetscLogEventEnd(DMPLEX_Preallocate,dm,0,0,0);CHKERRQ(ierr);
555   PetscFunctionReturn(0);
556 }
557 
558 #if 0
559 #undef __FUNCT__
560 #define __FUNCT__ "DMPlexPreallocateOperator_2"
561 PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
562 {
563   PetscInt       *tmpClosure,*tmpAdj,*visits;
564   PetscInt        c,cStart,cEnd,pStart,pEnd;
565   PetscErrorCode  ierr;
566 
567   PetscFunctionBegin;
568   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
569   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
570   ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
571 
572   maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1));
573 
574   ierr    = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
575   npoints = pEnd - pStart;
576 
577   ierr = PetscMalloc3(maxClosureSize,&tmpClosure,npoints,&lvisits,npoints,&visits);CHKERRQ(ierr);
578   ierr = PetscMemzero(lvisits,(pEnd-pStart)*sizeof(PetscInt));CHKERRQ(ierr);
579   ierr = PetscMemzero(visits,(pEnd-pStart)*sizeof(PetscInt));CHKERRQ(ierr);
580   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
581   for (c=cStart; c<cEnd; c++) {
582     PetscInt *support = tmpClosure;
583     ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support);CHKERRQ(ierr);
584     for (p=0; p<supportSize; p++) lvisits[support[p]]++;
585   }
586   ierr = PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
587   ierr = PetscSFReduceEnd  (sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
588   ierr = PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits);CHKERRQ(ierr);
589   ierr = PetscSFBcastEnd  (sf,MPIU_INT,visits,lvisits);CHKERRQ(ierr);
590 
591   ierr = PetscSFGetRanks();CHKERRQ(ierr);
592 
593 
594   ierr = PetscMalloc2(maxClosureSize*maxClosureSize,&cellmat,npoints,&owner);CHKERRQ(ierr);
595   for (c=cStart; c<cEnd; c++) {
596     ierr = PetscMemzero(cellmat,maxClosureSize*maxClosureSize*sizeof(PetscInt));CHKERRQ(ierr);
597     /*
598      Depth-first walk of transitive closure.
599      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.
600      This contribution is added to dnz if owning ranks of p and q match, to onz otherwise.
601      */
602   }
603 
604   ierr = PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM);CHKERRQ(ierr);
605   ierr = PetscSFReduceEnd  (sf,MPIU_INT,lonz,onz,MPI_SUM);CHKERRQ(ierr);
606   PetscFunctionReturn(0);
607 }
608 #endif
609