Actual source code: plexpreallocate.c

  1: #include <petsc/private/dmpleximpl.h>
  2: #include <petsc/private/isimpl.h>
  3: #include <petscsf.h>
  4: #include <petscds.h>

  6: /* get adjacencies due to point-to-point constraints that can't be found with DMPlexGetAdjacency() */
  7: static PetscErrorCode DMPlexComputeAnchorAdjacencies(DM dm, PetscBool useCone, PetscBool useClosure, PetscSection *anchorSectionAdj, PetscInt *anchorAdj[])
  8: {
  9:   PetscInt     pStart, pEnd;
 10:   PetscSection section, sectionGlobal, adjSec, aSec;
 11:   IS           aIS;

 13:   DMGetLocalSection(dm, &section);
 14:   DMGetGlobalSection(dm, &sectionGlobal);
 15:   PetscSectionCreate(PetscObjectComm((PetscObject)section), &adjSec);
 16:   PetscSectionGetChart(section, &pStart, &pEnd);
 17:   PetscSectionSetChart(adjSec, pStart, pEnd);

 19:   DMPlexGetAnchors(dm, &aSec, &aIS);
 20:   if (aSec) {
 21:     const PetscInt *anchors;
 22:     PetscInt        p, q, a, aSize, *offsets, aStart, aEnd, *inverse, iSize, *adj, adjSize;
 23:     PetscInt       *tmpAdjP = NULL, *tmpAdjQ = NULL;
 24:     PetscSection    inverseSec;

 26:     /* invert the constraint-to-anchor map */
 27:     PetscSectionCreate(PetscObjectComm((PetscObject)aSec), &inverseSec);
 28:     PetscSectionSetChart(inverseSec, pStart, pEnd);
 29:     ISGetLocalSize(aIS, &aSize);
 30:     ISGetIndices(aIS, &anchors);

 32:     for (p = 0; p < aSize; p++) {
 33:       PetscInt a = anchors[p];

 35:       PetscSectionAddDof(inverseSec, a, 1);
 36:     }
 37:     PetscSectionSetUp(inverseSec);
 38:     PetscSectionGetStorageSize(inverseSec, &iSize);
 39:     PetscMalloc1(iSize, &inverse);
 40:     PetscCalloc1(pEnd - pStart, &offsets);
 41:     PetscSectionGetChart(aSec, &aStart, &aEnd);
 42:     for (p = aStart; p < aEnd; p++) {
 43:       PetscInt dof, off;

 45:       PetscSectionGetDof(aSec, p, &dof);
 46:       PetscSectionGetOffset(aSec, p, &off);

 48:       for (q = 0; q < dof; q++) {
 49:         PetscInt iOff;

 51:         a = anchors[off + q];
 52:         PetscSectionGetOffset(inverseSec, a, &iOff);
 53:         inverse[iOff + offsets[a - pStart]++] = p;
 54:       }
 55:     }
 56:     ISRestoreIndices(aIS, &anchors);
 57:     PetscFree(offsets);

 59:     /* construct anchorAdj and adjSec
 60:      *
 61:      * loop over anchors:
 62:      *   construct anchor adjacency
 63:      *   loop over constrained:
 64:      *     construct constrained adjacency
 65:      *     if not in anchor adjacency, add to dofs
 66:      * setup adjSec, allocate anchorAdj
 67:      * loop over anchors:
 68:      *   construct anchor adjacency
 69:      *   loop over constrained:
 70:      *     construct constrained adjacency
 71:      *     if not in anchor adjacency
 72:      *       if not already in list, put in list
 73:      *   sort, unique, reduce dof count
 74:      * optional: compactify
 75:      */
 76:     for (p = pStart; p < pEnd; p++) {
 77:       PetscInt iDof, iOff, i, r, s, numAdjP = PETSC_DETERMINE;

 79:       PetscSectionGetDof(inverseSec, p, &iDof);
 80:       if (!iDof) continue;
 81:       PetscSectionGetOffset(inverseSec, p, &iOff);
 82:       DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, PETSC_TRUE, &numAdjP, &tmpAdjP);
 83:       for (i = 0; i < iDof; i++) {
 84:         PetscInt iNew = 0, qAdj, qAdjDof, qAdjCDof, numAdjQ = PETSC_DETERMINE;

 86:         q = inverse[iOff + i];
 87:         DMPlexGetAdjacency_Internal(dm, q, useCone, useClosure, PETSC_TRUE, &numAdjQ, &tmpAdjQ);
 88:         for (r = 0; r < numAdjQ; r++) {
 89:           qAdj = tmpAdjQ[r];
 90:           if ((qAdj < pStart) || (qAdj >= pEnd)) continue;
 91:           for (s = 0; s < numAdjP; s++) {
 92:             if (qAdj == tmpAdjP[s]) break;
 93:           }
 94:           if (s < numAdjP) continue;
 95:           PetscSectionGetDof(section, qAdj, &qAdjDof);
 96:           PetscSectionGetConstraintDof(section, qAdj, &qAdjCDof);
 97:           iNew += qAdjDof - qAdjCDof;
 98:         }
 99:         PetscSectionAddDof(adjSec, p, iNew);
100:       }
101:     }

103:     PetscSectionSetUp(adjSec);
104:     PetscSectionGetStorageSize(adjSec, &adjSize);
105:     PetscMalloc1(adjSize, &adj);

107:     for (p = pStart; p < pEnd; p++) {
108:       PetscInt iDof, iOff, i, r, s, aOff, aOffOrig, aDof, numAdjP = PETSC_DETERMINE;

110:       PetscSectionGetDof(inverseSec, p, &iDof);
111:       if (!iDof) continue;
112:       PetscSectionGetOffset(inverseSec, p, &iOff);
113:       DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, PETSC_TRUE, &numAdjP, &tmpAdjP);
114:       PetscSectionGetDof(adjSec, p, &aDof);
115:       PetscSectionGetOffset(adjSec, p, &aOff);
116:       aOffOrig = aOff;
117:       for (i = 0; i < iDof; i++) {
118:         PetscInt qAdj, qAdjDof, qAdjCDof, qAdjOff, nd, numAdjQ = PETSC_DETERMINE;

120:         q = inverse[iOff + i];
121:         DMPlexGetAdjacency_Internal(dm, q, useCone, useClosure, PETSC_TRUE, &numAdjQ, &tmpAdjQ);
122:         for (r = 0; r < numAdjQ; r++) {
123:           qAdj = tmpAdjQ[r];
124:           if ((qAdj < pStart) || (qAdj >= pEnd)) continue;
125:           for (s = 0; s < numAdjP; s++) {
126:             if (qAdj == tmpAdjP[s]) break;
127:           }
128:           if (s < numAdjP) continue;
129:           PetscSectionGetDof(section, qAdj, &qAdjDof);
130:           PetscSectionGetConstraintDof(section, qAdj, &qAdjCDof);
131:           PetscSectionGetOffset(sectionGlobal, qAdj, &qAdjOff);
132:           for (nd = 0; nd < qAdjDof - qAdjCDof; ++nd) adj[aOff++] = (qAdjOff < 0 ? -(qAdjOff + 1) : qAdjOff) + nd;
133:         }
134:       }
135:       PetscSortRemoveDupsInt(&aDof, &adj[aOffOrig]);
136:       PetscSectionSetDof(adjSec, p, aDof);
137:     }
138:     *anchorAdj = adj;

