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