Actual source code: sfutils.c

      1: #include <petsc/private/sfimpl.h>
      2: #include <petsc/private/sectionimpl.h>

      4: /*@C
      5:    PetscSFSetGraphLayout - Set a parallel star forest via global indices and a PetscLayout

      7:    Collective

      9:    Input Arguments:
     10: +  sf - star forest
     11: .  layout - PetscLayout defining the global space
     12: .  nleaves - number of leaf vertices on the current process, each of these references a root on any process
     13: .  ilocal - locations of leaves in leafdata buffers, pass NULL for contiguous storage
     14: .  localmode - copy mode for ilocal
     15: -  iremote - remote locations of root vertices for each leaf on the current process

     17:    Level: intermediate

     19:    Developers Note: Local indices which are the identity permutation in the range [0,nleaves) are discarded as they
     20:    encode contiguous storage. In such case, if localmode is PETSC_OWN_POINTER, the memory is deallocated as it is not
     21:    needed

     23: .seealso: PetscSFCreate(), PetscSFView(), PetscSFSetGraph(), PetscSFGetGraph()
     24: @*/
     25: PetscErrorCode PetscSFSetGraphLayout(PetscSF sf,PetscLayout layout,PetscInt nleaves,const PetscInt *ilocal,PetscCopyMode localmode,const PetscInt *iremote)
     26: {
     28:   const PetscInt *range;
     29:   PetscInt       i, nroots, ls = -1, ln = -1;
     30:   PetscMPIInt    lr = -1;
     31:   PetscSFNode    *remote;

     34:   PetscLayoutGetLocalSize(layout,&nroots);
     35:   PetscLayoutGetRanges(layout,&range);
     36:   PetscMalloc1(nleaves,&remote);
     37:   if (nleaves) { ls = iremote[0] + 1; }
     38:   for (i=0; i<nleaves; i++) {
     39:     const PetscInt idx = iremote[i] - ls;
     40:     if (idx < 0 || idx >= ln) { /* short-circuit the search */
     41:       PetscLayoutFindOwnerIndex(layout,iremote[i],&lr,&remote[i].index);
     42:       remote[i].rank = lr;
     43:       ls = range[lr];
     44:       ln = range[lr+1] - ls;
     45:     } else {
     46:       remote[i].rank  = lr;
     47:       remote[i].index = idx;
     48:     }
     49:   }
     50:   PetscSFSetGraph(sf,nroots,nleaves,ilocal,localmode,remote,PETSC_OWN_POINTER);
     51:   return(0);
     52: }

     54: /*@
     55:   PetscSFSetGraphSection - Sets the PetscSF graph encoding the parallel dof overlap based upon the PetscSections
     56:   describing the data layout.

     58:   Input Parameters:
     59: + sf - The SF
     60: . localSection - PetscSection describing the local data layout
     61: - globalSection - PetscSection describing the global data layout

     63:   Level: developer

     65: .seealso: PetscSFSetGraph(), PetscSFSetGraphLayout()
     66: @*/
     67: PetscErrorCode PetscSFSetGraphSection(PetscSF sf, PetscSection localSection, PetscSection globalSection)
     68: {
     69:   MPI_Comm       comm;
     70:   PetscLayout    layout;
     71:   const PetscInt *ranges;
     72:   PetscInt       *local;
     73:   PetscSFNode    *remote;
     74:   PetscInt       pStart, pEnd, p, nroots, nleaves = 0, l;
     75:   PetscMPIInt    size, rank;

     83:   PetscObjectGetComm((PetscObject)sf,&comm);
     84:   MPI_Comm_size(comm, &size);
     85:   MPI_Comm_rank(comm, &rank);
     86:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
     87:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
     88:   PetscLayoutCreate(comm, &layout);
     89:   PetscLayoutSetBlockSize(layout, 1);
     90:   PetscLayoutSetLocalSize(layout, nroots);
     91:   PetscLayoutSetUp(layout);
     92:   PetscLayoutGetRanges(layout, &ranges);
     93:   for (p = pStart; p < pEnd; ++p) {
     94:     PetscInt gdof, gcdof;

     96:     PetscSectionGetDof(globalSection, p, &gdof);
     97:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
     98:     if (gcdof > (gdof < 0 ? -(gdof+1) : gdof)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %D has %D constraints > %D dof", p, gcdof, (gdof < 0 ? -(gdof+1) : gdof));
     99:     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
    100:   }
    101:   PetscMalloc1(nleaves, &local);
    102:   PetscMalloc1(nleaves, &remote);
    103:   for (p = pStart, l = 0; p < pEnd; ++p) {
    104:     const PetscInt *cind;
    105:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

    107:     PetscSectionGetDof(localSection, p, &dof);
    108:     PetscSectionGetOffset(localSection, p, &off);
    109:     PetscSectionGetConstraintDof(localSection, p, &cdof);
    110:     PetscSectionGetConstraintIndices(localSection, p, &cind);
    111:     PetscSectionGetDof(globalSection, p, &gdof);
    112:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
    113:     PetscSectionGetOffset(globalSection, p, &goff);
    114:     if (!gdof) continue; /* Censored point */
    115:     gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
    116:     if (gsize != dof-cdof) {
    117:       if (gsize != dof) SETERRQ4(comm, PETSC_ERR_ARG_WRONG, "Global dof %D for point %D is neither the constrained size %D, nor the unconstrained %D", gsize, p, dof-cdof, dof);
    118:       cdof = 0; /* Ignore constraints */
    119:     }
    120:     for (d = 0, c = 0; d < dof; ++d) {
    121:       if ((c < cdof) && (cind[c] == d)) {++c; continue;}
    122:       local[l+d-c] = off+d;
    123:     }
    124:     if (gdof < 0) {
    125:       for (d = 0; d < gsize; ++d, ++l) {
    126:         PetscInt offset = -(goff+1) + d, r;

    128:         PetscFindInt(offset,size+1,ranges,&r);
    129:         if (r < 0) r = -(r+2);
    130:         if ((r < 0) || (r >= size)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %D mapped to invalid process %D (%D, %D)", p, r, gdof, goff);
    131:         remote[l].rank  = r;
    132:         remote[l].index = offset - ranges[r];
    133:       }
    134:     } else {
    135:       for (d = 0; d < gsize; ++d, ++l) {
    136:         remote[l].rank  = rank;
    137:         remote[l].index = goff+d - ranges[rank];
    138:       }
    139:     }
    140:   }
    141:   if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %D != nleaves %D", l, nleaves);
    142:   PetscLayoutDestroy(&layout);
    143:   PetscSFSetGraph(sf, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
    144:   return(0);
    145: }

    147: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
    148: {

    152:   if (!s->bc) {
    153:     PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
    154:     PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
    155:   }
    156:   return(0);
    157: }

    159: /*@C
    160:   PetscSFDistributeSection - Create a new PetscSection reorganized, moving from the root to the leaves of the SF

    162:   Collective on sf

    164:   Input Parameters:
    165: + sf - The SF
    166: - rootSection - Section defined on root space

    168:   Output Parameters:
    169: + remoteOffsets - root offsets in leaf storage, or NULL
    170: - leafSection - Section defined on the leaf space

    172:   Level: advanced

    174: .seealso: PetscSFCreate()
    175: @*/
    176: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
    177: {
    178:   PetscSF        embedSF;
    179:   const PetscInt *indices;
    180:   IS             selected;
    181:   PetscInt       numFields, nroots, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, f, c;
    182:   PetscBool      *sub, hasc;

    186:   PetscLogEventBegin(PETSCSF_DistSect,sf,0,0,0);
    187:   PetscSectionGetNumFields(rootSection, &numFields);
    188:   if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
    189:   PetscMalloc1(numFields+2, &sub);
    190:   sub[1] = rootSection->bc ? PETSC_TRUE : PETSC_FALSE;
    191:   for (f = 0; f < numFields; ++f) {
    192:     PetscSectionSym sym;
    193:     const char      *name   = NULL;
    194:     PetscInt        numComp = 0;

    196:     sub[2 + f] = rootSection->field[f]->bc ? PETSC_TRUE : PETSC_FALSE;
    197:     PetscSectionGetFieldComponents(rootSection, f, &numComp);
    198:     PetscSectionGetFieldName(rootSection, f, &name);
    199:     PetscSectionGetFieldSym(rootSection, f, &sym);
    200:     PetscSectionSetFieldComponents(leafSection, f, numComp);
    201:     PetscSectionSetFieldName(leafSection, f, name);
    202:     PetscSectionSetFieldSym(leafSection, f, sym);
    203:     for (c = 0; c < rootSection->numFieldComponents[f]; ++c) {
    204:       PetscSectionGetComponentName(rootSection, f, c, &name);
    205:       PetscSectionSetComponentName(leafSection, f, c, name);
    206:     }
    207:   }
    208:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
    209:   PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);
    210:   rpEnd = PetscMin(rpEnd,nroots);
    211:   rpEnd = PetscMax(rpStart,rpEnd);
    212:   /* see if we can avoid creating the embedded SF, since it can cost more than an allreduce */
    213:   sub[0] = (PetscBool)(nroots != rpEnd - rpStart);
    214:   MPIU_Allreduce(MPI_IN_PLACE, sub, 2+numFields, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)sf));
    215:   if (sub[0]) {
    216:     ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
    217:     ISGetIndices(selected, &indices);
    218:     PetscSFCreateEmbeddedRootSF(sf, rpEnd - rpStart, indices, &embedSF);
    219:     ISRestoreIndices(selected, &indices);
    220:     ISDestroy(&selected);
    221:   } else {
    222:     PetscObjectReference((PetscObject)sf);
    223:     embedSF = sf;
    224:   }
    225:   PetscSFGetLeafRange(embedSF, &lpStart, &lpEnd);
    226:   lpEnd++;

    228:   PetscSectionSetChart(leafSection, lpStart, lpEnd);

    230:   /* Constrained dof section */
    231:   hasc = sub[1];
    232:   for (f = 0; f < numFields; ++f) hasc = (PetscBool)(hasc || sub[2+f]);

    234:   /* Could fuse these at the cost of copies and extra allocation */
    235:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart],MPI_REPLACE);
    236:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart],MPI_REPLACE);
    237:   if (sub[1]) {
    238:     PetscSectionCheckConstraints_Static(rootSection);
    239:     PetscSectionCheckConstraints_Static(leafSection);
    240:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart],MPI_REPLACE);
    241:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart],MPI_REPLACE);
    242:   }
    243:   for (f = 0; f < numFields; ++f) {
    244:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart],MPI_REPLACE);
    245:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart],MPI_REPLACE);
    246:     if (sub[2+f]) {
    247:       PetscSectionCheckConstraints_Static(rootSection->field[f]);
    248:       PetscSectionCheckConstraints_Static(leafSection->field[f]);
    249:       PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart],MPI_REPLACE);
    250:       PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart],MPI_REPLACE);
    251:     }
    252:   }
    253:   if (remoteOffsets) {
    254:     PetscMalloc1(lpEnd - lpStart, remoteOffsets);
    255:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart],MPI_REPLACE);
    256:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart],MPI_REPLACE);
    257:   }
    258:   PetscSectionSetUp(leafSection);
    259:   if (hasc) { /* need to communicate bcIndices */
    260:     PetscSF  bcSF;
    261:     PetscInt *rOffBc;

    263:     PetscMalloc1(lpEnd - lpStart, &rOffBc);
    264:     if (sub[1]) {
    265:       PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart],MPI_REPLACE);
    266:       PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart],MPI_REPLACE);
    267:       PetscSFCreateSectionSF(embedSF, rootSection->bc, rOffBc, leafSection->bc, &bcSF);
    268:       PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices,MPI_REPLACE);
    269:       PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices,MPI_REPLACE);
    270:       PetscSFDestroy(&bcSF);
    271:     }
    272:     for (f = 0; f < numFields; ++f) {
    273:       if (sub[2+f]) {
    274:         PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart],MPI_REPLACE);
    275:         PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart],MPI_REPLACE);
    276:         PetscSFCreateSectionSF(embedSF, rootSection->field[f]->bc, rOffBc, leafSection->field[f]->bc, &bcSF);
    277:         PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices,MPI_REPLACE);
    278:         PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices,MPI_REPLACE);
    279:         PetscSFDestroy(&bcSF);
    280:       }
    281:     }
    282:     PetscFree(rOffBc);
    283:   }
    284:   PetscSFDestroy(&embedSF);
    285:   PetscFree(sub);
    286:   PetscLogEventEnd(PETSCSF_DistSect,sf,0,0,0);
    287:   return(0);
    288: }

    290: /*@C
    291:   PetscSFCreateRemoteOffsets - Create offsets for point data on remote processes

    293:   Collective on sf

    295:   Input Parameters:
    296: + sf          - The SF
    297: . rootSection - Data layout of remote points for outgoing data (this is layout for SF roots)
    298: - leafSection - Data layout of local points for incoming data  (this is layout for SF leaves)

    300:   Output Parameter:
    301: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL

    303:   Level: developer

    305: .seealso: PetscSFCreate()
    306: @*/
    307: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
    308: {
    309:   PetscSF         embedSF;
    310:   const PetscInt *indices;
    311:   IS              selected;
    312:   PetscInt        numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
    313:   PetscErrorCode  ierr;

    316:   *remoteOffsets = NULL;
    317:   PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
    318:   if (numRoots < 0) return(0);
    319:   PetscLogEventBegin(PETSCSF_RemoteOff,sf,0,0,0);
    320:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
    321:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
    322:   ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
    323:   ISGetIndices(selected, &indices);
    324:   PetscSFCreateEmbeddedRootSF(sf, rpEnd - rpStart, indices, &embedSF);
    325:   ISRestoreIndices(selected, &indices);
    326:   ISDestroy(&selected);
    327:   PetscCalloc1(lpEnd - lpStart, remoteOffsets);
    328:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart],MPI_REPLACE);
    329:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart],MPI_REPLACE);
    330:   PetscSFDestroy(&embedSF);
    331:   PetscLogEventEnd(PETSCSF_RemoteOff,sf,0,0,0);
    332:   return(0);
    333: }

    335: /*@C
    336:   PetscSFCreateSectionSF - Create an expanded SF of dofs, assuming the input SF relates points

    338:   Collective on sf

    340:   Input Parameters:
    341: + sf - The SF
    342: . rootSection - Data layout of remote points for outgoing data (this is usually the serial section)
    343: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
    344: - leafSection - Data layout of local points for incoming data  (this is the distributed section)

    346:   Output Parameters:
    347: - sectionSF - The new SF

    349:   Note: Either rootSection or remoteOffsets can be specified

    351:   Level: advanced

    353: .seealso: PetscSFCreate()
    354: @*/
    355: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
    356: {
    357:   MPI_Comm          comm;
    358:   const PetscInt    *localPoints;
    359:   const PetscSFNode *remotePoints;
    360:   PetscInt          lpStart, lpEnd;
    361:   PetscInt          numRoots, numSectionRoots, numPoints, numIndices = 0;
    362:   PetscInt          *localIndices;
    363:   PetscSFNode       *remoteIndices;
    364:   PetscInt          i, ind;
    365:   PetscErrorCode    ierr;

    373:   PetscObjectGetComm((PetscObject)sf,&comm);
    374:   PetscSFCreate(comm, sectionSF);
    375:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
    376:   PetscSectionGetStorageSize(rootSection, &numSectionRoots);
    377:   PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
    378:   if (numRoots < 0) return(0);
    379:   PetscLogEventBegin(PETSCSF_SectSF,sf,0,0,0);
    380:   for (i = 0; i < numPoints; ++i) {
    381:     PetscInt localPoint = localPoints ? localPoints[i] : i;
    382:     PetscInt dof;

    384:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
    385:       PetscSectionGetDof(leafSection, localPoint, &dof);
    386:       numIndices += dof;
    387:     }
    388:   }
    389:   PetscMalloc1(numIndices, &localIndices);
    390:   PetscMalloc1(numIndices, &remoteIndices);
    391:   /* Create new index graph */
    392:   for (i = 0, ind = 0; i < numPoints; ++i) {
    393:     PetscInt localPoint = localPoints ? localPoints[i] : i;
    394:     PetscInt rank       = remotePoints[i].rank;

    396:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
    397:       PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
    398:       PetscInt loff, dof, d;

    400:       PetscSectionGetOffset(leafSection, localPoint, &loff);
    401:       PetscSectionGetDof(leafSection, localPoint, &dof);
    402:       for (d = 0; d < dof; ++d, ++ind) {
    403:         localIndices[ind]        = loff+d;
    404:         remoteIndices[ind].rank  = rank;
    405:         remoteIndices[ind].index = remoteOffset+d;
    406:       }
    407:     }
    408:   }
    409:   if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %D should be %D", ind, numIndices);
    410:   PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
    411:   PetscSFSetUp(*sectionSF);
    412:   PetscLogEventEnd(PETSCSF_SectSF,sf,0,0,0);
    413:   return(0);
    414: }

    416: /*@C
    417:    PetscSFCreateFromLayouts - Creates a parallel star forest mapping two PetscLayout objects

    419:    Collective

    421:    Input Arguments:
    422: +  rmap - PetscLayout defining the global root space
    423: -  lmap - PetscLayout defining the global leaf space

    425:    Output Arguments:
    426: .  sf - The parallel star forest

    428:    Level: intermediate

    430: .seealso: PetscSFCreate(), PetscLayoutCreate(), PetscSFSetGraphLayout()
    431: @*/
    432: PetscErrorCode PetscSFCreateFromLayouts(PetscLayout rmap, PetscLayout lmap, PetscSF* sf)
    433: {
    435:   PetscInt       i,nroots,nleaves = 0;
    436:   PetscInt       rN, lst, len;
    437:   PetscMPIInt    owner = -1;
    438:   PetscSFNode    *remote;
    439:   MPI_Comm       rcomm = rmap->comm;
    440:   MPI_Comm       lcomm = lmap->comm;
    441:   PetscMPIInt    flg;

    445:   if (!rmap->setupcalled) SETERRQ(rcomm,PETSC_ERR_ARG_WRONGSTATE,"Root layout not setup");
    446:   if (!lmap->setupcalled) SETERRQ(lcomm,PETSC_ERR_ARG_WRONGSTATE,"Leaf layout not setup");
    447:   MPI_Comm_compare(rcomm,lcomm,&flg);
    448:   if (flg != MPI_CONGRUENT && flg != MPI_IDENT) SETERRQ(rcomm,PETSC_ERR_SUP,"cannot map two layouts with non-matching communicators");
    449:   PetscSFCreate(rcomm,sf);
    450:   PetscLayoutGetLocalSize(rmap,&nroots);
    451:   PetscLayoutGetSize(rmap,&rN);
    452:   PetscLayoutGetRange(lmap,&lst,&len);
    453:   PetscMalloc1(len-lst,&remote);
    454:   for (i = lst; i < len && i < rN; i++) {
    455:     if (owner < -1 || i >= rmap->range[owner+1]) {
    456:       PetscLayoutFindOwner(rmap,i,&owner);
    457:     }
    458:     remote[nleaves].rank  = owner;
    459:     remote[nleaves].index = i - rmap->range[owner];
    460:     nleaves++;
    461:   }
    462:   PetscSFSetGraph(*sf,nroots,nleaves,NULL,PETSC_OWN_POINTER,remote,PETSC_COPY_VALUES);
    463:   PetscFree(remote);
    464:   return(0);
    465: }

    467: /* TODO: handle nooffprocentries like MatZeroRowsMapLocal_Private, since this code is the same */
    468: PetscErrorCode PetscLayoutMapLocal(PetscLayout map,PetscInt N,const PetscInt idxs[], PetscInt *on,PetscInt **oidxs,PetscInt **ogidxs)
    469: {
    470:   PetscInt      *owners = map->range;
    471:   PetscInt       n      = map->n;
    472:   PetscSF        sf;
    473:   PetscInt      *lidxs,*work = NULL;
    474:   PetscSFNode   *ridxs;
    475:   PetscMPIInt    rank, p = 0;
    476:   PetscInt       r, len = 0;

    480:   if (on) *on = 0;              /* squelch -Wmaybe-uninitialized */
    481:   /* Create SF where leaves are input idxs and roots are owned idxs */
    482:   MPI_Comm_rank(map->comm,&rank);
    483:   PetscMalloc1(n,&lidxs);
    484:   for (r = 0; r < n; ++r) lidxs[r] = -1;
    485:   PetscMalloc1(N,&ridxs);
    486:   for (r = 0; r < N; ++r) {
    487:     const PetscInt idx = idxs[r];
    488:     if (idx < 0 || map->N <= idx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Index %D out of range [0,%D)",idx,map->N);
    489:     if (idx < owners[p] || owners[p+1] <= idx) { /* short-circuit the search if the last p owns this idx too */
    490:       PetscLayoutFindOwner(map,idx,&p);
    491:     }
    492:     ridxs[r].rank = p;
    493:     ridxs[r].index = idxs[r] - owners[p];
    494:   }
    495:   PetscSFCreate(map->comm,&sf);
    496:   PetscSFSetGraph(sf,n,N,NULL,PETSC_OWN_POINTER,ridxs,PETSC_OWN_POINTER);
    497:   PetscSFReduceBegin(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);
    498:   PetscSFReduceEnd(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);
    499:   if (ogidxs) { /* communicate global idxs */
    500:     PetscInt cum = 0,start,*work2;

    502:     PetscMalloc1(n,&work);
    503:     PetscCalloc1(N,&work2);
    504:     for (r = 0; r < N; ++r) if (idxs[r] >=0) cum++;
    505:     MPI_Scan(&cum,&start,1,MPIU_INT,MPI_SUM,map->comm);
    506:     start -= cum;
    507:     cum = 0;
    508:     for (r = 0; r < N; ++r) if (idxs[r] >=0) work2[r] = start+cum++;
    509:     PetscSFReduceBegin(sf,MPIU_INT,work2,work,MPI_REPLACE);
    510:     PetscSFReduceEnd(sf,MPIU_INT,work2,work,MPI_REPLACE);
    511:     PetscFree(work2);
    512:   }
    513:   PetscSFDestroy(&sf);
    514:   /* Compress and put in indices */
    515:   for (r = 0; r < n; ++r)
    516:     if (lidxs[r] >= 0) {
    517:       if (work) work[len] = work[r];
    518:       lidxs[len++] = r;
    519:     }
    520:   if (on) *on = len;
    521:   if (oidxs) *oidxs = lidxs;
    522:   if (ogidxs) *ogidxs = work;
    523:   return(0);
    524: }

    526: /*@
    527:   PetscSFCreateByMatchingIndices - Create SF by matching root and leaf indices

    529:   Collective

    531:   Input Parameters:
    532: + layout - PetscLayout defining the global index space and the rank that brokers each index
    533: . numRootIndices - size of rootIndices
    534: . rootIndices - PetscInt array of global indices of which this process requests ownership
    535: . rootLocalIndices - root local index permutation (NULL if no permutation)
    536: . rootLocalOffset - offset to be added to root local indices
    537: . numLeafIndices - size of leafIndices
    538: . leafIndices - PetscInt array of global indices with which this process requires data associated
    539: . leafLocalIndices - leaf local index permutation (NULL if no permutation)
    540: - leafLocalOffset - offset to be added to leaf local indices

    542:   Output Parameter:
    543: + sfA - star forest representing the communication pattern from the layout space to the leaf space (NULL if not needed)
    544: - sf - star forest representing the communication pattern from the root space to the leaf space

    546:   Example 1:
    547: $
    548: $  rank             : 0            1            2
    549: $  rootIndices      : [1 0 2]      [3]          [3]
    550: $  rootLocalOffset  : 100          200          300
    551: $  layout           : [0 1]        [2]          [3]
    552: $  leafIndices      : [0]          [2]          [0 3]
    553: $  leafLocalOffset  : 400          500          600
    554: $
    555: would build the following SF:
    556: $
    557: $  [0] 400 <- (0,101)
    558: $  [1] 500 <- (0,102)
    559: $  [2] 600 <- (0,101)
    560: $  [2] 601 <- (2,300)
    561: $
    562:   Example 2:
    563: $
    564: $  rank             : 0               1               2
    565: $  rootIndices      : [1 0 2]         [3]             [3]
    566: $  rootLocalOffset  : 100             200             300
    567: $  layout           : [0 1]           [2]             [3]
    568: $  leafIndices      : rootIndices     rootIndices     rootIndices
    569: $  leafLocalOffset  : rootLocalOffset rootLocalOffset rootLocalOffset
    570: $
    571: would build the following SF:
    572: $
    573: $  [1] 200 <- (2,300)
    574: $
    575:   Example 3:
    576: $
    577: $  No process requests ownership of global index 1, but no process needs it.
    578: $
    579: $  rank             : 0            1            2
    580: $  numRootIndices   : 2            1            1
    581: $  rootIndices      : [0 2]        [3]          [3]
    582: $  rootLocalOffset  : 100          200          300
    583: $  layout           : [0 1]        [2]          [3]
    584: $  numLeafIndices   : 1            1            2
    585: $  leafIndices      : [0]          [2]          [0 3]
    586: $  leafLocalOffset  : 400          500          600
    587: $
    588: would build the following SF:
    589: $
    590: $  [0] 400 <- (0,100)
    591: $  [1] 500 <- (0,101)
    592: $  [2] 600 <- (0,100)
    593: $  [2] 601 <- (2,300)
    594: $

    596:   Notes:
    597:   The layout parameter represents any partitioning of [0, N), where N is the total number of global indices, and its
    598:   local size can be set to PETSC_DECIDE.
    599:   If a global index x lies in the partition owned by process i, each process whose rootIndices contains x requests
    600:   ownership of x and sends its own rank and the local index of x to process i.
    601:   If multiple processes request ownership of x, the one with the highest rank is to own x.
    602:   Process i then broadcasts the ownership information, so that each process whose leafIndices contains x knows the
    603:   ownership information of x.
    604:   The output sf is constructed by associating each leaf point to a root point in this way.

    606:   Suppose there is point data ordered according to the global indices and partitioned according to the given layout.
    607:   The optional output PetscSF object sfA can be used to push such data to leaf points.

    609:   All indices in rootIndices and leafIndices must lie in the layout range. The union (over all processes) of rootIndices
    610:   must cover that of leafIndices, but need not cover the entire layout.

    612:   If (leafIndices, leafLocalIndices, leafLocalOffset) == (rootIndices, rootLocalIndices, rootLocalOffset), the output
    613:   star forest is almost identity, so will only include non-trivial part of the map.

    615:   Developer Notes:
    616:   Current approach of a process of the highest rank gaining the ownership may cause load imbalance; consider using
    617:   hash(rank, root_local_index) as the bid for the ownership determination.

    619:   Level: advanced

    621: .seealso: PetscSFCreate()
    622: @*/
    623: PetscErrorCode PetscSFCreateByMatchingIndices(PetscLayout layout, PetscInt numRootIndices, const PetscInt *rootIndices, const PetscInt *rootLocalIndices, PetscInt rootLocalOffset, PetscInt numLeafIndices, const PetscInt *leafIndices, const PetscInt *leafLocalIndices, PetscInt leafLocalOffset, PetscSF *sfA, PetscSF *sf)
    624: {
    625:   MPI_Comm        comm = layout->comm;
    626:   PetscMPIInt     size, rank;
    627:   PetscSF         sf1;
    628:   PetscSFNode    *owners, *buffer, *iremote;
    629:   PetscInt       *ilocal, nleaves, N, n, i;
    630: #if defined(PETSC_USE_DEBUG)
    631:   PetscInt        N1;
    632: #endif
    633:   PetscBool       flag;
    634:   PetscErrorCode  ierr;

    643:   if (numRootIndices < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "numRootIndices (%D) must be non-negative", numRootIndices);
    644:   if (numLeafIndices < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "numLeafIndices (%D) must be non-negative", numLeafIndices);
    645:   MPI_Comm_size(comm, &size);
    646:   MPI_Comm_rank(comm, &rank);
    647:   PetscLayoutSetUp(layout);
    648:   PetscLayoutGetSize(layout, &N);
    649:   PetscLayoutGetLocalSize(layout, &n);
    650:   flag = (PetscBool)(leafIndices == rootIndices);
    651:   MPIU_Allreduce(MPI_IN_PLACE, &flag, 1, MPIU_BOOL, MPI_LAND, comm);
    652:   if (flag && numLeafIndices != numRootIndices) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "leafIndices == rootIndices, but numLeafIndices (%D) != numRootIndices(%D)", numLeafIndices, numRootIndices);
    653: #if defined(PETSC_USE_DEBUG)
    654:   N1 = PETSC_MIN_INT;
    655:   for (i = 0; i < numRootIndices; i++) if (rootIndices[i] > N1) N1 = rootIndices[i];
    656:   MPI_Allreduce(MPI_IN_PLACE, &N1, 1, MPIU_INT, MPI_MAX, comm);
    657:   if (N1 >= N) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Max. root index (%D) out of layout range [0,%D)", N1, N);
    658:   if (!flag) {
    659:     N1 = PETSC_MIN_INT;
    660:     for (i = 0; i < numLeafIndices; i++) if (leafIndices[i] > N1) N1 = leafIndices[i];
    661:     MPI_Allreduce(MPI_IN_PLACE, &N1, 1, MPIU_INT, MPI_MAX, comm);
    662:     if (N1 >= N) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Max. leaf index (%D) out of layout range [0,%D)", N1, N);
    663:   }
    664: #endif
    665:   /* Reduce: owners -> buffer */
    666:   PetscMalloc1(n, &buffer);
    667:   PetscSFCreate(comm, &sf1);
    668:   PetscSFSetFromOptions(sf1);
    669:   PetscSFSetGraphLayout(sf1, layout, numRootIndices, NULL, PETSC_OWN_POINTER, rootIndices);
    670:   PetscMalloc1(numRootIndices, &owners);
    671:   for (i = 0; i < numRootIndices; ++i) {
    672:     owners[i].rank = rank;
    673:     owners[i].index = rootLocalOffset + (rootLocalIndices ? rootLocalIndices[i] : i);
    674:   }
    675:   for (i = 0; i < n; ++i) {
    676:     buffer[i].index = -1;
    677:     buffer[i].rank = -1;
    678:   }
    679:   PetscSFReduceBegin(sf1, MPIU_2INT, owners, buffer, MPI_MAXLOC);
    680:   PetscSFReduceEnd(sf1, MPIU_2INT, owners, buffer, MPI_MAXLOC);
    681:   /* Bcast: buffer -> owners */
    682:   if (!flag) {
    683:     /* leafIndices is different from rootIndices */
    684:     PetscFree(owners);
    685:     PetscSFSetGraphLayout(sf1, layout, numLeafIndices, NULL, PETSC_OWN_POINTER, leafIndices);
    686:     PetscMalloc1(numLeafIndices, &owners);
    687:   }
    688:   PetscSFBcastBegin(sf1, MPIU_2INT, buffer, owners, MPI_REPLACE);
    689:   PetscSFBcastEnd(sf1, MPIU_2INT, buffer, owners, MPI_REPLACE);
    690:   for (i = 0; i < numLeafIndices; ++i) if (owners[i].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Global point %D was unclaimed", leafIndices[i]);
    691:   PetscFree(buffer);
    692:   if (sfA) {*sfA = sf1;}
    693:   else     {PetscSFDestroy(&sf1);}
    694:   /* Create sf */
    695:   if (flag && rootLocalIndices == leafLocalIndices && leafLocalOffset == rootLocalOffset) {
    696:     /* leaf space == root space */
    697:     for (i = 0, nleaves = 0; i < numLeafIndices; ++i) if (owners[i].rank != rank) ++nleaves;
    698:     PetscMalloc1(nleaves, &ilocal);
    699:     PetscMalloc1(nleaves, &iremote);
    700:     for (i = 0, nleaves = 0; i < numLeafIndices; ++i) {
    701:       if (owners[i].rank != rank) {
    702:         ilocal[nleaves]        = leafLocalOffset + i;
    703:         iremote[nleaves].rank  = owners[i].rank;
    704:         iremote[nleaves].index = owners[i].index;
    705:         ++nleaves;
    706:       }
    707:     }
    708:     PetscFree(owners);
    709:   } else {
    710:     nleaves = numLeafIndices;
    711:     PetscMalloc1(nleaves, &ilocal);
    712:     for (i = 0; i < nleaves; ++i) {ilocal[i] = leafLocalOffset + (leafLocalIndices ? leafLocalIndices[i] : i);}
    713:     iremote = owners;
    714:   }
    715:   PetscSFCreate(comm, sf);
    716:   PetscSFSetFromOptions(*sf);
    717:   PetscSFSetGraph(*sf, rootLocalOffset + numRootIndices, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
    718:   return(0);
    719: }