140:     /* clean up */
141:     PetscSectionDestroy(&inverseSec);
142:     PetscFree(inverse);
143:     PetscFree(tmpAdjP);
144:     PetscFree(tmpAdjQ);
145:   } else {
146:     *anchorAdj = NULL;
147:     PetscSectionSetUp(adjSec);
148:   }
149:   *anchorSectionAdj = adjSec;
150:   return 0;
151: }

153: static PetscErrorCode DMPlexCreateAdjacencySection_Static(DM dm, PetscInt bs, PetscSF sfDof, PetscBool useCone, PetscBool useClosure, PetscBool useAnchors, PetscSection *sA, PetscInt **colIdx)
154: {
155:   MPI_Comm           comm;
156:   PetscMPIInt        size;
157:   PetscBool          doCommLocal, doComm, debug = PETSC_FALSE;
158:   PetscSF            sf, sfAdj;
159:   PetscSection       section, sectionGlobal, leafSectionAdj, rootSectionAdj, sectionAdj, anchorSectionAdj;
160:   PetscInt           nroots, nleaves, l, p, r;
161:   const PetscInt    *leaves;
162:   const PetscSFNode *remotes;
163:   PetscInt           dim, pStart, pEnd, numDof, globalOffStart, globalOffEnd, numCols;
164:   PetscInt          *tmpAdj = NULL, *adj, *rootAdj, *anchorAdj = NULL, *cols, *remoteOffsets;
165:   PetscInt           adjSize;

167:   PetscObjectGetComm((PetscObject)dm, &comm);
168:   PetscOptionsGetBool(NULL, NULL, "-dm_view_preallocation", &debug, NULL);
169:   MPI_Comm_size(comm, &size);
170:   DMGetDimension(dm, &dim);
171:   DMGetPointSF(dm, &sf);
172:   DMGetLocalSection(dm, &section);
173:   DMGetGlobalSection(dm, &sectionGlobal);
174:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
175:   doCommLocal = (size > 1) && (nroots >= 0) ? PETSC_TRUE : PETSC_FALSE;
176:   MPIU_Allreduce(&doCommLocal, &doComm, 1, MPIU_BOOL, MPI_LAND, comm);
177:   /* Create section for dof adjacency (dof ==> # adj dof) */
178:   PetscSectionGetChart(section, &pStart, &pEnd);
179:   PetscSectionGetStorageSize(section, &numDof);
180:   PetscSectionCreate(comm, &leafSectionAdj);
181:   PetscSectionSetChart(leafSectionAdj, 0, numDof);
182:   PetscSectionCreate(comm, &rootSectionAdj);
183:   PetscSectionSetChart(rootSectionAdj, 0, numDof);
184:   /*   Fill in the ghost dofs on the interface */
185:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, &remotes);
186:   /*
187:    section        - maps points to (# dofs, local dofs)
188:    sectionGlobal  - maps points to (# dofs, global dofs)
189:    leafSectionAdj - maps unowned local dofs to # adj dofs
190:    rootSectionAdj - maps   owned local dofs to # adj dofs
191:    adj            - adj global dofs indexed by leafSectionAdj
192:    rootAdj        - adj global dofs indexed by rootSectionAdj
193:    sf    - describes shared points across procs
194:    sfDof - describes shared dofs across procs
195:    sfAdj - describes shared adjacent dofs across procs
196:    ** The bootstrapping process involves six rounds with similar structure of visiting neighbors of each point.
197:   (0). If there are point-to-point constraints, add the adjacencies of constrained points to anchors in anchorAdj
198:        (This is done in DMPlexComputeAnchorAdjacencies())
199:     1. Visit unowned points on interface, count adjacencies placing in leafSectionAdj
200:        Reduce those counts to rootSectionAdj (now redundantly counting some interface points)
201:     2. Visit owned points on interface, count adjacencies placing in rootSectionAdj
202:        Create sfAdj connecting rootSectionAdj and leafSectionAdj
203:     3. Visit unowned points on interface, write adjacencies to adj
204:        Gather adj to rootAdj (note that there is redundancy in rootAdj when multiple procs find the same adjacencies)
205:     4. Visit owned points on interface, write adjacencies to rootAdj
206:        Remove redundancy in rootAdj
207:    ** The last two traversals use transitive closure
208:     5. Visit all owned points in the subdomain, count dofs for each point (sectionAdj)
209:        Allocate memory addressed by sectionAdj (cols)
210:     6. Visit all owned points in the subdomain, insert dof adjacencies into cols
211:    ** Knowing all the column adjacencies, check ownership and sum into dnz and onz
212:   */
213:   DMPlexComputeAnchorAdjacencies(dm, useCone, useClosure, &anchorSectionAdj, &anchorAdj);
214:   for (l = 0; l < nleaves; ++l) {
215:     PetscInt dof, off, d, q, anDof;
216:     PetscInt p = leaves[l], numAdj = PETSC_DETERMINE;

218:     if ((p < pStart) || (p >= pEnd)) continue;
219:     PetscSectionGetDof(section, p, &dof);
220:     PetscSectionGetOffset(section, p, &off);
221:     DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
222:     for (q = 0; q < numAdj; ++q) {
223:       const PetscInt padj = tmpAdj[q];
224:       PetscInt       ndof, ncdof;

226:       if ((padj < pStart) || (padj >= pEnd)) continue;
227:       PetscSectionGetDof(section, padj, &ndof);
228:       PetscSectionGetConstraintDof(section, padj, &ncdof);
229:       for (d = off; d < off + dof; ++d) PetscSectionAddDof(leafSectionAdj, d, ndof - ncdof);
230:     }
231:     PetscSectionGetDof(anchorSectionAdj, p, &anDof);
232:     if (anDof) {
233:       for (d = off; d < off + dof; ++d) PetscSectionAddDof(leafSectionAdj, d, anDof);
234:     }
235:   }
236:   PetscSectionSetUp(leafSectionAdj);
237:   if (debug) {
238:     PetscPrintf(comm, "Adjacency Section for Preallocation on Leaves:\n");
239:     PetscSectionView(leafSectionAdj, NULL);
240:   }
241:   /* Get maximum remote adjacency sizes for owned dofs on interface (roots) */
242:   if (doComm) {
243:     PetscSFReduceBegin(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);
244:     PetscSFReduceEnd(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);
245:     PetscSectionInvalidateMaxDof_Internal(rootSectionAdj);
246:   }
247:   if (debug) {
248:     PetscPrintf(comm, "Adjancency Section for Preallocation on Roots:\n");
249:     PetscSectionView(rootSectionAdj, NULL);
250:   }
251:   /* Add in local adjacency sizes for owned dofs on interface (roots) */
252:   for (p = pStart; p < pEnd; ++p) {
253:     PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof;

255:     PetscSectionGetDof(section, p, &dof);
256:     PetscSectionGetOffset(section, p, &off);
257:     if (!dof) continue;
258:     PetscSectionGetDof(rootSectionAdj, off, &adof);
259:     if (adof <= 0) continue;
260:     DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
261:     for (q = 0; q < numAdj; ++q) {
262:       const PetscInt padj = tmpAdj[q];
263:       PetscInt       ndof, ncdof;

265:       if ((padj < pStart) || (padj >= pEnd)) continue;
266:       PetscSectionGetDof(section, padj, &ndof);
267:       PetscSectionGetConstraintDof(section, padj, &ncdof);
268:       for (d = off; d < off + dof; ++d) PetscSectionAddDof(rootSectionAdj, d, ndof - ncdof);
269:     }
270:     PetscSectionGetDof(anchorSectionAdj, p, &anDof);
271:     if (anDof) {
272:       for (d = off; d < off + dof; ++d) PetscSectionAddDof(rootSectionAdj, d, anDof);
273:     }
274:   }
275:   PetscSectionSetUp(rootSectionAdj);
276:   if (debug) {
277:     PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after local additions:\n");
278:     PetscSectionView(rootSectionAdj, NULL);
279:   }
280:   /* Create adj SF based on dof SF */
281:   PetscSFCreateRemoteOffsets(sfDof, rootSectionAdj, leafSectionAdj, &remoteOffsets);
282:   PetscSFCreateSectionSF(sfDof, rootSectionAdj, remoteOffsets, leafSectionAdj, &sfAdj);
283:   PetscFree(remoteOffsets);
284:   if (debug && size > 1) {
285:     PetscPrintf(comm, "Adjacency SF for Preallocation:\n");
286:     PetscSFView(sfAdj, NULL);
287:   }
288:   /* Create leaf adjacency */
289:   PetscSectionSetUp(leafSectionAdj);
290:   PetscSectionGetStorageSize(leafSectionAdj, &adjSize);
291:   PetscCalloc1(adjSize, &adj);
292:   for (l = 0; l < nleaves; ++l) {
293:     PetscInt dof, off, d, q, anDof, anOff;
294:     PetscInt p = leaves[l], numAdj = PETSC_DETERMINE;

296:     if ((p < pStart) || (p >= pEnd)) continue;
297:     PetscSectionGetDof(section, p, &dof);
298:     PetscSectionGetOffset(section, p, &off);
299:     DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
300:     PetscSectionGetDof(anchorSectionAdj, p, &anDof);
301:     PetscSectionGetOffset(anchorSectionAdj, p, &anOff);
302:     for (d = off; d < off + dof; ++d) {
303:       PetscInt aoff, i = 0;

305:       PetscSectionGetOffset(leafSectionAdj, d, &aoff);
306:       for (q = 0; q < numAdj; ++q) {
307:         const PetscInt padj = tmpAdj[q];
308:         PetscInt       ndof, ncdof, ngoff, nd;

310:         if ((padj < pStart) || (padj >= pEnd)) continue;
311:         PetscSectionGetDof(section, padj, &ndof);
312:         PetscSectionGetConstraintDof(section, padj, &ncdof);
313:         PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
314:         for (nd = 0; nd < ndof - ncdof; ++nd) {
315:           adj[aoff + i] = (ngoff < 0 ? -(ngoff + 1) : ngoff) + nd;
316:           ++i;
317:         }
318:       }
319:       for (q = 0; q < anDof; q++) {
320:         adj[aoff + i] = anchorAdj[anOff + q];
321:         ++i;
322:       }
323:     }
324:   }
325:   /* Debugging */
326:   if (debug) {
327:     IS tmp;
328:     PetscPrintf(comm, "Leaf adjacency indices\n");
329:     ISCreateGeneral(comm, adjSize, adj, PETSC_USE_POINTER, &tmp);
330:     ISView(tmp, NULL);
331:     ISDestroy(&tmp);
332:   }
333:   /* Gather adjacent indices to root */
334:   PetscSectionGetStorageSize(rootSectionAdj, &adjSize);
335:   PetscMalloc1(adjSize, &rootAdj);
336:   for (r = 0; r < adjSize; ++r) rootAdj[r] = -1;
337:   if (doComm) {
338:     const PetscInt *indegree;
339:     PetscInt       *remoteadj, radjsize = 0;

341:     PetscSFComputeDegreeBegin(sfAdj, &indegree);
342:     PetscSFComputeDegreeEnd(sfAdj, &indegree);
343:     for (p = 0; p < adjSize; ++p) radjsize += indegree[p];
344:     PetscMalloc1(radjsize, &remoteadj);
345:     PetscSFGatherBegin(sfAdj, MPIU_INT, adj, remoteadj);
346:     PetscSFGatherEnd(sfAdj, MPIU_INT, adj, remoteadj);
347:     for (p = 0, l = 0, r = 0; p < adjSize; ++p, l = PetscMax(p, l + indegree[p - 1])) {
348:       PetscInt s;
349:       for (s = 0; s < indegree[p]; ++s, ++r) rootAdj[l + s] = remoteadj[r];
350:     }
353:     PetscFree(remoteadj);
354:   }
355:   PetscSFDestroy(&sfAdj);
356:   PetscFree(adj);
357:   /* Debugging */
358:   if (debug) {
359:     IS tmp;
360:     PetscPrintf(comm, "Root adjacency indices after gather\n");
361:     ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
362:     ISView(tmp, NULL);
363:     ISDestroy(&tmp);
364:   }
365:   /* Add in local adjacency indices for owned dofs on interface (roots) */
366:   for (p = pStart; p < pEnd; ++p) {
367:     PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof, anOff;

369:     PetscSectionGetDof(section, p, &dof);
370:     PetscSectionGetOffset(section, p, &off);
371:     if (!dof) continue;
372:     PetscSectionGetDof(rootSectionAdj, off, &adof);
373:     if (adof <= 0) continue;
374:     DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
375:     PetscSectionGetDof(anchorSectionAdj, p, &anDof);
376:     PetscSectionGetOffset(anchorSectionAdj, p, &anOff);
377:     for (d = off; d < off + dof; ++d) {
378:       PetscInt adof, aoff, i;

380:       PetscSectionGetDof(rootSectionAdj, d, &adof);
381:       PetscSectionGetOffset(rootSectionAdj, d, &aoff);
382:       i = adof - 1;
383:       for (q = 0; q < anDof; q++) {
384:         rootAdj[aoff + i] = anchorAdj[anOff + q];
385:         --i;
386:       }
387:       for (q = 0; q < numAdj; ++q) {
388:         const PetscInt padj = tmpAdj[q];
389:         PetscInt       ndof, ncdof, ngoff, nd;

391:         if ((padj < pStart) || (padj >= pEnd)) continue;
392:         PetscSectionGetDof(section, padj, &ndof);
393:         PetscSectionGetConstraintDof(section, padj, &ncdof);
394:         PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
395:         for (nd = 0; nd < ndof - ncdof; ++nd) {
396:           rootAdj[aoff + i] = ngoff < 0 ? -(ngoff + 1) + nd : ngoff + nd;
397:           --i;
398:         }
399:       }
400:     }
401:   }
402:   /* Debugging */
403:   if (debug) {
404:     IS tmp;
405:     PetscPrintf(comm, "Root adjacency indices\n");
406:     ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
407:     ISView(tmp, NULL);
408:     ISDestroy(&tmp);
409:   }
410:   /* Compress indices */
411:   PetscSectionSetUp(rootSectionAdj);
412:   for (p = pStart; p < pEnd; ++p) {
413:     PetscInt dof, cdof, off, d;
414:     PetscInt adof, aoff;

416:     PetscSectionGetDof(section, p, &dof);
417:     PetscSectionGetConstraintDof(section, p, &cdof);
418:     PetscSectionGetOffset(section, p, &off);
419:     if (!dof) continue;
420:     PetscSectionGetDof(rootSectionAdj, off, &adof);
421:     if (adof <= 0) continue;
422:     for (d = off; d < off + dof - cdof; ++d) {
423:       PetscSectionGetDof(rootSectionAdj, d, &adof);
424:       PetscSectionGetOffset(rootSectionAdj, d, &aoff);
425:       PetscSortRemoveDupsInt(&adof, &rootAdj[aoff]);
426:       PetscSectionSetDof(rootSectionAdj, d, adof);
427:     }
428:   }
429:   /* Debugging */
430:   if (debug) {
431:     IS tmp;
432:     PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after compression:\n");
433:     PetscSectionView(rootSectionAdj, NULL);
434:     PetscPrintf(comm, "Root adjacency indices after compression\n");
435:     ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
436:     ISView(tmp, NULL);
437:     ISDestroy(&tmp);
438:   }
439:   /* Build adjacency section: Maps global indices to sets of adjacent global indices */
440:   PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd);
441:   PetscSectionCreate(comm, &sectionAdj);
442:   PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd);
443:   for (p = pStart; p < pEnd; ++p) {
444:     PetscInt  numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof;
445:     PetscBool found  = PETSC_TRUE;

447:     PetscSectionGetDof(section, p, &dof);
448:     PetscSectionGetConstraintDof(section, p, &cdof);
449:     PetscSectionGetOffset(section, p, &off);
450:     PetscSectionGetOffset(sectionGlobal, p, &goff);
451:     for (d = 0; d < dof - cdof; ++d) {
452:       PetscInt ldof, rdof;

454:       PetscSectionGetDof(leafSectionAdj, off + d, &ldof);
455:       PetscSectionGetDof(rootSectionAdj, off + d, &rdof);
456:       if (ldof > 0) {
457:         /* We do not own this point */
458:       } else if (rdof > 0) {
459:         PetscSectionSetDof(sectionAdj, goff + d, rdof);
460:       } else {
461:         found = PETSC_FALSE;
462:       }
463:     }
464:     if (found) continue;
465:     PetscSectionGetDof(section, p, &dof);
466:     PetscSectionGetOffset(sectionGlobal, p, &goff);
467:     DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
468:     for (q = 0; q < numAdj; ++q) {
469:       const PetscInt padj = tmpAdj[q];
470:       PetscInt       ndof, ncdof, noff;

472:       if ((padj < pStart) || (padj >= pEnd)) continue;
473:       PetscSectionGetDof(section, padj, &ndof);
474:       PetscSectionGetConstraintDof(section, padj, &ncdof);
475:       PetscSectionGetOffset(section, padj, &noff);
476:       for (d = goff; d < goff + dof - cdof; ++d) PetscSectionAddDof(sectionAdj, d, ndof - ncdof);
477:     }
478:     PetscSectionGetDof(anchorSectionAdj, p, &anDof);
479:     if (anDof) {
480:       for (d = goff; d < goff + dof - cdof; ++d) PetscSectionAddDof(sectionAdj, d, anDof);
481:     }
482:   }
483:   PetscSectionSetUp(sectionAdj);
484:   if (debug) {
485:     PetscPrintf(comm, "Adjacency Section for Preallocation:\n");
486:     PetscSectionView(sectionAdj, NULL);
487:   }
488:   /* Get adjacent indices */
489:   PetscSectionGetStorageSize(sectionAdj, &numCols);
490:   PetscMalloc1(numCols, &cols);
491:   for (p = pStart; p < pEnd; ++p) {
492:     PetscInt  numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof, anOff;
493:     PetscBool found  = PETSC_TRUE;

495:     PetscSectionGetDof(section, p, &dof);
496:     PetscSectionGetConstraintDof(section, p, &cdof);
497:     PetscSectionGetOffset(section, p, &off);
498:     PetscSectionGetOffset(sectionGlobal, p, &goff);
499:     for (d = 0; d < dof - cdof; ++d) {
500:       PetscInt ldof, rdof;

502:       PetscSectionGetDof(leafSectionAdj, off + d, &ldof);
503:       PetscSectionGetDof(rootSectionAdj, off + d, &rdof);
504:       if (ldof > 0) {
505:         /* We do not own this point */
506:       } else if (rdof > 0) {
507:         PetscInt aoff, roff;

509:         PetscSectionGetOffset(sectionAdj, goff + d, &aoff);
510:         PetscSectionGetOffset(rootSectionAdj, off + d, &roff);
511:         PetscArraycpy(&cols[aoff], &rootAdj[roff], rdof);
512:       } else {
513:         found = PETSC_FALSE;
514:       }
515:     }
516:     if (found) continue;
517:     DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
518:     PetscSectionGetDof(anchorSectionAdj, p, &anDof);
519:     PetscSectionGetOffset(anchorSectionAdj, p, &anOff);
520:     for (d = goff; d < goff + dof - cdof; ++d) {
521:       PetscInt adof, aoff, i = 0;

523:       PetscSectionGetDof(sectionAdj, d, &adof);
524:       PetscSectionGetOffset(sectionAdj, d, &aoff);
525:       for (q = 0; q < numAdj; ++q) {
526:         const PetscInt  padj = tmpAdj[q];
527:         PetscInt        ndof, ncdof, ngoff, nd;
528:         const PetscInt *ncind;

530:         /* Adjacent points may not be in the section chart */
531:         if ((padj < pStart) || (padj >= pEnd)) continue;
532:         PetscSectionGetDof(section, padj, &ndof);
533:         PetscSectionGetConstraintDof(section, padj, &ncdof);
534:         PetscSectionGetConstraintIndices(section, padj, &ncind);
535:         PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
536:         for (nd = 0; nd < ndof - ncdof; ++nd, ++i) cols[aoff + i] = ngoff < 0 ? -(ngoff + 1) + nd : ngoff + nd;
537:       }
538:       for (q = 0; q < anDof; q++, i++) cols[aoff + i] = anchorAdj[anOff + q];
540:     }
541:   }
542:   PetscSectionDestroy(&anchorSectionAdj);
543:   PetscSectionDestroy(&leafSectionAdj);
544:   PetscSectionDestroy(&rootSectionAdj);
545:   PetscFree(anchorAdj);
546:   PetscFree(rootAdj);
547:   PetscFree(tmpAdj);
548:   /* Debugging */
549:   if (debug) {
550:     IS tmp;
551:     PetscPrintf(comm, "Column indices\n");
552:     ISCreateGeneral(comm, numCols, cols, PETSC_USE_POINTER, &tmp);
553:     ISView(tmp, NULL);
554:     ISDestroy(&tmp);
555:   }

