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