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