Actual source code: sfutils.c

petsc-3.15.0 2021-04-05
Report Typos and Errors
  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: }