557:   *sA     = sectionAdj;
558:   *colIdx = cols;
559:   return 0;
560: }

562: static PetscErrorCode DMPlexUpdateAllocation_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[])
563: {
564:   PetscSection section;
565:   PetscInt     rStart, rEnd, r, pStart, pEnd, p;

567:   /* This loop needs to change to a loop over points, then field dofs, which means we need to look both sections */
568:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
570:   if (f >= 0 && bs == 1) {
571:     DMGetLocalSection(dm, &section);
572:     PetscSectionGetChart(section, &pStart, &pEnd);
573:     for (p = pStart; p < pEnd; ++p) {
574:       PetscInt rS, rE;

576:       DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE);
577:       for (r = rS; r < rE; ++r) {
578:         PetscInt numCols, cStart, c;

580:         PetscSectionGetDof(sectionAdj, r, &numCols);
581:         PetscSectionGetOffset(sectionAdj, r, &cStart);
582:         for (c = cStart; c < cStart + numCols; ++c) {
583:           if ((cols[c] >= rStart) && (cols[c] < rEnd)) {
584:             ++dnz[r - rStart];
585:             if (cols[c] >= r) ++dnzu[r - rStart];
586:           } else {
587:             ++onz[r - rStart];
588:             if (cols[c] >= r) ++onzu[r - rStart];
589:           }
590:         }
591:       }
592:     }
593:   } else {
594:     /* Only loop over blocks of rows */
595:     for (r = rStart / bs; r < rEnd / bs; ++r) {
596:       const PetscInt row = r * bs;
597:       PetscInt       numCols, cStart, c;

599:       PetscSectionGetDof(sectionAdj, row, &numCols);
600:       PetscSectionGetOffset(sectionAdj, row, &cStart);
601:       for (c = cStart; c < cStart + numCols; ++c) {
602:         if ((cols[c] >= rStart) && (cols[c] < rEnd)) {
603:           ++dnz[r - rStart / bs];
604:           if (cols[c] >= row) ++dnzu[r - rStart / bs];
605:         } else {
606:           ++onz[r - rStart / bs];
607:           if (cols[c] >= row) ++onzu[r - rStart / bs];
608:         }
609:       }
610:     }
611:     for (r = 0; r < (rEnd - rStart) / bs; ++r) {
612:       dnz[r] /= bs;
613:       onz[r] /= bs;
614:       dnzu[r] /= bs;
615:       onzu[r] /= bs;
616:     }
617:   }
618:   return 0;
619: }

