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