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