Actual source code: plexpreallocate.c

petsc-master 2019-10-20
Report Typos and Errors
  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;

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

 21:   DMPlexGetAnchors(dm,&aSec,&aIS);
 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;

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

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

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

 47:       PetscSectionGetDof(aSec, p, &dof);
 48:       PetscSectionGetOffset(aSec, p, &off);

 50:       for (q = 0; q < dof; q++) {
 51:         PetscInt iOff;

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

 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;

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

 88:         q = inverse[iOff + i];
 89:         DMPlexGetAdjacency_Internal(dm,q,useCone,useClosure,PETSC_TRUE,&numAdjQ,&tmpAdjQ);
 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:           PetscSectionGetDof(section,qAdj,&qAdjDof);
 98:           PetscSectionGetConstraintDof(section,qAdj,&qAdjCDof);
 99:           iNew += qAdjDof - qAdjCDof;
100:         }
101:         PetscSectionAddDof(adjSec,p,iNew);
102:       }
103:     }

105:     PetscSectionSetUp(adjSec);
106:     PetscSectionGetStorageSize(adjSec,&adjSize);
107:     PetscMalloc1(adjSize,&adj);

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

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

122:         q = inverse[iOff + i];
123:         DMPlexGetAdjacency_Internal(dm,q,useCone,useClosure,PETSC_TRUE,&numAdjQ,&tmpAdjQ);
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:           PetscSectionGetDof(section,qAdj,&qAdjDof);
132:           PetscSectionGetConstraintDof(section,qAdj,&qAdjCDof);
133:           PetscSectionGetOffset(sectionGlobal,qAdj,&qAdjOff);
134:           for (nd = 0; nd < qAdjDof-qAdjCDof; ++nd) {
135:             adj[aOff++] = (qAdjOff < 0 ? -(qAdjOff+1) : qAdjOff) + nd;
136:           }
137:         }
138:       }
139:       PetscSortRemoveDupsInt(&aDof,&adj[aOffOrig]);
140:       PetscSectionSetDof(adjSec,p,aDof);
141:     }
142:     *anchorAdj = adj;

144:     /* clean up */
145:     PetscSectionDestroy(&inverseSec);
146:     PetscFree(inverse);
147:     PetscFree(tmpAdjP);
148:     PetscFree(tmpAdjQ);
149:   }
150:   else {
151:     *anchorAdj = NULL;
152:     PetscSectionSetUp(adjSec);
153:   }
154:   *anchorSectionAdj = adjSec;
155:   return(0);
156: }

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

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

