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