621: static PetscErrorCode DMPlexFillMatrix_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], Mat A)
622: {
623:   PetscSection section;
624:   PetscScalar *values;
625:   PetscInt     rStart, rEnd, r, pStart, pEnd, p, len, maxRowLen = 0;

627:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
628:   for (r = rStart; r < rEnd; ++r) {
629:     PetscSectionGetDof(sectionAdj, r, &len);
630:     maxRowLen = PetscMax(maxRowLen, len);
631:   }
632:   PetscCalloc1(maxRowLen, &values);
633:   if (f >= 0 && bs == 1) {
634:     DMGetLocalSection(dm, &section);
635:     PetscSectionGetChart(section, &pStart, &pEnd);
636:     for (p = pStart; p < pEnd; ++p) {
637:       PetscInt rS, rE;

639:       DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE);
640:       for (r = rS; r < rE; ++r) {
641:         PetscInt numCols, cStart;

643:         PetscSectionGetDof(sectionAdj, r, &numCols);
644:         PetscSectionGetOffset(sectionAdj, r, &cStart);
645:         MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);
646:       }
647:     }
648:   } else {
649:     for (r = rStart; r < rEnd; ++r) {
650:       PetscInt numCols, cStart;

652:       PetscSectionGetDof(sectionAdj, r, &numCols);
653:       PetscSectionGetOffset(sectionAdj, r, &cStart);
654:       MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);
655:     }
656:   }
657:   PetscFree(values);
658:   return 0;
659: }

