xref: /petsc/src/dm/impls/plex/plexpreallocate.c (revision fe998a80077c9ee0917a39496df43fc256e1b478) 
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__ "DMPlexComputeAnchorAdjacencies"
7 /* get adjacencies due to point-to-point constraints that can't be found with DMPlexGetAdjacency() */
8 static PetscErrorCode DMPlexComputeAnchorAdjacencies(DM dm, PetscSection section, PetscSection sectionGlobal, PetscSection *anchorSectionAdj, PetscInt *anchorAdj[])
9 {
10   PetscInt       pStart, pEnd;
11   PetscSection   adjSec, aSec;
12   IS             aIS;
13   PetscErrorCode ierr;
14 
15   PetscFunctionBegin;
16 
17   ierr = PetscSectionCreate(PetscObjectComm((PetscObject)section),&adjSec);CHKERRQ(ierr);
18   ierr = PetscSectionGetChart(section,&pStart,&pEnd);CHKERRQ(ierr);
19   ierr = PetscSectionSetChart(adjSec,pStart,pEnd);CHKERRQ(ierr);
20 
21   ierr = DMPlexGetAnchors(dm,&aSec,&aIS);CHKERRQ(ierr);
22   if (aSec) {
23     const PetscInt *anchors;
24     PetscInt       p, q, a, aSize, *offsets, aStart, aEnd, *inverse, iSize, *adj, adjSize;
25     PetscInt       *tmpAdjP = NULL, *tmpAdjQ = NULL;
26     PetscSection   inverseSec;
27 
28     /* invert the constraint-to-anchor map */
29     ierr = PetscSectionCreate(PetscObjectComm((PetscObject)aSec),&inverseSec);CHKERRQ(ierr);
30     ierr = PetscSectionSetChart(inverseSec,pStart,pEnd);CHKERRQ(ierr);
31     ierr = ISGetLocalSize(aIS, &aSize);CHKERRQ(ierr);
32     ierr = ISGetIndices(aIS, &anchors);CHKERRQ(ierr);
33 
34     for (p = 0; p < aSize; p++) {
35       PetscInt a = anchors[p];
36 
37       ierr = PetscSectionAddDof(inverseSec,a,1);CHKERRQ(ierr);
38     }
39     ierr = PetscSectionSetUp(inverseSec);CHKERRQ(ierr);
40     ierr = PetscSectionGetStorageSize(inverseSec,&iSize);CHKERRQ(ierr);
41     ierr = PetscMalloc1(iSize,&inverse);CHKERRQ(ierr);
42     ierr = PetscCalloc1(pEnd-pStart,&offsets);CHKERRQ(ierr);
43     ierr = PetscSectionGetChart(aSec,&aStart,&aEnd);CHKERRQ(ierr);
44     for (p = aStart; p < aEnd; p++) {
45       PetscInt dof, off;
46 
47       ierr = PetscSectionGetDof(aSec, p, &dof);CHKERRQ(ierr);
48       ierr = PetscSectionGetOffset(aSec, p, &off);CHKERRQ(ierr);
49 
50       for (q = 0; q < dof; q++) {
51         PetscInt iOff;
52 
53         a = anchors[off + q];
54         ierr = PetscSectionGetOffset(inverseSec, a, &iOff);CHKERRQ(ierr);
55         inverse[iOff + offsets[a-pStart]++] = p;
56       }
57     }
58     ierr = ISRestoreIndices(aIS, &anchors);CHKERRQ(ierr);
59     ierr = PetscFree(offsets);CHKERRQ(ierr);
60 
61     /* construct anchorAdj and adjSec
62      *
63      * loop over anchors:
64      *   construct anchor adjacency
65      *   loop over constrained:
66      *     construct constrained adjacency
67      *     if not in anchor adjacency, add to dofs
68      * setup adjSec, allocate anchorAdj
69      * loop over anchors:
70      *   construct anchor adjacency
71      *   loop over constrained:
72      *     construct constrained adjacency
73      *     if not in anchor adjacency
74      *       if not already in list, put in list
75      *   sort, unique, reduce dof count
76      * optional: compactify
77      */
78     for (p = pStart; p < pEnd; p++) {
79       PetscInt iDof, iOff, i, r, s, numAdjP = PETSC_DETERMINE;
80 
81       ierr = PetscSectionGetDof(inverseSec,p,&iDof);CHKERRQ(ierr);
82       if (!iDof) continue;
83       ierr = PetscSectionGetOffset(inverseSec,p,&iOff);CHKERRQ(ierr);
84       ierr = DMPlexGetAdjacency(dm,p,&numAdjP,&tmpAdjP);CHKERRQ(ierr);
85       for (i = 0; i < iDof; i++) {
86         PetscInt iNew = 0, qAdj, qAdjDof, qAdjCDof, numAdjQ = PETSC_DETERMINE;
87 
88         q = inverse[iOff + i];
89         ierr = DMPlexGetAdjacency(dm,q,&numAdjQ,&tmpAdjQ);CHKERRQ(ierr);
90         for (r = 0; r < numAdjQ; r++) {
91           qAdj = tmpAdjQ[r];
92           if ((qAdj < pStart) || (qAdj >= pEnd)) continue;
93           for (s = 0; s < numAdjP; s++) {
94             if (qAdj == tmpAdjP[s]) break;
95           }
96           if (s < numAdjP) continue;
97           ierr  = PetscSectionGetDof(section,qAdj,&qAdjDof);CHKERRQ(ierr);
98           ierr  = PetscSectionGetConstraintDof(section,qAdj,&qAdjCDof);CHKERRQ(ierr);
99           iNew += qAdjDof - qAdjCDof;
100         }
101         ierr = PetscSectionAddDof(adjSec,p,iNew);CHKERRQ(ierr);
102       }
103     }
104 
105     ierr = PetscSectionSetUp(adjSec);CHKERRQ(ierr);
106     ierr = PetscSectionGetStorageSize(adjSec,&adjSize);CHKERRQ(ierr);
107     ierr = PetscMalloc1(adjSize,&adj);CHKERRQ(ierr);
108 
109     for (p = pStart; p < pEnd; p++) {
110       PetscInt iDof, iOff, i, r, s, aOff, aOffOrig, aDof, numAdjP = PETSC_DETERMINE;
111 
112       ierr = PetscSectionGetDof(inverseSec,p,&iDof);CHKERRQ(ierr);
113       if (!iDof) continue;
114       ierr = PetscSectionGetOffset(inverseSec,p,&iOff);CHKERRQ(ierr);
115       ierr = DMPlexGetAdjacency(dm,p,&numAdjP,&tmpAdjP);CHKERRQ(ierr);
116       ierr = PetscSectionGetDof(adjSec,p,&aDof);CHKERRQ(ierr);
117       ierr = PetscSectionGetOffset(adjSec,p,&aOff);CHKERRQ(ierr);
118       aOffOrig = aOff;
119       for (i = 0; i < iDof; i++) {
120         PetscInt qAdj, qAdjDof, qAdjCDof, qAdjOff, nd, numAdjQ = PETSC_DETERMINE;
121 
122         q = inverse[iOff + i];
123         ierr = DMPlexGetAdjacency(dm,q,&numAdjQ,&tmpAdjQ);CHKERRQ(ierr);
124         for (r = 0; r < numAdjQ; r++) {
125           qAdj = tmpAdjQ[r];
126           if ((qAdj < pStart) || (qAdj >= pEnd)) continue;
127           for (s = 0; s < numAdjP; s++) {
128             if (qAdj == tmpAdjP[s]) break;
129           }
130           if (s < numAdjP) continue;
131           ierr  = PetscSectionGetDof(section,qAdj,&qAdjDof);CHKERRQ(ierr);
132           ierr  = PetscSectionGetConstraintDof(section,qAdj,&qAdjCDof);CHKERRQ(ierr);
133           ierr  = PetscSectionGetOffset(sectionGlobal,qAdj,&qAdjOff);CHKERRQ(ierr);
134           for (nd = 0; nd < qAdjDof-qAdjCDof; ++nd) {
135             adj[aOff++] = (qAdjOff < 0 ? -(qAdjOff+1) : qAdjOff) + nd;
136           }
137         }
138       }
139       ierr = PetscSortRemoveDupsInt(&aDof,&adj[aOffOrig]);CHKERRQ(ierr);
140       ierr = PetscSectionSetDof(adjSec,p,aDof);CHKERRQ(ierr);
141     }
142     *anchorAdj = adj;
143 
144     /* clean up */
145     ierr = PetscSectionDestroy(&inverseSec);CHKERRQ(ierr);
146     ierr = PetscFree(inverse);CHKERRQ(ierr);
147     ierr = PetscFree(tmpAdjP);CHKERRQ(ierr);
148     ierr = PetscFree(tmpAdjQ);CHKERRQ(ierr);
149   }
150   else {
151     *anchorAdj = NULL;
152     ierr = PetscSectionSetUp(adjSec);CHKERRQ(ierr);
153   }
154   *anchorSectionAdj = adjSec;
155   PetscFunctionReturn(0);
156 }
157 
158 #undef __FUNCT__
159 #define __FUNCT__ "DMPlexPreallocateOperator"
160 PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
161 {
162   MPI_Comm           comm;
163   MatType            mtype;
164   PetscSF            sf, sfDof, sfAdj;
165   PetscSection       leafSectionAdj, rootSectionAdj, sectionAdj, anchorSectionAdj;
166   PetscInt           nroots, nleaves, l, p;
167   const PetscInt    *leaves;
168   const PetscSFNode *remotes;
169   PetscInt           dim, pStart, pEnd, numDof, globalOffStart, globalOffEnd, numCols;
170   PetscInt          *tmpAdj = NULL, *adj, *rootAdj, *anchorAdj = NULL, *cols, *remoteOffsets;
171   PetscInt           adjSize;
172   PetscLayout        rLayout;
173   PetscInt           locRows, rStart, rEnd, r;
174   PetscMPIInt        size;
175   PetscBool          doCommLocal, doComm, debug = PETSC_FALSE, isSymBlock, isSymSeqBlock, isSymMPIBlock;
176   PetscBool          useAnchors;
177   PetscErrorCode     ierr;
178 
179   PetscFunctionBegin;
180   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
181   PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 3);
182   PetscValidHeaderSpecific(sectionGlobal, PETSC_SECTION_CLASSID, 4);
183   PetscValidHeaderSpecific(A, MAT_CLASSID, 9);
184   if (dnz)  PetscValidPointer(dnz,5);
185   if (onz)  PetscValidPointer(onz,6);
186   if (dnzu) PetscValidPointer(dnzu,7);
187   if (onzu) PetscValidPointer(onzu,8);
188   ierr = PetscLogEventBegin(DMPLEX_Preallocate,dm,0,0,0);CHKERRQ(ierr);
189   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
190   ierr = PetscOptionsGetBool(NULL, "-dm_view_preallocation", &debug, NULL);CHKERRQ(ierr);
191   ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
192   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
193   ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
194   ierr = PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);CHKERRQ(ierr);
195   doCommLocal = (size > 1) && (nroots >= 0) ? PETSC_TRUE : PETSC_FALSE;
196   ierr = MPI_Allreduce(&doCommLocal, &doComm, 1, MPIU_BOOL, MPI_LAND, comm);CHKERRQ(ierr);
197   /* Create dof SF based on point SF */
198   if (debug) {
199     ierr = PetscPrintf(comm, "Input Section for Preallocation:\n");CHKERRQ(ierr);
200     ierr = PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
201     ierr = PetscPrintf(comm, "Input Global Section for Preallocation:\n");CHKERRQ(ierr);
202     ierr = PetscSectionView(sectionGlobal, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
203     ierr = PetscPrintf(comm, "Input SF for Preallocation:\n");CHKERRQ(ierr);
204     ierr = PetscSFView(sf, NULL);CHKERRQ(ierr);
205   }
206   ierr = PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets);CHKERRQ(ierr);
207   ierr = PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof);CHKERRQ(ierr);
208   if (debug) {
209     ierr = PetscPrintf(comm, "Dof SF for Preallocation:\n");CHKERRQ(ierr);
210     ierr = PetscSFView(sfDof, NULL);CHKERRQ(ierr);
211   }
212   /* Create section for dof adjacency (dof ==> # adj dof) */
213   ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
214   ierr = PetscSectionGetStorageSize(section, &numDof);CHKERRQ(ierr);
215   ierr = PetscSectionCreate(comm, &leafSectionAdj);CHKERRQ(ierr);
216   ierr = PetscSectionSetChart(leafSectionAdj, 0, numDof);CHKERRQ(ierr);
217   ierr = PetscSectionCreate(comm, &rootSectionAdj);CHKERRQ(ierr);
218   ierr = PetscSectionSetChart(rootSectionAdj, 0, numDof);CHKERRQ(ierr);
219   /*   Fill in the ghost dofs on the interface */
220   ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, &remotes);CHKERRQ(ierr);
221   /* use constraints in finding adjacency in this routine */
222   ierr = DMPlexGetAdjacencyUseAnchors(dm,&useAnchors);CHKERRQ(ierr);
223   ierr = DMPlexSetAdjacencyUseAnchors(dm,PETSC_TRUE);CHKERRQ(ierr);
224 
225   /*
226    section        - maps points to (# dofs, local dofs)
227    sectionGlobal  - maps points to (# dofs, global dofs)
228    leafSectionAdj - maps unowned local dofs to # adj dofs
229    rootSectionAdj - maps   owned local dofs to # adj dofs
230    adj            - adj global dofs indexed by leafSectionAdj
231    rootAdj        - adj global dofs indexed by rootSectionAdj
232    sf    - describes shared points across procs
233    sfDof - describes shared dofs across procs
234    sfAdj - describes shared adjacent dofs across procs
235    ** The bootstrapping process involves six rounds with similar structure of visiting neighbors of each point.
236   (0). If there are point-to-point constraints, add the adjacencies of constrained points to anchors in anchorAdj
237        (This is done in DMPlexComputeAnchorAdjacencies())
238     1. Visit unowned points on interface, count adjacencies placing in leafSectionAdj
239        Reduce those counts to rootSectionAdj (now redundantly counting some interface points)
240     2. Visit owned points on interface, count adjacencies placing in rootSectionAdj
241        Create sfAdj connecting rootSectionAdj and leafSectionAdj
242     3. Visit unowned points on interface, write adjacencies to adj
243        Gather adj to rootAdj (note that there is redundancy in rootAdj when multiple procs find the same adjacencies)
244     4. Visit owned points on interface, write adjacencies to rootAdj
245        Remove redundancy in rootAdj
246    ** The last two traversals use transitive closure
247     5. Visit all owned points in the subdomain, count dofs for each point (sectionAdj)
248        Allocate memory addressed by sectionAdj (cols)
249     6. Visit all owned points in the subdomain, insert dof adjacencies into cols
250    ** Knowing all the column adjacencies, check ownership and sum into dnz and onz
251   */
252 
253   ierr = DMPlexComputeAnchorAdjacencies(dm,section,sectionGlobal,&anchorSectionAdj,&anchorAdj);CHKERRQ(ierr);
254 
255   for (l = 0; l < nleaves; ++l) {
256     PetscInt dof, off, d, q, anDof;
257     PetscInt p = leaves[l], numAdj = PETSC_DETERMINE;
258 
259     if ((p < pStart) || (p >= pEnd)) continue;
260     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
261     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
262     ierr = DMPlexGetAdjacency(dm, p, &numAdj, &tmpAdj);CHKERRQ(ierr);
263     for (q = 0; q < numAdj; ++q) {
264       const PetscInt padj = tmpAdj[q];
265       PetscInt ndof, ncdof;
266 
267       if ((padj < pStart) || (padj >= pEnd)) continue;
268       ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
269       ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
270       for (d = off; d < off+dof; ++d) {
271         ierr = PetscSectionAddDof(leafSectionAdj, d, ndof-ncdof);CHKERRQ(ierr);
272       }
273     }
274     ierr = PetscSectionGetDof(anchorSectionAdj, p, &anDof);CHKERRQ(ierr);
275     if (anDof) {
276       for (d = off; d < off+dof; ++d) {
277         ierr = PetscSectionAddDof(leafSectionAdj, d, anDof);CHKERRQ(ierr);
278       }
279     }
280   }
281   ierr = PetscSectionSetUp(leafSectionAdj);CHKERRQ(ierr);
282   if (debug) {
283     ierr = PetscPrintf(comm, "Adjacency Section for Preallocation on Leaves:\n");CHKERRQ(ierr);
284     ierr = PetscSectionView(leafSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
285   }
286   /* Get maximum remote adjacency sizes for owned dofs on interface (roots) */
287   if (doComm) {
288     ierr = PetscSFReduceBegin(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);CHKERRQ(ierr);
289     ierr = PetscSFReduceEnd(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);CHKERRQ(ierr);
290   }
291   if (debug) {
292     ierr = PetscPrintf(comm, "Adjancency Section for Preallocation on Roots:\n");CHKERRQ(ierr);
293     ierr = PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
294   }
295   /* Add in local adjacency sizes for owned dofs on interface (roots) */
296   for (p = pStart; p < pEnd; ++p) {
297     PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof;
298 
299     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
300     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
301     if (!dof) continue;
302     ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
303     if (adof <= 0) continue;
304     ierr = DMPlexGetAdjacency(dm, p, &numAdj, &tmpAdj);CHKERRQ(ierr);
305     for (q = 0; q < numAdj; ++q) {
306       const PetscInt padj = tmpAdj[q];
307       PetscInt ndof, ncdof;
308 
309       if ((padj < pStart) || (padj >= pEnd)) continue;
310       ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
311       ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
312       for (d = off; d < off+dof; ++d) {
313         ierr = PetscSectionAddDof(rootSectionAdj, d, ndof-ncdof);CHKERRQ(ierr);
314       }
315     }
316     ierr = PetscSectionGetDof(anchorSectionAdj, p, &anDof);CHKERRQ(ierr);
317     if (anDof) {
318       for (d = off; d < off+dof; ++d) {
319         ierr = PetscSectionAddDof(rootSectionAdj, d, anDof);CHKERRQ(ierr);
320       }
321     }
322   }
323   ierr = PetscSectionSetUp(rootSectionAdj);CHKERRQ(ierr);
324   if (debug) {
325     ierr = PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after local additions:\n");CHKERRQ(ierr);
326     ierr = PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
327   }
328   /* Create adj SF based on dof SF */
329   ierr = PetscSFCreateRemoteOffsets(sfDof, rootSectionAdj, leafSectionAdj, &remoteOffsets);CHKERRQ(ierr);
330   ierr = PetscSFCreateSectionSF(sfDof, rootSectionAdj, remoteOffsets, leafSectionAdj, &sfAdj);CHKERRQ(ierr);
331   if (debug) {
332     ierr = PetscPrintf(comm, "Adjacency SF for Preallocation:\n");CHKERRQ(ierr);
333     ierr = PetscSFView(sfAdj, NULL);CHKERRQ(ierr);
334   }
335   ierr = PetscSFDestroy(&sfDof);CHKERRQ(ierr);
336   /* Create leaf adjacency */
337   ierr = PetscSectionSetUp(leafSectionAdj);CHKERRQ(ierr);
338   ierr = PetscSectionGetStorageSize(leafSectionAdj, &adjSize);CHKERRQ(ierr);
339   ierr = PetscCalloc1(adjSize, &adj);CHKERRQ(ierr);
340   for (l = 0; l < nleaves; ++l) {
341     PetscInt dof, off, d, q, anDof, anOff;
342     PetscInt p = leaves[l], numAdj = PETSC_DETERMINE;
343 
344     if ((p < pStart) || (p >= pEnd)) continue;
345     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
346     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
347     ierr = DMPlexGetAdjacency(dm, p, &numAdj, &tmpAdj);CHKERRQ(ierr);
348     ierr = PetscSectionGetDof(anchorSectionAdj, p, &anDof);CHKERRQ(ierr);
349     ierr = PetscSectionGetOffset(anchorSectionAdj, p, &anOff);CHKERRQ(ierr);
350     for (d = off; d < off+dof; ++d) {
351       PetscInt aoff, i = 0;
352 
353       ierr = PetscSectionGetOffset(leafSectionAdj, d, &aoff);CHKERRQ(ierr);
354       for (q = 0; q < numAdj; ++q) {
355         const PetscInt padj = tmpAdj[q];
356         PetscInt ndof, ncdof, ngoff, nd;
357 
358         if ((padj < pStart) || (padj >= pEnd)) continue;
359         ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
360         ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
361         ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
362         for (nd = 0; nd < ndof-ncdof; ++nd) {
363           adj[aoff+i] = (ngoff < 0 ? -(ngoff+1) : ngoff) + nd;
364           ++i;
365         }
366       }
367       for (q = 0; q < anDof; q++) {
368         adj[aoff+i] = anchorAdj[anOff+q];
369         ++i;
370       }
371     }
372   }
373   /* Debugging */
374   if (debug) {
375     IS tmp;
376     ierr = PetscPrintf(comm, "Leaf adjacency indices\n");CHKERRQ(ierr);
377     ierr = ISCreateGeneral(comm, adjSize, adj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
378     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
379     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
380   }
381   /* Gather adjacent indices to root */
382   ierr = PetscSectionGetStorageSize(rootSectionAdj, &adjSize);CHKERRQ(ierr);
383   ierr = PetscMalloc1(adjSize, &rootAdj);CHKERRQ(ierr);
384   for (r = 0; r < adjSize; ++r) rootAdj[r] = -1;
385   if (doComm) {
386     const PetscInt *indegree;
387     PetscInt       *remoteadj, radjsize = 0;
388 
389     ierr = PetscSFComputeDegreeBegin(sfAdj, &indegree);CHKERRQ(ierr);
390     ierr = PetscSFComputeDegreeEnd(sfAdj, &indegree);CHKERRQ(ierr);
391     for (p = 0; p < adjSize; ++p) radjsize += indegree[p];
392     ierr = PetscMalloc1(radjsize, &remoteadj);CHKERRQ(ierr);
393     ierr = PetscSFGatherBegin(sfAdj, MPIU_INT, adj, remoteadj);CHKERRQ(ierr);
394     ierr = PetscSFGatherEnd(sfAdj, MPIU_INT, adj, remoteadj);CHKERRQ(ierr);
395     for (p = 0, l = 0, r = 0; p < adjSize; ++p, l = PetscMax(p, l + indegree[p-1])) {
396       PetscInt s;
397       for (s = 0; s < indegree[p]; ++s, ++r) rootAdj[l+s] = remoteadj[r];
398     }
399     if (r != radjsize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %d != %d", r, radjsize);
400     if (l != adjSize)  SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %d != %d", l, adjSize);
401     ierr = PetscFree(remoteadj);CHKERRQ(ierr);
402   }
403   ierr = PetscSFDestroy(&sfAdj);CHKERRQ(ierr);
404   ierr = PetscFree(adj);CHKERRQ(ierr);
405   /* Debugging */
406   if (debug) {
407     IS tmp;
408     ierr = PetscPrintf(comm, "Root adjacency indices after gather\n");CHKERRQ(ierr);
409     ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
410     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
411     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
412   }
413   /* Add in local adjacency indices for owned dofs on interface (roots) */
414   for (p = pStart; p < pEnd; ++p) {
415     PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof, anOff;
416 
417     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
418     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
419     if (!dof) continue;
420     ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
421     if (adof <= 0) continue;
422     ierr = DMPlexGetAdjacency(dm, p, &numAdj, &tmpAdj);CHKERRQ(ierr);
423     ierr = PetscSectionGetDof(anchorSectionAdj, p, &anDof);CHKERRQ(ierr);
424     ierr = PetscSectionGetOffset(anchorSectionAdj, p, &anOff);CHKERRQ(ierr);
425     for (d = off; d < off+dof; ++d) {
426       PetscInt adof, aoff, i;
427 
428       ierr = PetscSectionGetDof(rootSectionAdj, d, &adof);CHKERRQ(ierr);
429       ierr = PetscSectionGetOffset(rootSectionAdj, d, &aoff);CHKERRQ(ierr);
430       i    = adof-1;
431       for (q = 0; q < anDof; q++) {
432         rootAdj[aoff+i] = anchorAdj[anOff+q];
433         --i;
434       }
435       for (q = 0; q < numAdj; ++q) {
436         const PetscInt padj = tmpAdj[q];
437         PetscInt ndof, ncdof, ngoff, nd;
438 
439         if ((padj < pStart) || (padj >= pEnd)) continue;
440         ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
441         ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
442         ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
443         for (nd = 0; nd < ndof-ncdof; ++nd) {
444           rootAdj[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
445           --i;
446         }
447       }
448     }
449   }
450   /* Debugging */
451   if (debug) {
452     IS tmp;
453     ierr = PetscPrintf(comm, "Root adjacency indices\n");CHKERRQ(ierr);
454     ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
455     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
456     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
457   }
458   /* Compress indices */
459   ierr = PetscSectionSetUp(rootSectionAdj);CHKERRQ(ierr);
460   for (p = pStart; p < pEnd; ++p) {
461     PetscInt dof, cdof, off, d;
462     PetscInt adof, aoff;
463 
464     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
465     ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
466     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
467     if (!dof) continue;
468     ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
469     if (adof <= 0) continue;
470     for (d = off; d < off+dof-cdof; ++d) {
471       ierr = PetscSectionGetDof(rootSectionAdj, d, &adof);CHKERRQ(ierr);
472       ierr = PetscSectionGetOffset(rootSectionAdj, d, &aoff);CHKERRQ(ierr);
473       ierr = PetscSortRemoveDupsInt(&adof, &rootAdj[aoff]);CHKERRQ(ierr);
474       ierr = PetscSectionSetDof(rootSectionAdj, d, adof);CHKERRQ(ierr);
475     }
476   }
477   /* Debugging */
478   if (debug) {
479     IS tmp;
480     ierr = PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after compression:\n");CHKERRQ(ierr);
481     ierr = PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
482     ierr = PetscPrintf(comm, "Root adjacency indices after compression\n");CHKERRQ(ierr);
483     ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
484     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
485     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
486   }
487   /* Build adjacency section: Maps global indices to sets of adjacent global indices */
488   ierr = PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd);CHKERRQ(ierr);
489   ierr = PetscSectionCreate(comm, &sectionAdj);CHKERRQ(ierr);
490   ierr = PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd);CHKERRQ(ierr);
491   for (p = pStart; p < pEnd; ++p) {
492     PetscInt  numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof;
493     PetscBool found  = PETSC_TRUE;
494 
495     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
496     ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
497     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
498     ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
499     for (d = 0; d < dof-cdof; ++d) {
500       PetscInt ldof, rdof;
501 
502       ierr = PetscSectionGetDof(leafSectionAdj, off+d, &ldof);CHKERRQ(ierr);
503       ierr = PetscSectionGetDof(rootSectionAdj, off+d, &rdof);CHKERRQ(ierr);
504       if (ldof > 0) {
505         /* We do not own this point */
506       } else if (rdof > 0) {
507         ierr = PetscSectionSetDof(sectionAdj, goff+d, rdof);CHKERRQ(ierr);
508       } else {
509         found = PETSC_FALSE;
510       }
511     }
512     if (found) continue;
513     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
514     ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
515     ierr = DMPlexGetAdjacency(dm, p, &numAdj, &tmpAdj);CHKERRQ(ierr);
516     for (q = 0; q < numAdj; ++q) {
517       const PetscInt padj = tmpAdj[q];
518       PetscInt ndof, ncdof, noff;
519 
520       if ((padj < pStart) || (padj >= pEnd)) continue;
521       ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
522       ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
523       ierr = PetscSectionGetOffset(section, padj, &noff);CHKERRQ(ierr);
524       for (d = goff; d < goff+dof-cdof; ++d) {
525         ierr = PetscSectionAddDof(sectionAdj, d, ndof-ncdof);CHKERRQ(ierr);
526       }
527     }
528     ierr = PetscSectionGetDof(anchorSectionAdj, p, &anDof);CHKERRQ(ierr);
529     if (anDof) {
530       for (d = goff; d < goff+dof-cdof; ++d) {
531         ierr = PetscSectionAddDof(sectionAdj, d, anDof);CHKERRQ(ierr);
532       }
533     }
534   }
535   ierr = PetscSectionSetUp(sectionAdj);CHKERRQ(ierr);
536   if (debug) {
537     ierr = PetscPrintf(comm, "Adjacency Section for Preallocation:\n");CHKERRQ(ierr);
538     ierr = PetscSectionView(sectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
539   }
540   /* Get adjacent indices */
541   ierr = PetscSectionGetStorageSize(sectionAdj, &numCols);CHKERRQ(ierr);
542   ierr = PetscMalloc1(numCols, &cols);CHKERRQ(ierr);
543   for (p = pStart; p < pEnd; ++p) {
544     PetscInt  numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof, anOff;
545     PetscBool found  = PETSC_TRUE;
546 
547     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
548     ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
549     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
550     ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
551     for (d = 0; d < dof-cdof; ++d) {
552       PetscInt ldof, rdof;
553 
554       ierr = PetscSectionGetDof(leafSectionAdj, off+d, &ldof);CHKERRQ(ierr);
555       ierr = PetscSectionGetDof(rootSectionAdj, off+d, &rdof);CHKERRQ(ierr);
556       if (ldof > 0) {
557         /* We do not own this point */
558       } else if (rdof > 0) {
559         PetscInt aoff, roff;
560 
561         ierr = PetscSectionGetOffset(sectionAdj, goff+d, &aoff);CHKERRQ(ierr);
562         ierr = PetscSectionGetOffset(rootSectionAdj, off+d, &roff);CHKERRQ(ierr);
563         ierr = PetscMemcpy(&cols[aoff], &rootAdj[roff], rdof * sizeof(PetscInt));CHKERRQ(ierr);
564       } else {
565         found = PETSC_FALSE;
566       }
567     }
568     if (found) continue;
569     ierr = DMPlexGetAdjacency(dm, p, &numAdj, &tmpAdj);CHKERRQ(ierr);
570     ierr = PetscSectionGetDof(anchorSectionAdj, p, &anDof);CHKERRQ(ierr);
571     ierr = PetscSectionGetOffset(anchorSectionAdj, p, &anOff);CHKERRQ(ierr);
572     for (d = goff; d < goff+dof-cdof; ++d) {
573       PetscInt adof, aoff, i = 0;
574 
575       ierr = PetscSectionGetDof(sectionAdj, d, &adof);CHKERRQ(ierr);
576       ierr = PetscSectionGetOffset(sectionAdj, d, &aoff);CHKERRQ(ierr);
577       for (q = 0; q < numAdj; ++q) {
578         const PetscInt  padj = tmpAdj[q];
579         PetscInt        ndof, ncdof, ngoff, nd;
580         const PetscInt *ncind;
581 
582         /* Adjacent points may not be in the section chart */
583         if ((padj < pStart) || (padj >= pEnd)) continue;
584         ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
585         ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
586         ierr = PetscSectionGetConstraintIndices(section, padj, &ncind);CHKERRQ(ierr);
587         ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
588         for (nd = 0; nd < ndof-ncdof; ++nd, ++i) {
589           cols[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
590         }
591       }
592       for (q = 0; q < anDof; q++, i++) {
593         cols[aoff+i] = anchorAdj[anOff + q];
594       }
595       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);
596     }
597   }
598   ierr = PetscSectionDestroy(&anchorSectionAdj);CHKERRQ(ierr);
599   ierr = PetscSectionDestroy(&leafSectionAdj);CHKERRQ(ierr);
600   ierr = PetscSectionDestroy(&rootSectionAdj);CHKERRQ(ierr);
601   ierr = PetscFree(anchorAdj);CHKERRQ(ierr);
602   ierr = PetscFree(rootAdj);CHKERRQ(ierr);
603   ierr = PetscFree(tmpAdj);CHKERRQ(ierr);
604   /* Debugging */
605   if (debug) {
606     IS tmp;
607     ierr = PetscPrintf(comm, "Column indices\n");CHKERRQ(ierr);
608     ierr = ISCreateGeneral(comm, numCols, cols, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
609     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
610     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
611   }
612   /* Create allocation vectors from adjacency graph */
613   ierr = MatGetLocalSize(A, &locRows, NULL);CHKERRQ(ierr);
614   ierr = PetscLayoutCreate(PetscObjectComm((PetscObject)A), &rLayout);CHKERRQ(ierr);
615   ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr);
616   ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr);
617   ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr);
618   ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr);
619   ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr);
620   /* Only loop over blocks of rows */
621   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);
622   for (r = rStart/bs; r < rEnd/bs; ++r) {
623     const PetscInt row = r*bs;
624     PetscInt       numCols, cStart, c;
625 
626     ierr = PetscSectionGetDof(sectionAdj, row, &numCols);CHKERRQ(ierr);
627     ierr = PetscSectionGetOffset(sectionAdj, row, &cStart);CHKERRQ(ierr);
628     for (c = cStart; c < cStart+numCols; ++c) {
629       if ((cols[c] >= rStart*bs) && (cols[c] < rEnd*bs)) {
630         ++dnz[r-rStart];
631         if (cols[c] >= row) ++dnzu[r-rStart];
632       } else {
633         ++onz[r-rStart];
634         if (cols[c] >= row) ++onzu[r-rStart];
635       }
636     }
637   }
638   if (bs > 1) {
639     for (r = 0; r < locRows/bs; ++r) {
640       dnz[r]  /= bs;
641       onz[r]  /= bs;
642       dnzu[r] /= bs;
643       onzu[r] /= bs;
644     }
645   }
646   /* Set matrix pattern */
647   ierr = MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu);CHKERRQ(ierr);
648   ierr = MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
649   /* Check for symmetric storage */
650   ierr = MatGetType(A, &mtype);CHKERRQ(ierr);
651   ierr = PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);CHKERRQ(ierr);
652   ierr = PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);CHKERRQ(ierr);
653   ierr = PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);CHKERRQ(ierr);
654   if (isSymBlock || isSymSeqBlock || isSymMPIBlock) {ierr = MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);CHKERRQ(ierr);}
655   /* Fill matrix with zeros */
656   if (fillMatrix) {
657     PetscScalar *values;
658     PetscInt     maxRowLen = 0;
659 
660     for (r = rStart; r < rEnd; ++r) {
661       PetscInt len;
662 
663       ierr      = PetscSectionGetDof(sectionAdj, r, &len);CHKERRQ(ierr);
664       maxRowLen = PetscMax(maxRowLen, len);
665     }
666     ierr = PetscCalloc1(maxRowLen, &values);CHKERRQ(ierr);
667     for (r = rStart; r < rEnd; ++r) {
668       PetscInt numCols, cStart;
669 
670       ierr = PetscSectionGetDof(sectionAdj, r, &numCols);CHKERRQ(ierr);
671       ierr = PetscSectionGetOffset(sectionAdj, r, &cStart);CHKERRQ(ierr);
672       ierr = MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);CHKERRQ(ierr);
673     }
674     ierr = PetscFree(values);CHKERRQ(ierr);
675     ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
676     ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
677   }
678   /* restore original useAnchors */
679   ierr = DMPlexSetAdjacencyUseAnchors(dm,useAnchors);CHKERRQ(ierr);
680   ierr = PetscSectionDestroy(&sectionAdj);CHKERRQ(ierr);
681   ierr = PetscFree(cols);CHKERRQ(ierr);
682   ierr = PetscLogEventEnd(DMPLEX_Preallocate,dm,0,0,0);CHKERRQ(ierr);
683   PetscFunctionReturn(0);
684 }
685 
686 #if 0
687 #undef __FUNCT__
688 #define __FUNCT__ "DMPlexPreallocateOperator_2"
689 PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
690 {
691   PetscInt       *tmpClosure,*tmpAdj,*visits;
692   PetscInt        c,cStart,cEnd,pStart,pEnd;
693   PetscErrorCode  ierr;
694 
695   PetscFunctionBegin;
696   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
697   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
698   ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
699 
700   maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1));
701 
702   ierr    = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
703   npoints = pEnd - pStart;
704 
705   ierr = PetscMalloc3(maxClosureSize,&tmpClosure,npoints,&lvisits,npoints,&visits);CHKERRQ(ierr);
706   ierr = PetscMemzero(lvisits,(pEnd-pStart)*sizeof(PetscInt));CHKERRQ(ierr);
707   ierr = PetscMemzero(visits,(pEnd-pStart)*sizeof(PetscInt));CHKERRQ(ierr);
708   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
709   for (c=cStart; c<cEnd; c++) {
710     PetscInt *support = tmpClosure;
711     ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support);CHKERRQ(ierr);
712     for (p=0; p<supportSize; p++) lvisits[support[p]]++;
713   }
714   ierr = PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
715   ierr = PetscSFReduceEnd  (sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
716   ierr = PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits);CHKERRQ(ierr);
717   ierr = PetscSFBcastEnd  (sf,MPIU_INT,visits,lvisits);CHKERRQ(ierr);
718 
719   ierr = PetscSFGetRanks();CHKERRQ(ierr);
720 
721 
722   ierr = PetscMalloc2(maxClosureSize*maxClosureSize,&cellmat,npoints,&owner);CHKERRQ(ierr);
723   for (c=cStart; c<cEnd; c++) {
724     ierr = PetscMemzero(cellmat,maxClosureSize*maxClosureSize*sizeof(PetscInt));CHKERRQ(ierr);
725     /*
726      Depth-first walk of transitive closure.
727      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.
728      This contribution is added to dnz if owning ranks of p and q match, to onz otherwise.
729      */
730   }
731 
732   ierr = PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM);CHKERRQ(ierr);
733   ierr = PetscSFReduceEnd  (sf,MPIU_INT,lonz,onz,MPI_SUM);CHKERRQ(ierr);
734   PetscFunctionReturn(0);
735 }
736 #endif
737