Actual source code: plexpreallocate.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  1: #include <petsc/private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
  2: #include <petsc/private/isimpl.h>
  3: #include <petscsf.h>
  4: #include <petscds.h>

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

 17:   DMGetDefaultSection(dm, &section);
 18:   DMGetDefaultGlobalSection(dm, &sectionGlobal);
 19:   PetscSectionCreate(PetscObjectComm((PetscObject) section), &adjSec);
 20:   PetscSectionGetChart(section,&pStart,&pEnd);
 21:   PetscSectionSetChart(adjSec,pStart,pEnd);

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

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

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

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

 49:       PetscSectionGetDof(aSec, p, &dof);
 50:       PetscSectionGetOffset(aSec, p, &off);

 52:       for (q = 0; q < dof; q++) {
 53:         PetscInt iOff;

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

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

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

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

107:     PetscSectionSetUp(adjSec);
108:     PetscSectionGetStorageSize(adjSec,&adjSize);
109:     PetscMalloc1(adjSize,&adj);

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

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

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

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

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

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

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

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

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

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

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

323:       PetscSectionGetOffset(leafSectionAdj, d, &aoff);
324:       for (q = 0; q < numAdj; ++q) {
325:         const PetscInt padj = tmpAdj[q];
326:         PetscInt ndof, ncdof, ngoff, nd;

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

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

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

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

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

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

465:     PetscSectionGetDof(section, p, &dof);
466:     PetscSectionGetConstraintDof(section, p, &cdof);
467:     PetscSectionGetOffset(section, p, &off);
468:     PetscSectionGetOffset(sectionGlobal, p, &goff);
469:     for (d = 0; d < dof-cdof; ++d) {
470:       PetscInt ldof, rdof;

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

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

517:     PetscSectionGetDof(section, p, &dof);
518:     PetscSectionGetConstraintDof(section, p, &cdof);
519:     PetscSectionGetOffset(section, p, &off);
520:     PetscSectionGetOffset(sectionGlobal, p, &goff);
521:     for (d = 0; d < dof-cdof; ++d) {
522:       PetscInt ldof, rdof;

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

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

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

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

583:   *sA     = sectionAdj;
584:   *colIdx = cols;
585:   return(0);
586: }

590: PetscErrorCode DMPlexUpdateAllocation_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[])
591: {
592:   PetscSection   section;
593:   PetscInt       rStart, rEnd, r, pStart, pEnd, p;

597:   /* This loop needs to change to a loop over points, then field dofs, which means we need to look both sections */
598:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
599:   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);
600:   if (f >= 0 && bs == 1) {
601:     DMGetDefaultSection(dm, &section);
602:     PetscSectionGetChart(section, &pStart, &pEnd);
603:     for (p = pStart; p < pEnd; ++p) {
604:       PetscInt rS, rE;

606:       DMPlexGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE);
607:       for (r = rS; r < rE; ++r) {
608:         PetscInt numCols, cStart, c;

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

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

653: PetscErrorCode DMPlexFillMatrix_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], Mat A)
654: {
655:   PetscSection   section;
656:   PetscScalar   *values;
657:   PetscInt       rStart, rEnd, r, pStart, pEnd, p, len, maxRowLen = 0;

661:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
662:   for (r = rStart; r < rEnd; ++r) {
663:     PetscSectionGetDof(sectionAdj, r, &len);
664:     maxRowLen = PetscMax(maxRowLen, len);
665:   }
666:   PetscCalloc1(maxRowLen, &values);
667:   if (f >=0 && bs == 1) {
668:     DMGetDefaultSection(dm, &section);
669:     PetscSectionGetChart(section, &pStart, &pEnd);
670:     for (p = pStart; p < pEnd; ++p) {
671:       PetscInt rS, rE;

673:       DMPlexGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE);
674:       for (r = rS; r < rE; ++r) {
675:         PetscInt numCols, cStart;

677:         PetscSectionGetDof(sectionAdj, r, &numCols);
678:         PetscSectionGetOffset(sectionAdj, r, &cStart);
679:         MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);
680:       }
681:     }
682:   } else {
683:     for (r = rStart; r < rEnd; ++r) {
684:       PetscInt numCols, cStart;

686:       PetscSectionGetDof(sectionAdj, r, &numCols);
687:       PetscSectionGetOffset(sectionAdj, r, &cStart);
688:       MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);
689:     }
690:   }
691:   PetscFree(values);
692:   return(0);
693: }

697: PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
698: {
699:   MPI_Comm       comm;
700:   PetscDS        prob;
701:   MatType        mtype;
702:   PetscSF        sf, sfDof;
703:   PetscSection   section;
704:   PetscInt      *remoteOffsets;
705:   PetscSection   sectionAdj[4] = {NULL, NULL, NULL, NULL};
706:   PetscInt      *cols[4]       = {NULL, NULL, NULL, NULL};
707:   PetscBool      useCone, useClosure;
708:   PetscInt       Nf, f, idx, locRows;
709:   PetscLayout    rLayout;
710:   PetscBool      isSymBlock, isSymSeqBlock, isSymMPIBlock, debug = PETSC_FALSE;

718:   DMGetDS(dm, &prob);
719:   DMGetPointSF(dm, &sf);
720:   DMGetDefaultSection(dm, &section);
721:   PetscOptionsGetBool(NULL,NULL, "-dm_view_preallocation", &debug, NULL);
722:   PetscObjectGetComm((PetscObject) dm, &comm);
723:   PetscLogEventBegin(DMPLEX_Preallocate,dm,0,0,0);
724:   /* Create dof SF based on point SF */
725:   if (debug) {
726:     PetscSection section, sectionGlobal;
727:     PetscSF      sf;

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

801: #if 0
804: PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
805: {
806:   PetscInt       *tmpClosure,*tmpAdj,*visits;
807:   PetscInt        c,cStart,cEnd,pStart,pEnd;
808:   PetscErrorCode  ierr;

811:   DMGetDimension(dm, &dim);
812:   DMPlexGetDepth(dm, &depth);
813:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);

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

817:   PetscSectionGetChart(section, &pStart, &pEnd);
818:   npoints = pEnd - pStart;

820:   PetscMalloc3(maxClosureSize,&tmpClosure,npoints,&lvisits,npoints,&visits);
821:   PetscMemzero(lvisits,(pEnd-pStart)*sizeof(PetscInt));
822:   PetscMemzero(visits,(pEnd-pStart)*sizeof(PetscInt));
823:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
824:   for (c=cStart; c<cEnd; c++) {
825:     PetscInt *support = tmpClosure;
826:     DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support);
827:     for (p=0; p<supportSize; p++) lvisits[support[p]]++;
828:   }
829:   PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM);
830:   PetscSFReduceEnd  (sf,MPIU_INT,lvisits,visits,MPI_SUM);
831:   PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits);
832:   PetscSFBcastEnd  (sf,MPIU_INT,visits,lvisits);

834:   PetscSFGetRanks();


837:   PetscMalloc2(maxClosureSize*maxClosureSize,&cellmat,npoints,&owner);
838:   for (c=cStart; c<cEnd; c++) {
839:     PetscMemzero(cellmat,maxClosureSize*maxClosureSize*sizeof(PetscInt));
840:     /*
841:      Depth-first walk of transitive closure.
842:      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.
843:      This contribution is added to dnz if owning ranks of p and q match, to onz otherwise.
844:      */
845:   }

847:   PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM);
848:   PetscSFReduceEnd  (sf,MPIU_INT,lonz,onz,MPI_SUM);
849:   return(0);
850: }
851: #endif