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