225:     if ((p < pStart) || (p >= pEnd)) continue;
226:     PetscSectionGetDof(section, p, &dof);
227:     PetscSectionGetOffset(section, p, &off);
228:     DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
229:     for (q = 0; q < numAdj; ++q) {
230:       const PetscInt padj = tmpAdj[q];
231:       PetscInt ndof, ncdof;

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

265:     PetscSectionGetDof(section, p, &dof);
266:     PetscSectionGetOffset(section, p, &off);
267:     if (!dof) continue;
268:     PetscSectionGetDof(rootSectionAdj, off, &adof);
269:     if (adof <= 0) continue;
270:     DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
271:     for (q = 0; q < numAdj; ++q) {
272:       const PetscInt padj = tmpAdj[q];
273:       PetscInt ndof, ncdof;

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

310:     if ((p < pStart) || (p >= pEnd)) continue;
311:     PetscSectionGetDof(section, p, &dof);
312:     PetscSectionGetOffset(section, p, &off);
313:     DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
314:     PetscSectionGetDof(anchorSectionAdj, p, &anDof);
315:     PetscSectionGetOffset(anchorSectionAdj, p, &anOff);
316:     for (d = off; d < off+dof; ++d) {
317:       PetscInt aoff, i = 0;

319:       PetscSectionGetOffset(leafSectionAdj, d, &aoff);
320:       for (q = 0; q < numAdj; ++q) {
321:         const PetscInt padj = tmpAdj[q];
322:         PetscInt ndof, ncdof, ngoff, nd;

324:         if ((padj < pStart) || (padj >= pEnd)) continue;
325:         PetscSectionGetDof(section, padj, &ndof);
326:         PetscSectionGetConstraintDof(section, padj, &ncdof);
327:         PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
328:         for (nd = 0; nd < ndof-ncdof; ++nd) {
329:           adj[aoff+i] = (ngoff < 0 ? -(ngoff+1) : ngoff) + nd;
330:           ++i;
331:         }
332:       }
333:       for (q = 0; q < anDof; q++) {
334:         adj[aoff+i] = anchorAdj[anOff+q];
335:         ++i;
336:       }
337:     }
338:   }
339:   /* Debugging */
340:   if (debug) {
341:     IS tmp;
342:     PetscPrintf(comm, "Leaf adjacency indices\n");
343:     ISCreateGeneral(comm, adjSize, adj, PETSC_USE_POINTER, &tmp);
344:     ISView(tmp, NULL);
345:     ISDestroy(&tmp);
346:   }
347:   /* Gather adjacent indices to root */
348:   PetscSectionGetStorageSize(rootSectionAdj, &adjSize);
349:   PetscMalloc1(adjSize, &rootAdj);
350:   for (r = 0; r < adjSize; ++r) rootAdj[r] = -1;
351:   if (doComm) {
352:     const PetscInt *indegree;
353:     PetscInt       *remoteadj, radjsize = 0;

355:     PetscSFComputeDegreeBegin(sfAdj, &indegree);
356:     PetscSFComputeDegreeEnd(sfAdj, &indegree);
357:     for (p = 0; p < adjSize; ++p) radjsize += indegree[p];
358:     PetscMalloc1(radjsize, &remoteadj);
359:     PetscSFGatherBegin(sfAdj, MPIU_INT, adj, remoteadj);
360:     PetscSFGatherEnd(sfAdj, MPIU_INT, adj, remoteadj);
361:     for (p = 0, l = 0, r = 0; p < adjSize; ++p, l = PetscMax(p, l + indegree[p-1])) {
362:       PetscInt s;
363:       for (s = 0; s < indegree[p]; ++s, ++r) rootAdj[l+s] = remoteadj[r];
364:     }
365:     if (r != radjsize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %d != %d", r, radjsize);
366:     if (l != adjSize)  SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %d != %d", l, adjSize);
367:     PetscFree(remoteadj);
368:   }
369:   PetscSFDestroy(&sfAdj);
370:   PetscFree(adj);
371:   /* Debugging */
372:   if (debug) {
373:     IS tmp;
374:     PetscPrintf(comm, "Root adjacency indices after gather\n");
375:     ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
376:     ISView(tmp, NULL);
377:     ISDestroy(&tmp);
378:   }
379:   /* Add in local adjacency indices for owned dofs on interface (roots) */
380:   for (p = pStart; p < pEnd; ++p) {
381:     PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof, anOff;

383:     PetscSectionGetDof(section, p, &dof);
384:     PetscSectionGetOffset(section, p, &off);
385:     if (!dof) continue;
386:     PetscSectionGetDof(rootSectionAdj, off, &adof);
387:     if (adof <= 0) continue;
388:     DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
389:     PetscSectionGetDof(anchorSectionAdj, p, &anDof);
390:     PetscSectionGetOffset(anchorSectionAdj, p, &anOff);
391:     for (d = off; d < off+dof; ++d) {
392:       PetscInt adof, aoff, i;

394:       PetscSectionGetDof(rootSectionAdj, d, &adof);
395:       PetscSectionGetOffset(rootSectionAdj, d, &aoff);
396:       i    = adof-1;
397:       for (q = 0; q < anDof; q++) {
398:         rootAdj[aoff+i] = anchorAdj[anOff+q];
399:         --i;
400:       }
401:       for (q = 0; q < numAdj; ++q) {
402:         const PetscInt padj = tmpAdj[q];
403:         PetscInt ndof, ncdof, ngoff, nd;

405:         if ((padj < pStart) || (padj >= pEnd)) continue;
406:         PetscSectionGetDof(section, padj, &ndof);
407:         PetscSectionGetConstraintDof(section, padj, &ncdof);
408:         PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
409:         for (nd = 0; nd < ndof-ncdof; ++nd) {
410:           rootAdj[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
411:           --i;
412:         }
413:       }
414:     }
415:   }
416:   /* Debugging */
417:   if (debug) {
418:     IS tmp;
419:     PetscPrintf(comm, "Root adjacency indices\n");
420:     ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
421:     ISView(tmp, NULL);
422:     ISDestroy(&tmp);
423:   }
424:   /* Compress indices */
425:   PetscSectionSetUp(rootSectionAdj);
426:   for (p = pStart; p < pEnd; ++p) {
427:     PetscInt dof, cdof, off, d;
428:     PetscInt adof, aoff;

430:     PetscSectionGetDof(section, p, &dof);
431:     PetscSectionGetConstraintDof(section, p, &cdof);
432:     PetscSectionGetOffset(section, p, &off);
433:     if (!dof) continue;
434:     PetscSectionGetDof(rootSectionAdj, off, &adof);
435:     if (adof <= 0) continue;
436:     for (d = off; d < off+dof-cdof; ++d) {
437:       PetscSectionGetDof(rootSectionAdj, d, &adof);
438:       PetscSectionGetOffset(rootSectionAdj, d, &aoff);
439:       PetscSortRemoveDupsInt(&adof, &rootAdj[aoff]);
440:       PetscSectionSetDof(rootSectionAdj, d, adof);
441:     }
442:   }
443:   /* Debugging */
444:   if (debug) {
445:     IS tmp;
446:     PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after compression:\n");
447:     PetscSectionView(rootSectionAdj, NULL);
448:     PetscPrintf(comm, "Root adjacency indices after compression\n");
449:     ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
450:     ISView(tmp, NULL);
451:     ISDestroy(&tmp);
452:   }
453:   /* Build adjacency section: Maps global indices to sets of adjacent global indices */
454:   PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd);
455:   PetscSectionCreate(comm, &sectionAdj);
456:   PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd);
457:   for (p = pStart; p < pEnd; ++p) {
458:     PetscInt  numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof;
459:     PetscBool found  = PETSC_TRUE;

461:     PetscSectionGetDof(section, p, &dof);
462:     PetscSectionGetConstraintDof(section, p, &cdof);
463:     PetscSectionGetOffset(section, p, &off);
464:     PetscSectionGetOffset(sectionGlobal, p, &goff);
465:     for (d = 0; d < dof-cdof; ++d) {
466:       PetscInt ldof, rdof;

468:       PetscSectionGetDof(leafSectionAdj, off+d, &ldof);
469:       PetscSectionGetDof(rootSectionAdj, off+d, &rdof);
470:       if (ldof > 0) {
471:         /* We do not own this point */
472:       } else if (rdof > 0) {
473:         PetscSectionSetDof(sectionAdj, goff+d, rdof);
474:       } else {
475:         found = PETSC_FALSE;
476:       }
477:     }
478:     if (found) continue;
479:     PetscSectionGetDof(section, p, &dof);
480:     PetscSectionGetOffset(sectionGlobal, p, &goff);
481:     DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
482:     for (q = 0; q < numAdj; ++q) {
483:       const PetscInt padj = tmpAdj[q];
484:       PetscInt ndof, ncdof, noff;

486:       if ((padj < pStart) || (padj >= pEnd)) continue;
487:       PetscSectionGetDof(section, padj, &ndof);
488:       PetscSectionGetConstraintDof(section, padj, &ncdof);
489:       PetscSectionGetOffset(section, padj, &noff);
490:       for (d = goff; d < goff+dof-cdof; ++d) {
491:         PetscSectionAddDof(sectionAdj, d, ndof-ncdof);
492:       }
493:     }
494:     PetscSectionGetDof(anchorSectionAdj, p, &anDof);
495:     if (anDof) {
496:       for (d = goff; d < goff+dof-cdof; ++d) {
497:         PetscSectionAddDof(sectionAdj, d, anDof);
498:       }
499:     }
500:   }
501:   PetscSectionSetUp(sectionAdj);
502:   if (debug) {
503:     PetscPrintf(comm, "Adjacency Section for Preallocation:\n");
504:     PetscSectionView(sectionAdj, NULL);
505:   }
506:   /* Get adjacent indices */
507:   PetscSectionGetStorageSize(sectionAdj, &numCols);
508:   PetscMalloc1(numCols, &cols);
509:   for (p = pStart; p < pEnd; ++p) {
510:     PetscInt  numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof, anOff;
511:     PetscBool found  = PETSC_TRUE;

513:     PetscSectionGetDof(section, p, &dof);
514:     PetscSectionGetConstraintDof(section, p, &cdof);
515:     PetscSectionGetOffset(section, p, &off);
516:     PetscSectionGetOffset(sectionGlobal, p, &goff);
517:     for (d = 0; d < dof-cdof; ++d) {
518:       PetscInt ldof, rdof;

520:       PetscSectionGetDof(leafSectionAdj, off+d, &ldof);
521:       PetscSectionGetDof(rootSectionAdj, off+d, &rdof);
522:       if (ldof > 0) {
523:         /* We do not own this point */
524:       } else if (rdof > 0) {
525:         PetscInt aoff, roff;

527:         PetscSectionGetOffset(sectionAdj, goff+d, &aoff);
528:         PetscSectionGetOffset(rootSectionAdj, off+d, &roff);
529:         PetscArraycpy(&cols[aoff], &rootAdj[roff], rdof);
530:       } else {
531:         found = PETSC_FALSE;
532:       }
533:     }
534:     if (found) continue;
535:     DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj);
536:     PetscSectionGetDof(anchorSectionAdj, p, &anDof);
537:     PetscSectionGetOffset(anchorSectionAdj, p, &anOff);
538:     for (d = goff; d < goff+dof-cdof; ++d) {
539:       PetscInt adof, aoff, i = 0;

541:       PetscSectionGetDof(sectionAdj, d, &adof);
542:       PetscSectionGetOffset(sectionAdj, d, &aoff);
543:       for (q = 0; q < numAdj; ++q) {
544:         const PetscInt  padj = tmpAdj[q];
545:         PetscInt        ndof, ncdof, ngoff, nd;
546:         const PetscInt *ncind;

548:         /* Adjacent points may not be in the section chart */
549:         if ((padj < pStart) || (padj >= pEnd)) continue;
550:         PetscSectionGetDof(section, padj, &ndof);
551:         PetscSectionGetConstraintDof(section, padj, &ncdof);
552:         PetscSectionGetConstraintIndices(section, padj, &ncind);
553:         PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
554:         for (nd = 0; nd < ndof-ncdof; ++nd, ++i) {
555:           cols[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
556:         }
557:       }
558:       for (q = 0; q < anDof; q++, i++) {
559:         cols[aoff+i] = anchorAdj[anOff + q];
560:       }
561:       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);
562:     }
563:   }
564:   PetscSectionDestroy(&anchorSectionAdj);
565:   PetscSectionDestroy(&leafSectionAdj);
566:   PetscSectionDestroy(&rootSectionAdj);
567:   PetscFree(anchorAdj);
568:   PetscFree(rootAdj);
569:   PetscFree(tmpAdj);
570:   /* Debugging */
571:   if (debug) {
572:     IS tmp;
573:     PetscPrintf(comm, "Column indices\n");
574:     ISCreateGeneral(comm, numCols, cols, PETSC_USE_POINTER, &tmp);
575:     ISView(tmp, NULL);
576:     ISDestroy(&tmp);
577:   }

