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;

     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,MPI_REPLACE);
    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