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