579:   *sA     = sectionAdj;
580:   *colIdx = cols;
581:   return(0);
582: }

584: 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[])
585: {
586:   PetscSection   section;
587:   PetscInt       rStart, rEnd, r, pStart, pEnd, p;

591:   /* This loop needs to change to a loop over points, then field dofs, which means we need to look both sections */
592:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
593:   if (rStart%bs || rEnd%bs) SETERRQ3(PetscObjectComm((PetscObject) rLayout), PETSC_ERR_ARG_WRONG, "Invalid layout [%d, %d) for matrix, must be divisible by block size %d", rStart, rEnd, bs);
594:   if (f >= 0 && bs == 1) {
595:     DMGetLocalSection(dm, &section);
596:     PetscSectionGetChart(section, &pStart, &pEnd);
597:     for (p = pStart; p < pEnd; ++p) {
598:       PetscInt rS, rE;

600:       DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE);
601:       for (r = rS; r < rE; ++r) {
602:         PetscInt numCols, cStart, c;

604:         PetscSectionGetDof(sectionAdj, r, &numCols);
605:         PetscSectionGetOffset(sectionAdj, r, &cStart);
606:         for (c = cStart; c < cStart+numCols; ++c) {
607:           if ((cols[c] >= rStart) && (cols[c] < rEnd)) {
608:             ++dnz[r-rStart];
609:             if (cols[c] >= r) ++dnzu[r-rStart];
610:           } else {
611:             ++onz[r-rStart];
612:             if (cols[c] >= r) ++onzu[r-rStart];
613:           }
614:         }
615:       }
616:     }
617:   } else {
618:     /* Only loop over blocks of rows */
619:     for (r = rStart/bs; r < rEnd/bs; ++r) {
620:       const PetscInt row = r*bs;
621:       PetscInt       numCols, cStart, c;

623:       PetscSectionGetDof(sectionAdj, row, &numCols);
624:       PetscSectionGetOffset(sectionAdj, row, &cStart);
625:       for (c = cStart; c < cStart+numCols; ++c) {
626:         if ((cols[c] >= rStart) && (cols[c] < rEnd)) {
627:           ++dnz[r-rStart/bs];
628:           if (cols[c] >= row) ++dnzu[r-rStart/bs];
629:         } else {
630:           ++onz[r-rStart/bs];
631:           if (cols[c] >= row) ++onzu[r-rStart/bs];
632:         }
633:       }
634:     }
635:     for (r = 0; r < (rEnd - rStart)/bs; ++r) {
636:       dnz[r]  /= bs;
637:       onz[r]  /= bs;
638:       dnzu[r] /= bs;
639:       onzu[r] /= bs;
640:     }
641:   }
642:   return(0);
643: }

