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