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