645: static PetscErrorCode DMPlexFillMatrix_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], Mat A)
646: {
647:   PetscSection   section;
648:   PetscScalar   *values;
649:   PetscInt       rStart, rEnd, r, pStart, pEnd, p, len, maxRowLen = 0;

653:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
654:   for (r = rStart; r < rEnd; ++r) {
655:     PetscSectionGetDof(sectionAdj, r, &len);
656:     maxRowLen = PetscMax(maxRowLen, len);
657:   }
658:   PetscCalloc1(maxRowLen, &values);
659:   if (f >=0 && bs == 1) {
660:     DMGetLocalSection(dm, &section);
661:     PetscSectionGetChart(section, &pStart, &pEnd);
662:     for (p = pStart; p < pEnd; ++p) {
663:       PetscInt rS, rE;

665:       DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE);
666:       for (r = rS; r < rE; ++r) {
667:         PetscInt numCols, cStart;

669:         PetscSectionGetDof(sectionAdj, r, &numCols);
670:         PetscSectionGetOffset(sectionAdj, r, &cStart);
671:         MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);
672:       }
673:     }
674:   } else {
675:     for (r = rStart; r < rEnd; ++r) {
676:       PetscInt numCols, cStart;

678:       PetscSectionGetDof(sectionAdj, r, &numCols);
679:       PetscSectionGetOffset(sectionAdj, r, &cStart);
680:       MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);
681:     }
682:   }
683:   PetscFree(values);
684:   return(0);
685: }

