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