661: /*@C
662:   DMPlexPreallocateOperator - Calculate the matrix nonzero pattern based upon the information in the `DM`,
663:   the `PetscDS` it contains, and the default `PetscSection`.

665:   Collective

667:   Input Parameters:
668: + dm   - The `DMPLEX`
669: . bs   - The matrix blocksize
670: . dnz  - An array to hold the number of nonzeros in the diagonal block
671: . onz  - An array to hold the number of nonzeros in the off-diagonal block
672: . dnzu - An array to hold the number of nonzeros in the upper triangle of the diagonal block
673: . onzu - An array to hold the number of nonzeros in the upper triangle of the off-diagonal block
674: - fillMatrix - If `PETSC_TRUE`, fill the matrix with zeros

676:   Output Parameter:
677: . A - The preallocated matrix

679:   Level: advanced

681: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMCreateMatrix()`
682: @*/
683: PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
684: {
685:   MPI_Comm     comm;
686:   PetscDS      prob;
687:   MatType      mtype;
688:   PetscSF      sf, sfDof;
689:   PetscSection section;
690:   PetscInt    *remoteOffsets;
691:   PetscSection sectionAdj[4] = {NULL, NULL, NULL, NULL};
692:   PetscInt    *cols[4]       = {NULL, NULL, NULL, NULL};
693:   PetscBool    useCone, useClosure;
694:   PetscInt     Nf, f, idx, locRows;
695:   PetscLayout  rLayout;
696:   PetscBool    isSymBlock, isSymSeqBlock, isSymMPIBlock, debug = PETSC_FALSE;
697:   PetscMPIInt  size;

705:   DMGetDS(dm, &prob);
706:   DMGetPointSF(dm, &sf);
707:   DMGetLocalSection(dm, &section);
708:   PetscOptionsGetBool(NULL, NULL, "-dm_view_preallocation", &debug, NULL);
709:   PetscObjectGetComm((PetscObject)dm, &comm);
710:   MPI_Comm_size(comm, &size);
711:   PetscLogEventBegin(DMPLEX_Preallocate, dm, 0, 0, 0);
712:   /* Create dof SF based on point SF */
713:   if (debug) {
714:     PetscSection section, sectionGlobal;
715:     PetscSF      sf;

717:     DMGetPointSF(dm, &sf);
718:     DMGetLocalSection(dm, &section);
719:     DMGetGlobalSection(dm, &sectionGlobal);
720:     PetscPrintf(comm, "Input Section for Preallocation:\n");
721:     PetscSectionView(section, NULL);
722:     PetscPrintf(comm, "Input Global Section for Preallocation:\n");
723:     PetscSectionView(sectionGlobal, NULL);
724:     if (size > 1) {
725:       PetscPrintf(comm, "Input SF for Preallocation:\n");
726:       PetscSFView(sf, NULL);
727:     }
728:   }
729:   PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets);
730:   PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof);
731:   PetscFree(remoteOffsets);
732:   if (debug && size > 1) {
733:     PetscPrintf(comm, "Dof SF for Preallocation:\n");
734:     PetscSFView(sfDof, NULL);
735:   }
736:   /* Create allocation vectors from adjacency graph */
737:   MatGetLocalSize(A, &locRows, NULL);
738:   PetscLayoutCreate(comm, &rLayout);
739:   PetscLayoutSetLocalSize(rLayout, locRows);
740:   PetscLayoutSetBlockSize(rLayout, 1);
741:   PetscLayoutSetUp(rLayout);
742:   /* There are 4 types of adjacency */
743:   PetscSectionGetNumFields(section, &Nf);
744:   if (Nf < 1 || bs > 1) {
745:     DMGetBasicAdjacency(dm, &useCone, &useClosure);
746:     idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
747:     DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, &sectionAdj[idx], &cols[idx]);
748:     DMPlexUpdateAllocation_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu);
749:   } else {
750:     for (f = 0; f < Nf; ++f) {
751:       DMGetAdjacency(dm, f, &useCone, &useClosure);
752:       idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
753:       if (!sectionAdj[idx]) DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, &sectionAdj[idx], &cols[idx]);
754:       DMPlexUpdateAllocation_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu);
755:     }
756:   }
757:   PetscSFDestroy(&sfDof);
758:   /* Set matrix pattern */
759:   MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu);
760:   MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
761:   /* Check for symmetric storage */
762:   MatGetType(A, &mtype);
763:   PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
764:   PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
765:   PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
766:   if (isSymBlock || isSymSeqBlock || isSymMPIBlock) MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);
767:   /* Fill matrix with zeros */
768:   if (fillMatrix) {
769:     if (Nf < 1 || bs > 1) {
770:       DMGetBasicAdjacency(dm, &useCone, &useClosure);
771:       idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
772:       DMPlexFillMatrix_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], A);
773:     } else {
774:       for (f = 0; f < Nf; ++f) {
775:         DMGetAdjacency(dm, f, &useCone, &useClosure);
776:         idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
777:         DMPlexFillMatrix_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], A);
778:       }
779:     }
780:     MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
781:     MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
782:   }
783:   PetscLayoutDestroy(&rLayout);
784:   for (idx = 0; idx < 4; ++idx) {
785:     PetscSectionDestroy(&sectionAdj[idx]);
786:     PetscFree(cols[idx]);
787:   }
788:   PetscLogEventEnd(DMPLEX_Preallocate, dm, 0, 0, 0);
789:   return 0;
790: }