687: /*@C
688:   DMPlexPreallocateOperator - Calculate the matrix nonzero pattern based upon the information in the DM,
689:   the PetscDS it contains, and the default PetscSection.

691:   Collective

693:   Input Arguments:
694: + dm   - The DMPlex
695: . bs   - The matrix blocksize
696: . dnz  - An array to hold the number of nonzeros in the diagonal block
697: . onz  - An array to hold the number of nonzeros in the off-diagonal block
698: . dnzu - An array to hold the number of nonzeros in the upper triangle of the diagonal block
699: . onzu - An array to hold the number of nonzeros in the upper triangle of the off-diagonal block
700: - fillMatrix - If PETSC_TRUE, fill the matrix with zeros

702:   Ouput Argument:
703: . A - The preallocated matrix

705:   Level: advanced

707: .seealso: DMCreateMatrix()
708: @*/
709: PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
710: {
711:   MPI_Comm       comm;
712:   PetscDS        prob;
713:   MatType        mtype;
714:   PetscSF        sf, sfDof;
715:   PetscSection   section;
716:   PetscInt      *remoteOffsets;
717:   PetscSection   sectionAdj[4] = {NULL, NULL, NULL, NULL};
718:   PetscInt      *cols[4]       = {NULL, NULL, NULL, NULL};
719:   PetscBool      useCone, useClosure;
720:   PetscInt       Nf, f, idx, locRows;
721:   PetscLayout    rLayout;
722:   PetscBool      isSymBlock, isSymSeqBlock, isSymMPIBlock, debug = PETSC_FALSE;
723:   PetscMPIInt    size;

731:   DMGetDS(dm, &prob);
732:   DMGetPointSF(dm, &sf);
733:   DMGetLocalSection(dm, &section);
734:   PetscOptionsGetBool(NULL,NULL, "-dm_view_preallocation", &debug, NULL);
735:   PetscObjectGetComm((PetscObject) dm, &comm);
736:   MPI_Comm_size(comm, &size);
737:   PetscLogEventBegin(DMPLEX_Preallocate,dm,0,0,0);
738:   /* Create dof SF based on point SF */
739:   if (debug) {
740:     PetscSection section, sectionGlobal;
741:     PetscSF      sf;

743:     DMGetPointSF(dm, &sf);
744:     DMGetLocalSection(dm, &section);
745:     DMGetGlobalSection(dm, &sectionGlobal);
746:     PetscPrintf(comm, "Input Section for Preallocation:\n");
747:     PetscSectionView(section, NULL);
748:     PetscPrintf(comm, "Input Global Section for Preallocation:\n");
749:     PetscSectionView(sectionGlobal, NULL);
750:     if (size > 1) {
751:       PetscPrintf(comm, "Input SF for Preallocation:\n");
752:       PetscSFView(sf, NULL);
753:     }
754:   }
755:   PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets);
756:   PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof);
757:   PetscFree(remoteOffsets);
758:   if (debug && size > 1) {
759:     PetscPrintf(comm, "Dof SF for Preallocation:\n");
760:     PetscSFView(sfDof, NULL);
761:   }
762:   /* Create allocation vectors from adjacency graph */
763:   MatGetLocalSize(A, &locRows, NULL);
764:   PetscLayoutCreate(comm, &rLayout);
765:   PetscLayoutSetLocalSize(rLayout, locRows);
766:   PetscLayoutSetBlockSize(rLayout, 1);
767:   PetscLayoutSetUp(rLayout);
768:   /* There are 4 types of adjacency */
769:   PetscSectionGetNumFields(section, &Nf);
770:   if (Nf < 1 || bs > 1) {
771:     DMGetBasicAdjacency(dm, &useCone, &useClosure);
772:     idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
773:     DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, &sectionAdj[idx], &cols[idx]);
774:     DMPlexUpdateAllocation_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu);
775:   } else {
776:     for (f = 0; f < Nf; ++f) {
777:       DMGetAdjacency(dm, f, &useCone, &useClosure);
778:       idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
779:       if (!sectionAdj[idx]) {DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, &sectionAdj[idx], &cols[idx]);}
780:       DMPlexUpdateAllocation_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu);
781:     }
782:   }
783:   PetscSFDestroy(&sfDof);
784:   /* Set matrix pattern */
785:   MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu);
786:   MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
787:   /* Check for symmetric storage */
788:   MatGetType(A, &mtype);
789:   PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
790:   PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
791:   PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
792:   if (isSymBlock || isSymSeqBlock || isSymMPIBlock) {MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);}
793:   /* Fill matrix with zeros */
794:   if (fillMatrix) {
795:     if (Nf < 1 || bs > 1) {
796:       DMGetBasicAdjacency(dm, &useCone, &useClosure);
797:       idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
798:       DMPlexFillMatrix_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], A);
799:     } else {
800:       for (f = 0; f < Nf; ++f) {
801:         DMGetAdjacency(dm, f, &useCone, &useClosure);
802:         idx  = (useCone ? 1 : 0) + (useClosure ? 2 : 0);
803:         DMPlexFillMatrix_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], A);
804:       }
805:     }
806:     MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
807:     MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
808:   }
809:   PetscLayoutDestroy(&rLayout);
810:   for (idx = 0; idx < 4; ++idx) {PetscSectionDestroy(&sectionAdj[idx]); PetscFree(cols[idx]);}
811:   PetscLogEventEnd(DMPLEX_Preallocate,dm,0,0,0);
812:   return(0);
813: }

