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