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 Parameters:
 10: +  sf - star forest
 11: .  layout - PetscLayout defining the global space for roots
 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 - root vertices in global numbering corresponding to leaves in ilocal

 17:    Level: intermediate

 19:    Notes:
 20:    Global indices must lie in [0, N) where N is the global size of layout.
 21:    Leaf indices in ilocal get sorted; this means the user-provided array gets sorted if localmode is PETSC_OWN_POINTER.

 23:    Developer Notes:
 24:    Local indices which are the identity permutation in the range [0,nleaves) are discarded as they
 25:    encode contiguous storage. In such case, if localmode is PETSC_OWN_POINTER, the memory is deallocated as it is not
 26:    needed

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

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

 57: /*@
 58:   PetscSFSetGraphSection - Sets the PetscSF graph encoding the parallel dof overlap based upon the PetscSections
 59:   describing the data layout.

 61:   Input Parameters:
 62: + sf - The SF
 63: . localSection - PetscSection describing the local data layout
 64: - globalSection - PetscSection describing the global data layout

 66:   Level: developer

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


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

 97:     PetscSectionGetDof(globalSection, p, &gdof);
 98:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
100:     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
101:   }
102:   PetscMalloc1(nleaves, &local);
103:   PetscMalloc1(nleaves, &remote);
104:   for (p = pStart, l = 0; p < pEnd; ++p) {
105:     const PetscInt *cind;
106:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

108:     PetscSectionGetDof(localSection, p, &dof);
109:     PetscSectionGetOffset(localSection, p, &off);
110:     PetscSectionGetConstraintDof(localSection, p, &cdof);
111:     PetscSectionGetConstraintIndices(localSection, p, &cind);
112:     PetscSectionGetDof(globalSection, p, &gdof);
113:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
114:     PetscSectionGetOffset(globalSection, p, &goff);
115:     if (!gdof) continue; /* Censored point */
116:     gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
117:     if (gsize != dof-cdof) {
119:       cdof = 0; /* Ignore constraints */
120:     }
121:     for (d = 0, c = 0; d < dof; ++d) {
122:       if ((c < cdof) && (cind[c] == d)) {++c; continue;}
123:       local[l+d-c] = off+d;
124:     }
126:     if (gdof < 0) {
127:       for (d = 0; d < gsize; ++d, ++l) {
128:         PetscInt offset = -(goff+1) + d, r;

130:         PetscFindInt(offset,size+1,ranges,&r);
131:         if (r < 0) r = -(r+2);
133:         remote[l].rank  = r;
134:         remote[l].index = offset - ranges[r];
135:       }
136:     } else {
137:       for (d = 0; d < gsize; ++d, ++l) {
138:         remote[l].rank  = rank;
139:         remote[l].index = goff+d - ranges[rank];
140:       }
141:     }
142:   }
144:   PetscLayoutDestroy(&layout);
145:   PetscSFSetGraph(sf, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
146:   return 0;
147: }

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

152:   Collective on sf

154:   Input Parameters:
155: + sf - The SF
156: - rootSection - Section defined on root space

158:   Output Parameters:
159: + remoteOffsets - root offsets in leaf storage, or NULL
160: - leafSection - Section defined on the leaf space

162:   Level: advanced