815: #if 0
816: PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
817: {
818:   PetscInt       *tmpClosure,*tmpAdj,*visits;
819:   PetscInt        c,cStart,cEnd,pStart,pEnd;
820:   PetscErrorCode  ierr;

823:   DMGetDimension(dm, &dim);
824:   DMPlexGetDepth(dm, &depth);
825:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);

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

829:   PetscSectionGetChart(section, &pStart, &pEnd);
830:   npoints = pEnd - pStart;

832:   PetscMalloc3(maxClosureSize,&tmpClosure,npoints,&lvisits,npoints,&visits);
833:   PetscArrayzero(lvisits,pEnd-pStart);
834:   PetscArrayzero(visits,pEnd-pStart);
835:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
836:   for (c=cStart; c<cEnd; c++) {
837:     PetscInt *support = tmpClosure;
838:     DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support);
839:     for (p=0; p<supportSize; p++) lvisits[support[p]]++;
840:   }
841:   PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM);
842:   PetscSFReduceEnd  (sf,MPIU_INT,lvisits,visits,MPI_SUM);
843:   PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits);
844:   PetscSFBcastEnd  (sf,MPIU_INT,visits,lvisits);

846:   PetscSFGetRootRanks();


849:   PetscMalloc2(maxClosureSize*maxClosureSize,&cellmat,npoints,&owner);
850:   for (c=cStart; c<cEnd; c++) {
851:     PetscArrayzero(cellmat,maxClosureSize*maxClosureSize);
852:     /*
853:      Depth-first walk of transitive closure.
854:      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.
855:      This contribution is added to dnz if owning ranks of p and q match, to onz otherwise.
856:      */
857:   }

859:   PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM);
860:   PetscSFReduceEnd  (sf,MPIU_INT,lonz,onz,MPI_SUM);
861:   return(0);
862: }
863: #endif