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