164: .seealso: PetscSFCreate()
165: @*/
166: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
167: {
168:   PetscSF        embedSF;
169:   const PetscInt *indices;
170:   IS             selected;
171:   PetscInt       numFields, nroots, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, f, c;
172:   PetscBool      *sub, hasc;

174:   PetscLogEventBegin(PETSCSF_DistSect,sf,0,0,0);
175:   PetscSectionGetNumFields(rootSection, &numFields);
176:   if (numFields) {
177:     IS perm;

179:     /* PetscSectionSetNumFields() calls PetscSectionReset(), which destroys
180:        leafSection->perm. To keep this permutation set by the user, we grab
181:        the reference before calling PetscSectionSetNumFields() and set it
182:        back after. */
183:     PetscSectionGetPermutation(leafSection, &perm);
184:     PetscObjectReference((PetscObject)perm);
185:     PetscSectionSetNumFields(leafSection, numFields);
186:     PetscSectionSetPermutation(leafSection, perm);
187:     ISDestroy(&perm);
188:   }
189:   PetscMalloc1(numFields+2, &sub);
190:   sub[1] = rootSection->bc ? PETSC_TRUE : PETSC_FALSE;
191:   for (f = 0; f < numFields; ++f) {
192:     PetscSectionSym sym, dsym = NULL;
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:     if (sym) PetscSectionSymDistribute(sym, sf, &dsym);
201:     PetscSectionSetFieldComponents(leafSection, f, numComp);
202:     PetscSectionSetFieldName(leafSection, f, name);
203:     PetscSectionSetFieldSym(leafSection, f, dsym);
204:     PetscSectionSymDestroy(&dsym);
205:     for (c = 0; c < rootSection->numFieldComponents[f]; ++c) {
206:       PetscSectionGetComponentName(rootSection, f, c, &name);
207:       PetscSectionSetComponentName(leafSection, f, c, name);
208:     }
209:   }
210:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
211:   PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);
212:   rpEnd = PetscMin(rpEnd,nroots);
213:   rpEnd = PetscMax(rpStart,rpEnd);
214:   /* see if we can avoid creating the embedded SF, since it can cost more than an allreduce */
215:   sub[0] = (PetscBool)(nroots != rpEnd - rpStart);
216:   MPIU_Allreduce(MPI_IN_PLACE, sub, 2+numFields, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)sf));
217:   if (sub[0]) {
218:     ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
219:     ISGetIndices(selected, &indices);
220:     PetscSFCreateEmbeddedRootSF(sf, rpEnd - rpStart, indices, &embedSF);
221:     ISRestoreIndices(selected, &indices);
222:     ISDestroy(&selected);
223:   } else {
224:     PetscObjectReference((PetscObject)sf);
225:     embedSF = sf;
226:   }
227:   PetscSFGetLeafRange(embedSF, &lpStart, &lpEnd);
228:   lpEnd++;

230:   PetscSectionSetChart(leafSection, lpStart, lpEnd);

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

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

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

293: /*@C
294:   PetscSFCreateRemoteOffsets - Create offsets for point data on remote processes

296:   Collective on sf

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

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

306:   Level: developer

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

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

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

339:   Collective on sf

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

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

350:   Note: Either rootSection or remoteOffsets can be specified

352:   Level: advanced

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

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

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

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

399:       PetscSectionGetOffset(leafSection, localPoint, &loff);
400:       PetscSectionGetDof(leafSection, localPoint, &dof);
401:       for (d = 0; d < dof; ++d, ++ind) {
402:         localIndices[ind]        = loff+d;
403:         remoteIndices[ind].rank  = rank;
404:         remoteIndices[ind].index = remoteOffset+d;
405:       }
406:     }
407:   }
409:   PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
410:   PetscSFSetUp(*sectionSF);
411:   PetscLogEventEnd(PETSCSF_SectSF,sf,0,0,0);
412:   return 0;
413: }

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

418:    Collective

420:    Input Parameters:
421: +  rmap - PetscLayout defining the global root space
422: -  lmap - PetscLayout defining the global leaf space

424:    Output Parameter:
425: .  sf - The parallel star forest

427:    Level: intermediate

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

444:   MPI_Comm_compare(rcomm,lcomm,&flg);
446:   PetscSFCreate(rcomm,sf);
447:   PetscLayoutGetLocalSize(rmap,&nroots);
448:   PetscLayoutGetSize(rmap,&rN);
449:   PetscLayoutGetRange(lmap,&lst,&len);
450:   PetscMalloc1(len-lst,&remote);
451:   for (i = lst; i < len && i < rN; i++) {
452:     if (owner < -1 || i >= rmap->range[owner+1]) {
453:       PetscLayoutFindOwner(rmap,i,&owner);
454:     }
455:     remote[nleaves].rank  = owner;
456:     remote[nleaves].index = i - rmap->range[owner];
457:     nleaves++;
458:   }
459:   PetscSFSetGraph(*sf,nroots,nleaves,NULL,PETSC_OWN_POINTER,remote,PETSC_COPY_VALUES);
460:   PetscFree(remote);
461:   return 0;
462: }

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

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

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

521: /*@
522:   PetscSFCreateByMatchingIndices - Create SF by matching root and leaf indices

524:   Collective

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

537:   Output Parameters:
538: + sfA - star forest representing the communication pattern from the layout space to the leaf space (NULL if not needed)
539: - sf - star forest representing the communication pattern from the root space to the leaf space

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

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

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

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

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

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

614:   Level: advanced

616: .seealso: PetscSFCreate()
617: @*/
618: 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)
619: {
620:   MPI_Comm        comm = layout->comm;
621:   PetscMPIInt     size, rank;
622:   PetscSF         sf1;
623:   PetscSFNode    *owners, *buffer, *iremote;
624:   PetscInt       *ilocal, nleaves, N, n, i;
625: #if defined(PETSC_USE_DEBUG)
626:   PetscInt        N1;
627: #endif
628:   PetscBool       flag;

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