xref: /petsc/src/dm/impls/plex/plexpreallocate.c (revision fb47e2fe627f166e041d762227e6ecf669c80696) 
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   }
237   /* Gather adjacenct indices to root */
238   ierr = PetscSectionGetStorageSize(rootSectionAdj, &adjSize);CHKERRQ(ierr);
239   ierr = PetscMalloc(adjSize * sizeof(PetscInt), &rootAdj);CHKERRQ(ierr);
240   for (r = 0; r < adjSize; ++r) rootAdj[r] = -1;
241   if (size > 1) {
242     ierr = PetscSFGatherBegin(sfAdj, MPIU_INT, adj, rootAdj);CHKERRQ(ierr);
243     ierr = PetscSFGatherEnd(sfAdj, MPIU_INT, adj, rootAdj);CHKERRQ(ierr);
244   }
245   ierr = PetscSFDestroy(&sfAdj);CHKERRQ(ierr);
246   ierr = PetscFree(adj);CHKERRQ(ierr);
247   /* Debugging */
248   if (debug) {
249     IS tmp;
250     ierr = PetscPrintf(comm, "Root adjacency indices after gather\n");CHKERRQ(ierr);
251     ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
252     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
253   }
254   /* Add in local adjacency indices for owned dofs on interface (roots) */
255   for (p = pStart; p < pEnd; ++p) {
256     PetscInt numAdj = maxAdjSize, adof, dof, off, d, q;
257 
258     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
259     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
260     if (!dof) continue;
261     ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
262     if (adof <= 0) continue;
263     ierr = DMPlexGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);CHKERRQ(ierr);
264     for (d = off; d < off+dof; ++d) {
265       PetscInt adof, aoff, i;
266 
267       ierr = PetscSectionGetDof(rootSectionAdj, d, &adof);CHKERRQ(ierr);
268       ierr = PetscSectionGetOffset(rootSectionAdj, d, &aoff);CHKERRQ(ierr);
269       i    = adof-1;
270       for (q = 0; q < numAdj; ++q) {
271         const PetscInt padj = tmpAdj[q];
272         PetscInt ndof, ncdof, ngoff, nd;
273 
274         if ((padj < pStart) || (padj >= pEnd)) continue;
275         ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
276         ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
277         ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
278         for (nd = 0; nd < ndof-ncdof; ++nd) {
279           rootAdj[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
280           --i;
281         }
282       }
283     }
284   }
285   /* Debugging */
286   if (debug) {
287     IS tmp;
288     ierr = PetscPrintf(comm, "Root adjacency indices\n");CHKERRQ(ierr);
289     ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
290     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
291   }
292   /* Compress indices */
293   ierr = PetscSectionSetUp(rootSectionAdj);CHKERRQ(ierr);
294   for (p = pStart; p < pEnd; ++p) {
295     PetscInt dof, cdof, off, d;
296     PetscInt adof, aoff;
297 
298     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
299     ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
300     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
301     if (!dof) continue;
302     ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
303     if (adof <= 0) continue;
304     for (d = off; d < off+dof-cdof; ++d) {
305       ierr = PetscSectionGetDof(rootSectionAdj, d, &adof);CHKERRQ(ierr);
306       ierr = PetscSectionGetOffset(rootSectionAdj, d, &aoff);CHKERRQ(ierr);
307       ierr = PetscSortRemoveDupsInt(&adof, &rootAdj[aoff]);CHKERRQ(ierr);
308       ierr = PetscSectionSetDof(rootSectionAdj, d, adof);CHKERRQ(ierr);
309     }
310   }
311   /* Debugging */
312   if (debug) {
313     IS tmp;
314     ierr = PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after compression:\n");CHKERRQ(ierr);
315     ierr = PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
316     ierr = PetscPrintf(comm, "Root adjacency indices after compression\n");CHKERRQ(ierr);
317     ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
318     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
319   }
320   /* Build adjacency section: Maps global indices to sets of adjacent global indices */
321   ierr = PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd);CHKERRQ(ierr);
322   ierr = PetscSectionCreate(comm, &sectionAdj);CHKERRQ(ierr);
323   ierr = PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd);CHKERRQ(ierr);
324   for (p = pStart; p < pEnd; ++p) {
325     PetscInt  numAdj = maxAdjSize, dof, cdof, off, goff, d, q;
326     PetscBool found  = PETSC_TRUE;
327 
328     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
329     ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
330     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
331     ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
332     for (d = 0; d < dof-cdof; ++d) {
333       PetscInt ldof, rdof;
334 
335       ierr = PetscSectionGetDof(leafSectionAdj, off+d, &ldof);CHKERRQ(ierr);
336       ierr = PetscSectionGetDof(rootSectionAdj, off+d, &rdof);CHKERRQ(ierr);
337       if (ldof > 0) {
338         /* We do not own this point */
339       } else if (rdof > 0) {
340         ierr = PetscSectionSetDof(sectionAdj, goff+d, rdof);CHKERRQ(ierr);
341       } else {
342         found = PETSC_FALSE;
343       }
344     }
345     if (found) continue;
346     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
347     ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
348     ierr = DMPlexGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);CHKERRQ(ierr);
349     for (q = 0; q < numAdj; ++q) {
350       const PetscInt padj = tmpAdj[q];
351       PetscInt ndof, ncdof, noff;
352 
353       if ((padj < pStart) || (padj >= pEnd)) continue;
354       ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
355       ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
356       ierr = PetscSectionGetOffset(section, padj, &noff);CHKERRQ(ierr);
357       for (d = goff; d < goff+dof-cdof; ++d) {
358         ierr = PetscSectionAddDof(sectionAdj, d, ndof-ncdof);CHKERRQ(ierr);
359       }
360     }
361   }
362   ierr = PetscSectionSetUp(sectionAdj);CHKERRQ(ierr);
363   if (debug) {
364     ierr = PetscPrintf(comm, "Adjacency Section for Preallocation:\n");CHKERRQ(ierr);
365     ierr = PetscSectionView(sectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
366   }
367   /* Get adjacent indices */
368   ierr = PetscSectionGetStorageSize(sectionAdj, &numCols);CHKERRQ(ierr);
369   ierr = PetscMalloc(numCols * sizeof(PetscInt), &cols);CHKERRQ(ierr);
370   for (p = pStart; p < pEnd; ++p) {
371     PetscInt  numAdj = maxAdjSize, dof, cdof, off, goff, d, q;
372     PetscBool found  = PETSC_TRUE;
373 
374     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
375     ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
376     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
377     ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
378     for (d = 0; d < dof-cdof; ++d) {
379       PetscInt ldof, rdof;
380 
381       ierr = PetscSectionGetDof(leafSectionAdj, off+d, &ldof);CHKERRQ(ierr);
382       ierr = PetscSectionGetDof(rootSectionAdj, off+d, &rdof);CHKERRQ(ierr);
383       if (ldof > 0) {
384         /* We do not own this point */
385       } else if (rdof > 0) {
386         PetscInt aoff, roff;
387 
388         ierr = PetscSectionGetOffset(sectionAdj, goff+d, &aoff);CHKERRQ(ierr);
389         ierr = PetscSectionGetOffset(rootSectionAdj, off+d, &roff);CHKERRQ(ierr);
390         ierr = PetscMemcpy(&cols[aoff], &rootAdj[roff], rdof * sizeof(PetscInt));CHKERRQ(ierr);
391       } else {
392         found = PETSC_FALSE;
393       }
394     }
395     if (found) continue;
396     ierr = DMPlexGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);CHKERRQ(ierr);
397     for (d = goff; d < goff+dof-cdof; ++d) {
398       PetscInt adof, aoff, i = 0;
399 
400       ierr = PetscSectionGetDof(sectionAdj, d, &adof);CHKERRQ(ierr);
401       ierr = PetscSectionGetOffset(sectionAdj, d, &aoff);CHKERRQ(ierr);
402       for (q = 0; q < numAdj; ++q) {
403         const PetscInt  padj = tmpAdj[q];
404         PetscInt        ndof, ncdof, ngoff, nd;
405         const PetscInt *ncind;
406 
407         /* Adjacent points may not be in the section chart */
408         if ((padj < pStart) || (padj >= pEnd)) continue;
409         ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
410         ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
411         ierr = PetscSectionGetConstraintIndices(section, padj, &ncind);CHKERRQ(ierr);
412         ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
413         for (nd = 0; nd < ndof-ncdof; ++nd, ++i) {
414           cols[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
415         }
416       }
417       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);
418     }
419   }
420   ierr = PetscSectionDestroy(&leafSectionAdj);CHKERRQ(ierr);
421   ierr = PetscSectionDestroy(&rootSectionAdj);CHKERRQ(ierr);
422   ierr = PetscFree(rootAdj);CHKERRQ(ierr);
423   ierr = PetscFree2(tmpClosure, tmpAdj);CHKERRQ(ierr);
424   /* Debugging */
425   if (debug) {
426     IS tmp;
427     ierr = PetscPrintf(comm, "Column indices\n");CHKERRQ(ierr);
428     ierr = ISCreateGeneral(comm, numCols, cols, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
429     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
430   }
431   /* Create allocation vectors from adjacency graph */
432   ierr = MatGetLocalSize(A, &locRows, NULL);CHKERRQ(ierr);
433   ierr = PetscLayoutCreate(PetscObjectComm((PetscObject)A), &rLayout);CHKERRQ(ierr);
434   ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr);
435   ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr);
436   ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr);
437   ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr);
438   ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr);
439   /* Only loop over blocks of rows */
440   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);
441   for (r = rStart/bs; r < rEnd/bs; ++r) {
442     const PetscInt row = r*bs;
443     PetscInt       numCols, cStart, c;
444 
445     ierr = PetscSectionGetDof(sectionAdj, row, &numCols);CHKERRQ(ierr);
446     ierr = PetscSectionGetOffset(sectionAdj, row, &cStart);CHKERRQ(ierr);
447     for (c = cStart; c < cStart+numCols; ++c) {
448       if ((cols[c] >= rStart*bs) && (cols[c] < rEnd*bs)) {
449         ++dnz[r-rStart];
450         if (cols[c] >= row) ++dnzu[r-rStart];
451       } else {
452         ++onz[r-rStart];
453         if (cols[c] >= row) ++onzu[r-rStart];
454       }
455     }
456   }
457   if (bs > 1) {
458     for (r = 0; r < locRows/bs; ++r) {
459       dnz[r]  /= bs;
460       onz[r]  /= bs;
461       dnzu[r] /= bs;
462       onzu[r] /= bs;
463     }
464   }
465   /* Set matrix pattern */
466   ierr = MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu);CHKERRQ(ierr);
467   ierr = MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
468   /* Check for symmetric storage */
469   ierr = MatGetType(A, &mtype);CHKERRQ(ierr);
470   ierr = PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);CHKERRQ(ierr);
471   ierr = PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);CHKERRQ(ierr);
472   ierr = PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);CHKERRQ(ierr);
473   if (isSymBlock || isSymSeqBlock || isSymMPIBlock) {ierr = MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);CHKERRQ(ierr);}
474   /* Fill matrix with zeros */
475   if (fillMatrix) {
476     PetscScalar *values;
477     PetscInt     maxRowLen = 0;
478 
479     for (r = rStart; r < rEnd; ++r) {
480       PetscInt len;
481 
482       ierr      = PetscSectionGetDof(sectionAdj, r, &len);CHKERRQ(ierr);
483       maxRowLen = PetscMax(maxRowLen, len);
484     }
485     ierr = PetscMalloc(maxRowLen * sizeof(PetscScalar), &values);CHKERRQ(ierr);
486     ierr = PetscMemzero(values, maxRowLen * sizeof(PetscScalar));CHKERRQ(ierr);
487     for (r = rStart; r < rEnd; ++r) {
488       PetscInt numCols, cStart;
489 
490       ierr = PetscSectionGetDof(sectionAdj, r, &numCols);CHKERRQ(ierr);
491       ierr = PetscSectionGetOffset(sectionAdj, r, &cStart);CHKERRQ(ierr);
492       ierr = MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);CHKERRQ(ierr);
493     }
494     ierr = PetscFree(values);CHKERRQ(ierr);
495     ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
496     ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
497   }
498   ierr = PetscSectionDestroy(&sectionAdj);CHKERRQ(ierr);
499   ierr = PetscFree(cols);CHKERRQ(ierr);
500   PetscFunctionReturn(0);
501 }
502 
503 #if 0
504 #undef __FUNCT__
505 #define __FUNCT__ "DMPlexPreallocateOperator_2"
506 PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
507 {
508   PetscInt       *tmpClosure,*tmpAdj,*visits;
509   PetscInt        c,cStart,cEnd,pStart,pEnd;
510   PetscErrorCode  ierr;
511 
512   PetscFunctionBegin;
513   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
514   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
515   ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
516 
517   maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1));
518 
519   ierr    = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
520   npoints = pEnd - pStart;
521 
522   ierr = PetscMalloc3(maxClosureSize,PetscInt,&tmpClosure,npoints,PetscInt,&lvisits,npoints,PetscInt,&visits);CHKERRQ(ierr);
523   ierr = PetscMemzero(lvisits,(pEnd-pStart)*sizeof(PetscInt));CHKERRQ(ierr);
524   ierr = PetscMemzero(visits,(pEnd-pStart)*sizeof(PetscInt));CHKERRQ(ierr);
525   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
526   for (c=cStart; c<cEnd; c++) {
527     PetscInt *support = tmpClosure;
528     ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support);CHKERRQ(ierr);
529     for (p=0; p<supportSize; p++) lvisits[support[p]]++;
530   }
531   ierr = PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
532   ierr = PetscSFReduceEnd  (sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
533   ierr = PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits);CHKERRQ(ierr);
534   ierr = PetscSFBcastEnd  (sf,MPIU_INT,visits,lvisits);CHKERRQ(ierr);
535 
536   ierr = PetscSFGetRanks();CHKERRQ(ierr);
537 
538 
539   ierr = PetscMalloc2(maxClosureSize*maxClosureSize,PetscInt,&cellmat,npoints,PetscInt,&owner);CHKERRQ(ierr);
540   for (c=cStart; c<cEnd; c++) {
541     ierr = PetscMemzero(cellmat,maxClosureSize*maxClosureSize*sizeof(PetscInt));CHKERRQ(ierr);
542     /*
543      Depth-first walk of transitive closure.
544      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.
545      This contribution is added to dnz if owning ranks of p and q match, to onz otherwise.
546      */
547   }
548 
549   ierr = PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM);CHKERRQ(ierr);
550   ierr = PetscSFReduceEnd  (sf,MPIU_INT,lonz,onz,MPI_SUM);CHKERRQ(ierr);
551   PetscFunctionReturn(0);
552 }
553 #endif
554