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