xref: /petsc/src/dm/impls/plex/plexpreallocate.c (revision f5476a2b443e0eedb22d31cb1b2d85b9b918f73f) 
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 adjacenct 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     ierr = PetscSFGatherBegin(sfAdj, MPIU_INT, adj, rootAdj);CHKERRQ(ierr);
387     ierr = PetscSFGatherEnd(sfAdj, MPIU_INT, adj, rootAdj);CHKERRQ(ierr);
388   }
389   ierr = PetscSFDestroy(&sfAdj);CHKERRQ(ierr);
390   ierr = PetscFree(adj);CHKERRQ(ierr);
391   /* Debugging */
392   if (debug) {
393     IS tmp;
394     ierr = PetscPrintf(comm, "Root adjacency indices after gather\n");CHKERRQ(ierr);
395     ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
396     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
397     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
398   }
399   /* Add in local adjacency indices for owned dofs on interface (roots) */
400   for (p = pStart; p < pEnd; ++p) {
401     PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof, anOff;
402 
403     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
404     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
405     if (!dof) continue;
406     ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
407     if (adof <= 0) continue;
408     ierr = DMPlexGetAdjacency(dm, p, &numAdj, &tmpAdj);CHKERRQ(ierr);
409     ierr = PetscSectionGetDof(anchorSectionAdj, p, &anDof);CHKERRQ(ierr);
410     ierr = PetscSectionGetOffset(anchorSectionAdj, p, &anOff);CHKERRQ(ierr);
411     for (d = off; d < off+dof; ++d) {
412       PetscInt adof, aoff, i;
413 
414       ierr = PetscSectionGetDof(rootSectionAdj, d, &adof);CHKERRQ(ierr);
415       ierr = PetscSectionGetOffset(rootSectionAdj, d, &aoff);CHKERRQ(ierr);
416       i    = adof-1;
417       for (q = 0; q < anDof; q++) {
418         rootAdj[aoff+i] = anchorAdj[anOff+q];
419         --i;
420       }
421       for (q = 0; q < numAdj; ++q) {
422         const PetscInt padj = tmpAdj[q];
423         PetscInt ndof, ncdof, ngoff, nd;
424 
425         if ((padj < pStart) || (padj >= pEnd)) continue;
426         ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
427         ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
428         ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
429         for (nd = 0; nd < ndof-ncdof; ++nd) {
430           rootAdj[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
431           --i;
432         }
433       }
434     }
435   }
436   /* Debugging */
437   if (debug) {
438     IS tmp;
439     ierr = PetscPrintf(comm, "Root adjacency indices\n");CHKERRQ(ierr);
440     ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
441     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
442     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
443   }
444   /* Compress indices */
445   ierr = PetscSectionSetUp(rootSectionAdj);CHKERRQ(ierr);
446   for (p = pStart; p < pEnd; ++p) {
447     PetscInt dof, cdof, off, d;
448     PetscInt adof, aoff;
449 
450     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
451     ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
452     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
453     if (!dof) continue;
454     ierr = PetscSectionGetDof(rootSectionAdj, off, &adof);CHKERRQ(ierr);
455     if (adof <= 0) continue;
456     for (d = off; d < off+dof-cdof; ++d) {
457       ierr = PetscSectionGetDof(rootSectionAdj, d, &adof);CHKERRQ(ierr);
458       ierr = PetscSectionGetOffset(rootSectionAdj, d, &aoff);CHKERRQ(ierr);
459       ierr = PetscSortRemoveDupsInt(&adof, &rootAdj[aoff]);CHKERRQ(ierr);
460       ierr = PetscSectionSetDof(rootSectionAdj, d, adof);CHKERRQ(ierr);
461     }
462   }
463   /* Debugging */
464   if (debug) {
465     IS tmp;
466     ierr = PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after compression:\n");CHKERRQ(ierr);
467     ierr = PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
468     ierr = PetscPrintf(comm, "Root adjacency indices after compression\n");CHKERRQ(ierr);
469     ierr = ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
470     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
471     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
472   }
473   /* Build adjacency section: Maps global indices to sets of adjacent global indices */
474   ierr = PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd);CHKERRQ(ierr);
475   ierr = PetscSectionCreate(comm, &sectionAdj);CHKERRQ(ierr);
476   ierr = PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd);CHKERRQ(ierr);
477   for (p = pStart; p < pEnd; ++p) {
478     PetscInt  numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof;
479     PetscBool found  = PETSC_TRUE;
480 
481     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
482     ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
483     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
484     ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
485     for (d = 0; d < dof-cdof; ++d) {
486       PetscInt ldof, rdof;
487 
488       ierr = PetscSectionGetDof(leafSectionAdj, off+d, &ldof);CHKERRQ(ierr);
489       ierr = PetscSectionGetDof(rootSectionAdj, off+d, &rdof);CHKERRQ(ierr);
490       if (ldof > 0) {
491         /* We do not own this point */
492       } else if (rdof > 0) {
493         ierr = PetscSectionSetDof(sectionAdj, goff+d, rdof);CHKERRQ(ierr);
494       } else {
495         found = PETSC_FALSE;
496       }
497     }
498     if (found) continue;
499     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
500     ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
501     ierr = DMPlexGetAdjacency(dm, p, &numAdj, &tmpAdj);CHKERRQ(ierr);
502     for (q = 0; q < numAdj; ++q) {
503       const PetscInt padj = tmpAdj[q];
504       PetscInt ndof, ncdof, noff;
505 
506       if ((padj < pStart) || (padj >= pEnd)) continue;
507       ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
508       ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
509       ierr = PetscSectionGetOffset(section, padj, &noff);CHKERRQ(ierr);
510       for (d = goff; d < goff+dof-cdof; ++d) {
511         ierr = PetscSectionAddDof(sectionAdj, d, ndof-ncdof);CHKERRQ(ierr);
512       }
513     }
514     ierr = PetscSectionGetDof(anchorSectionAdj, p, &anDof);CHKERRQ(ierr);
515     if (anDof) {
516       for (d = goff; d < goff+dof-cdof; ++d) {
517         ierr = PetscSectionAddDof(sectionAdj, d, anDof);CHKERRQ(ierr);
518       }
519     }
520   }
521   ierr = PetscSectionSetUp(sectionAdj);CHKERRQ(ierr);
522   if (debug) {
523     ierr = PetscPrintf(comm, "Adjacency Section for Preallocation:\n");CHKERRQ(ierr);
524     ierr = PetscSectionView(sectionAdj, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
525   }
526   /* Get adjacent indices */
527   ierr = PetscSectionGetStorageSize(sectionAdj, &numCols);CHKERRQ(ierr);
528   ierr = PetscMalloc1(numCols, &cols);CHKERRQ(ierr);
529   for (p = pStart; p < pEnd; ++p) {
530     PetscInt  numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof, anOff;
531     PetscBool found  = PETSC_TRUE;
532 
533     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
534     ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr);
535     ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
536     ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
537     for (d = 0; d < dof-cdof; ++d) {
538       PetscInt ldof, rdof;
539 
540       ierr = PetscSectionGetDof(leafSectionAdj, off+d, &ldof);CHKERRQ(ierr);
541       ierr = PetscSectionGetDof(rootSectionAdj, off+d, &rdof);CHKERRQ(ierr);
542       if (ldof > 0) {
543         /* We do not own this point */
544       } else if (rdof > 0) {
545         PetscInt aoff, roff;
546 
547         ierr = PetscSectionGetOffset(sectionAdj, goff+d, &aoff);CHKERRQ(ierr);
548         ierr = PetscSectionGetOffset(rootSectionAdj, off+d, &roff);CHKERRQ(ierr);
549         ierr = PetscMemcpy(&cols[aoff], &rootAdj[roff], rdof * sizeof(PetscInt));CHKERRQ(ierr);
550       } else {
551         found = PETSC_FALSE;
552       }
553     }
554     if (found) continue;
555     ierr = DMPlexGetAdjacency(dm, p, &numAdj, &tmpAdj);CHKERRQ(ierr);
556     ierr = PetscSectionGetDof(anchorSectionAdj, p, &anDof);CHKERRQ(ierr);
557     ierr = PetscSectionGetOffset(anchorSectionAdj, p, &anOff);CHKERRQ(ierr);
558     for (d = goff; d < goff+dof-cdof; ++d) {
559       PetscInt adof, aoff, i = 0;
560 
561       ierr = PetscSectionGetDof(sectionAdj, d, &adof);CHKERRQ(ierr);
562       ierr = PetscSectionGetOffset(sectionAdj, d, &aoff);CHKERRQ(ierr);
563       for (q = 0; q < numAdj; ++q) {
564         const PetscInt  padj = tmpAdj[q];
565         PetscInt        ndof, ncdof, ngoff, nd;
566         const PetscInt *ncind;
567 
568         /* Adjacent points may not be in the section chart */
569         if ((padj < pStart) || (padj >= pEnd)) continue;
570         ierr = PetscSectionGetDof(section, padj, &ndof);CHKERRQ(ierr);
571         ierr = PetscSectionGetConstraintDof(section, padj, &ncdof);CHKERRQ(ierr);
572         ierr = PetscSectionGetConstraintIndices(section, padj, &ncind);CHKERRQ(ierr);
573         ierr = PetscSectionGetOffset(sectionGlobal, padj, &ngoff);CHKERRQ(ierr);
574         for (nd = 0; nd < ndof-ncdof; ++nd, ++i) {
575           cols[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
576         }
577       }
578       for (q = 0; q < anDof; q++, i++) {
579         cols[aoff+i] = anchorAdj[anOff + q];
580       }
581       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);
582     }
583   }
584   ierr = PetscSectionDestroy(&anchorSectionAdj);CHKERRQ(ierr);
585   ierr = PetscSectionDestroy(&leafSectionAdj);CHKERRQ(ierr);
586   ierr = PetscSectionDestroy(&rootSectionAdj);CHKERRQ(ierr);
587   ierr = PetscFree(anchorAdj);CHKERRQ(ierr);
588   ierr = PetscFree(rootAdj);CHKERRQ(ierr);
589   ierr = PetscFree(tmpAdj);CHKERRQ(ierr);
590   /* Debugging */
591   if (debug) {
592     IS tmp;
593     ierr = PetscPrintf(comm, "Column indices\n");CHKERRQ(ierr);
594     ierr = ISCreateGeneral(comm, numCols, cols, PETSC_USE_POINTER, &tmp);CHKERRQ(ierr);
595     ierr = ISView(tmp, NULL);CHKERRQ(ierr);
596     ierr = ISDestroy(&tmp);CHKERRQ(ierr);
597   }
598   /* Create allocation vectors from adjacency graph */
599   ierr = MatGetLocalSize(A, &locRows, NULL);CHKERRQ(ierr);
600   ierr = PetscLayoutCreate(PetscObjectComm((PetscObject)A), &rLayout);CHKERRQ(ierr);
601   ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr);
602   ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr);
603   ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr);
604   ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr);
605   ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr);
606   /* Only loop over blocks of rows */
607   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);
608   for (r = rStart/bs; r < rEnd/bs; ++r) {
609     const PetscInt row = r*bs;
610     PetscInt       numCols, cStart, c;
611 
612     ierr = PetscSectionGetDof(sectionAdj, row, &numCols);CHKERRQ(ierr);
613     ierr = PetscSectionGetOffset(sectionAdj, row, &cStart);CHKERRQ(ierr);
614     for (c = cStart; c < cStart+numCols; ++c) {
615       if ((cols[c] >= rStart*bs) && (cols[c] < rEnd*bs)) {
616         ++dnz[r-rStart];
617         if (cols[c] >= row) ++dnzu[r-rStart];
618       } else {
619         ++onz[r-rStart];
620         if (cols[c] >= row) ++onzu[r-rStart];
621       }
622     }
623   }
624   if (bs > 1) {
625     for (r = 0; r < locRows/bs; ++r) {
626       dnz[r]  /= bs;
627       onz[r]  /= bs;
628       dnzu[r] /= bs;
629       onzu[r] /= bs;
630     }
631   }
632   /* Set matrix pattern */
633   ierr = MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu);CHKERRQ(ierr);
634   ierr = MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
635   /* Check for symmetric storage */
636   ierr = MatGetType(A, &mtype);CHKERRQ(ierr);
637   ierr = PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);CHKERRQ(ierr);
638   ierr = PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);CHKERRQ(ierr);
639   ierr = PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);CHKERRQ(ierr);
640   if (isSymBlock || isSymSeqBlock || isSymMPIBlock) {ierr = MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);CHKERRQ(ierr);}
641   /* Fill matrix with zeros */
642   if (fillMatrix) {
643     PetscScalar *values;
644     PetscInt     maxRowLen = 0;
645 
646     for (r = rStart; r < rEnd; ++r) {
647       PetscInt len;
648 
649       ierr      = PetscSectionGetDof(sectionAdj, r, &len);CHKERRQ(ierr);
650       maxRowLen = PetscMax(maxRowLen, len);
651     }
652     ierr = PetscCalloc1(maxRowLen, &values);CHKERRQ(ierr);
653     for (r = rStart; r < rEnd; ++r) {
654       PetscInt numCols, cStart;
655 
656       ierr = PetscSectionGetDof(sectionAdj, r, &numCols);CHKERRQ(ierr);
657       ierr = PetscSectionGetOffset(sectionAdj, r, &cStart);CHKERRQ(ierr);
658       ierr = MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);CHKERRQ(ierr);
659     }
660     ierr = PetscFree(values);CHKERRQ(ierr);
661     ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
662     ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
663   }
664   /* restore original useAnchors */
665   ierr = DMPlexSetAdjacencyUseAnchors(dm,useAnchors);CHKERRQ(ierr);
666   ierr = PetscSectionDestroy(&sectionAdj);CHKERRQ(ierr);
667   ierr = PetscFree(cols);CHKERRQ(ierr);
668   ierr = PetscLogEventEnd(DMPLEX_Preallocate,dm,0,0,0);CHKERRQ(ierr);
669   PetscFunctionReturn(0);
670 }
671 
672 #if 0
673 #undef __FUNCT__
674 #define __FUNCT__ "DMPlexPreallocateOperator_2"
675 PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
676 {
677   PetscInt       *tmpClosure,*tmpAdj,*visits;
678   PetscInt        c,cStart,cEnd,pStart,pEnd;
679   PetscErrorCode  ierr;
680 
681   PetscFunctionBegin;
682   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
683   ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
684   ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
685 
686   maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1));
687 
688   ierr    = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
689   npoints = pEnd - pStart;
690 
691   ierr = PetscMalloc3(maxClosureSize,&tmpClosure,npoints,&lvisits,npoints,&visits);CHKERRQ(ierr);
692   ierr = PetscMemzero(lvisits,(pEnd-pStart)*sizeof(PetscInt));CHKERRQ(ierr);
693   ierr = PetscMemzero(visits,(pEnd-pStart)*sizeof(PetscInt));CHKERRQ(ierr);
694   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
695   for (c=cStart; c<cEnd; c++) {
696     PetscInt *support = tmpClosure;
697     ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support);CHKERRQ(ierr);
698     for (p=0; p<supportSize; p++) lvisits[support[p]]++;
699   }
700   ierr = PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
701   ierr = PetscSFReduceEnd  (sf,MPIU_INT,lvisits,visits,MPI_SUM);CHKERRQ(ierr);
702   ierr = PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits);CHKERRQ(ierr);
703   ierr = PetscSFBcastEnd  (sf,MPIU_INT,visits,lvisits);CHKERRQ(ierr);
704 
705   ierr = PetscSFGetRanks();CHKERRQ(ierr);
706 
707 
708   ierr = PetscMalloc2(maxClosureSize*maxClosureSize,&cellmat,npoints,&owner);CHKERRQ(ierr);
709   for (c=cStart; c<cEnd; c++) {
710     ierr = PetscMemzero(cellmat,maxClosureSize*maxClosureSize*sizeof(PetscInt));CHKERRQ(ierr);
711     /*
712      Depth-first walk of transitive closure.
713      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.
714      This contribution is added to dnz if owning ranks of p and q match, to onz otherwise.
715      */
716   }
717 
718   ierr = PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM);CHKERRQ(ierr);
719   ierr = PetscSFReduceEnd  (sf,MPIU_INT,lonz,onz,MPI_SUM);CHKERRQ(ierr);
720   PetscFunctionReturn(0);
721 }
722 #endif
723