792: #if 0
793: PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
794: {
795:   PetscInt       *tmpClosure,*tmpAdj,*visits;
796:   PetscInt        c,cStart,cEnd,pStart,pEnd;

798:   DMGetDimension(dm, &dim);
799:   DMPlexGetDepth(dm, &depth);
800:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);

802:   maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1));

804:   PetscSectionGetChart(section, &pStart, &pEnd);
805:   npoints = pEnd - pStart;

807:   PetscMalloc3(maxClosureSize,&tmpClosure,npoints,&lvisits,npoints,&visits);
808:   PetscArrayzero(lvisits,pEnd-pStart);
809:   PetscArrayzero(visits,pEnd-pStart);
810:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
811:   for (c=cStart; c<cEnd; c++) {
812:     PetscInt *support = tmpClosure;
813:     DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support);
814:     for (p=0; p<supportSize; p++) lvisits[support[p]]++;
815:   }
816:   PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM);
817:   PetscSFReduceEnd  (sf,MPIU_INT,lvisits,visits,MPI_SUM);
818:   PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits,MPI_REPLACE);
819:   PetscSFBcastEnd  (sf,MPIU_INT,visits,lvisits);

821:   PetscSFGetRootRanks();

823:   PetscMalloc2(maxClosureSize*maxClosureSize,&cellmat,npoints,&owner);
824:   for (c=cStart; c<cEnd; c++) {
825:     PetscArrayzero(cellmat,maxClosureSize*maxClosureSize);
826:     /*
827:      Depth-first walk of transitive closure.
828:      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.
829:      This contribution is added to dnz if owning ranks of p and q match, to onz otherwise.
830:      */
831:   }

833:   PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM);
834:   PetscSFReduceEnd  (sf,MPIU_INT,lonz,onz,MPI_SUM);
835:   return 0;
836: }
837: #endif