Actual source code: sfutils.c

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

  4: /*@
  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: -  gremote - root vertices in global numbering corresponding to leaves in ilocal

 17:    Level: intermediate

 19:    Note:
 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 Note:
 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: `PetscSF`, `PetscSFGetGraphLayout()`, `PetscSFCreate()`, `PetscSFView()`, `PetscSFSetGraph()`, `PetscSFGetGraph()`
 29: @*/
 30: PetscErrorCode PetscSFSetGraphLayout(PetscSF sf, PetscLayout layout, PetscInt nleaves, PetscInt *ilocal, PetscCopyMode localmode, const PetscInt *gremote)
 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 = gremote[0] + 1;
 41:   for (i = 0; i < nleaves; i++) {
 42:     const PetscInt idx = gremote[i] - ls;
 43:     if (idx < 0 || idx >= ln) { /* short-circuit the search */
 44:       PetscLayoutFindOwnerIndex(layout, gremote[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: /*@C
 58:    PetscSFGetGraphLayout - Get the global indices and `PetscLayout` that describe this star forest

 60:    Collective

 62:    Input Parameter:
 63: .  sf - star forest

 65:    Output Parameters:
 66: +  layout - `PetscLayout` defining the global space for roots
 67: .  nleaves - number of leaf vertices on the current process, each of these references a root on any process
 68: .  ilocal - locations of leaves in leafdata buffers, or NULL for contiguous storage
 69: -  gremote - root vertices in global numbering corresponding to leaves in ilocal

 71:    Level: intermediate

 73:    Notes:
 74:    The outputs are such that passing them as inputs to `PetscSFSetGraphLayout()` would lead to the same star forest.
 75:    The outputs layout and gremote are freshly created each time this function is called,
 76:    so they need to be freed by user and cannot be qualified as const.

 78: .seealso: `PetscSF`, `PetscSFSetGraphLayout()`, `PetscSFCreate()`, `PetscSFView()`, `PetscSFSetGraph()`, `PetscSFGetGraph()`
 79: @*/
 80: PetscErrorCode PetscSFGetGraphLayout(PetscSF sf, PetscLayout *layout, PetscInt *nleaves, const PetscInt *ilocal[], PetscInt *gremote[])
 81: {
 82:   PetscInt           nr, nl;
 83:   const PetscSFNode *ir;
 84:   PetscLayout        lt;

 86:   PetscSFGetGraph(sf, &nr, &nl, ilocal, &ir);
 87:   PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)sf), nr, PETSC_DECIDE, 1, &lt);
 88:   if (gremote) {
 89:     PetscInt        i;
 90:     const PetscInt *range;
 91:     PetscInt       *gr;

 93:     PetscLayoutGetRanges(lt, &range);
 94:     PetscMalloc1(nl, &gr);
 95:     for (i = 0; i < nl; i++) gr[i] = range[ir[i].rank] + ir[i].index;
 96:     *gremote = gr;
 97:   }
 98:   if (nleaves) *nleaves = nl;
 99:   if (layout) *layout = lt;
100:   else PetscLayoutDestroy(&lt);
101:   return 0;
102: }

104: /*@
105:   PetscSFSetGraphSection - Sets the `PetscSF` graph encoding the parallel dof overlap based upon the `PetscSection` describing the data layout.

107:   Input Parameters:
108: + sf - The `PetscSF`
109: . localSection - `PetscSection` describing the local data layout
110: - globalSection - `PetscSection` describing the global data layout

112:   Level: developer

114: .seealso: `PetscSF`, `PetscSFSetGraph()`, `PetscSFSetGraphLayout()`
115: @*/
116: PetscErrorCode PetscSFSetGraphSection(PetscSF sf, PetscSection localSection, PetscSection globalSection)
117: {
118:   MPI_Comm        comm;
119:   PetscLayout     layout;
120:   const PetscInt *ranges;
121:   PetscInt       *local;
122:   PetscSFNode    *remote;
123:   PetscInt        pStart, pEnd, p, nroots, nleaves = 0, l;
124:   PetscMPIInt     size, rank;


130:   PetscObjectGetComm((PetscObject)sf, &comm);
131:   MPI_Comm_size(comm, &size);
132:   MPI_Comm_rank(comm, &rank);
133:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
134:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
135:   PetscLayoutCreate(comm, &layout);
136:   PetscLayoutSetBlockSize(layout, 1);
137:   PetscLayoutSetLocalSize(layout, nroots);
138:   PetscLayoutSetUp(layout);
139:   PetscLayoutGetRanges(layout, &ranges);
140:   for (p = pStart; p < pEnd; ++p) {
141:     PetscInt gdof, gcdof;

143:     PetscSectionGetDof(globalSection, p, &gdof);
144:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
146:     nleaves += gdof < 0 ? -(gdof + 1) - gcdof : gdof - gcdof;
147:   }
148:   PetscMalloc1(nleaves, &local);
149:   PetscMalloc1(nleaves, &remote);
150:   for (p = pStart, l = 0; p < pEnd; ++p) {
151:     const PetscInt *cind;
152:     PetscInt        dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

154:     PetscSectionGetDof(localSection, p, &dof);
155:     PetscSectionGetOffset(localSection, p, &off);
156:     PetscSectionGetConstraintDof(localSection, p, &cdof);
157:     PetscSectionGetConstraintIndices(localSection, p, &cind);
158:     PetscSectionGetDof(globalSection, p, &gdof);
159:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
160:     PetscSectionGetOffset(globalSection, p, &goff);
161:     if (!gdof) continue; /* Censored point */
162:     gsize = gdof < 0 ? -(gdof + 1) - gcdof : gdof - gcdof;
163:     if (gsize != dof - cdof) {
165:       cdof = 0; /* Ignore constraints */
166:     }
167:     for (d = 0, c = 0; d < dof; ++d) {
168:       if ((c < cdof) && (cind[c] == d)) {
169:         ++c;
170:         continue;
171:       }
172:       local[l + d - c] = off + d;
173:     }
175:     if (gdof < 0) {
176:       for (d = 0; d < gsize; ++d, ++l) {
177:         PetscInt offset = -(goff + 1) + d, r;

179:         PetscFindInt(offset, size + 1, ranges, &r);
180:         if (r < 0) r = -(r + 2);
182:         remote[l].rank  = r;
183:         remote[l].index = offset - ranges[r];
184:       }
185:     } else {
186:       for (d = 0; d < gsize; ++d, ++l) {
187:         remote[l].rank  = rank;
188:         remote[l].index = goff + d - ranges[rank];
189:       }
190:     }
191:   }
193:   PetscLayoutDestroy(&layout);
194:   PetscSFSetGraph(sf, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
195:   return 0;
196: }

198: /*@C
199:   PetscSFDistributeSection - Create a new `PetscSection` reorganized, moving from the root to the leaves of the `PetscSF`

201:   Collective

203:   Input Parameters:
204: + sf - The `PetscSF`
205: - rootSection - Section defined on root space

207:   Output Parameters:
208: + remoteOffsets - root offsets in leaf storage, or NULL
209: - leafSection - Section defined on the leaf space

211:   Level: advanced

213: .seealso: `PetscSF`, `PetscSFCreate()`
214: @*/
215: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
216: {
217:   PetscSF         embedSF;
218:   const PetscInt *indices;
219:   IS              selected;
220:   PetscInt        numFields, nroots, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, f, c;
221:   PetscBool      *sub, hasc;

223:   PetscLogEventBegin(PETSCSF_DistSect, sf, 0, 0, 0);
224:   PetscSectionGetNumFields(rootSection, &numFields);
225:   if (numFields) {
226:     IS perm;

228:     /* PetscSectionSetNumFields() calls PetscSectionReset(), which destroys
229:        leafSection->perm. To keep this permutation set by the user, we grab
230:        the reference before calling PetscSectionSetNumFields() and set it
231:        back after. */
232:     PetscSectionGetPermutation(leafSection, &perm);
233:     PetscObjectReference((PetscObject)perm);
234:     PetscSectionSetNumFields(leafSection, numFields);
235:     PetscSectionSetPermutation(leafSection, perm);
236:     ISDestroy(&perm);
237:   }
238:   PetscMalloc1(numFields + 2, &sub);
239:   sub[1] = rootSection->bc ? PETSC_TRUE : PETSC_FALSE;
240:   for (f = 0; f < numFields; ++f) {
241:     PetscSectionSym sym, dsym = NULL;
242:     const char     *name    = NULL;
243:     PetscInt        numComp = 0;

245:     sub[2 + f] = rootSection->field[f]->bc ? PETSC_TRUE : PETSC_FALSE;
246:     PetscSectionGetFieldComponents(rootSection, f, &numComp);
247:     PetscSectionGetFieldName(rootSection, f, &name);
248:     PetscSectionGetFieldSym(rootSection, f, &sym);
249:     if (sym) PetscSectionSymDistribute(sym, sf, &dsym);
250:     PetscSectionSetFieldComponents(leafSection, f, numComp);
251:     PetscSectionSetFieldName(leafSection, f, name);
252:     PetscSectionSetFieldSym(leafSection, f, dsym);
253:     PetscSectionSymDestroy(&dsym);
254:     for (c = 0; c < rootSection->numFieldComponents[f]; ++c) {
255:       PetscSectionGetComponentName(rootSection, f, c, &name);
256:       PetscSectionSetComponentName(leafSection, f, c, name);
257:     }
258:   }
259:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
260:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
261:   rpEnd = PetscMin(rpEnd, nroots);
262:   rpEnd = PetscMax(rpStart, rpEnd);
263:   /* see if we can avoid creating the embedded SF, since it can cost more than an allreduce */
264:   sub[0] = (PetscBool)(nroots != rpEnd - rpStart);
265:   MPIU_Allreduce(MPI_IN_PLACE, sub, 2 + numFields, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)sf));
266:   if (sub[0]) {
267:     ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
268:     ISGetIndices(selected, &indices);
269:     PetscSFCreateEmbeddedRootSF(sf, rpEnd - rpStart, indices, &embedSF);
270:     ISRestoreIndices(selected, &indices);
271:     ISDestroy(&selected);
272:   } else {
273:     PetscObjectReference((PetscObject)sf);
274:     embedSF = sf;
275:   }
276:   PetscSFGetLeafRange(embedSF, &lpStart, &lpEnd);
277:   lpEnd++;

279:   PetscSectionSetChart(leafSection, lpStart, lpEnd);

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

285:   /* Could fuse these at the cost of copies and extra allocation */
286:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart], MPI_REPLACE);
287:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart], MPI_REPLACE);
288:   if (sub[1]) {
289:     PetscSectionCheckConstraints_Private(rootSection);
290:     PetscSectionCheckConstraints_Private(leafSection);
291:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart], MPI_REPLACE);
292:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart], MPI_REPLACE);
293:   }
294:   for (f = 0; f < numFields; ++f) {
295:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart], MPI_REPLACE);
296:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart], MPI_REPLACE);
297:     if (sub[2 + f]) {
298:       PetscSectionCheckConstraints_Private(rootSection->field[f]);
299:       PetscSectionCheckConstraints_Private(leafSection->field[f]);
300:       PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart], MPI_REPLACE);
301:       PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart], MPI_REPLACE);
302:     }
303:   }
304:   if (remoteOffsets) {
305:     PetscMalloc1(lpEnd - lpStart, remoteOffsets);
306:     PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart], MPI_REPLACE);
307:     PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart], MPI_REPLACE);
308:   }
309:   PetscSectionInvalidateMaxDof_Internal(leafSection);
310:   PetscSectionSetUp(leafSection);
311:   if (hasc) { /* need to communicate bcIndices */
312:     PetscSF   bcSF;
313:     PetscInt *rOffBc;

315:     PetscMalloc1(lpEnd - lpStart, &rOffBc);
316:     if (sub[1]) {
317:       PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart], MPI_REPLACE);
318:       PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart], MPI_REPLACE);
319:       PetscSFCreateSectionSF(embedSF, rootSection->bc, rOffBc, leafSection->bc, &bcSF);
320:       PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices, MPI_REPLACE);
321:       PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices, MPI_REPLACE);
322:       PetscSFDestroy(&bcSF);
323:     }
324:     for (f = 0; f < numFields; ++f) {
325:       if (sub[2 + f]) {
326:         PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart], MPI_REPLACE);
327:         PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart], MPI_REPLACE);
328:         PetscSFCreateSectionSF(embedSF, rootSection->field[f]->bc, rOffBc, leafSection->field[f]->bc, &bcSF);
329:         PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices, MPI_REPLACE);
330:         PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices, MPI_REPLACE);
331:         PetscSFDestroy(&bcSF);
332:       }
333:     }
334:     PetscFree(rOffBc);
335:   }
336:   PetscSFDestroy(&embedSF);
337:   PetscFree(sub);
338:   PetscLogEventEnd(PETSCSF_DistSect, sf, 0, 0, 0);
339:   return 0;
340: }

342: /*@C
343:   PetscSFCreateRemoteOffsets - Create offsets for point data on remote processes

345:   Collective

347:   Input Parameters:
348: + sf          - The `PetscSF`
349: . rootSection - Data layout of remote points for outgoing data (this is layout for SF roots)
350: - leafSection - Data layout of local points for incoming data  (this is layout for SF leaves)

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

355:   Level: developer

357: .seealso: `PetscSF`, `PetscSFCreate()`
358: @*/
359: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
360: {
361:   PetscSF         embedSF;
362:   const PetscInt *indices;
363:   IS              selected;
364:   PetscInt        numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;

366:   *remoteOffsets = NULL;
367:   PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
368:   if (numRoots < 0) return 0;
369:   PetscLogEventBegin(PETSCSF_RemoteOff, sf, 0, 0, 0);
370:   PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
371:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
372:   ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
373:   ISGetIndices(selected, &indices);
374:   PetscSFCreateEmbeddedRootSF(sf, rpEnd - rpStart, indices, &embedSF);
375:   ISRestoreIndices(selected, &indices);
376:   ISDestroy(&selected);
377:   PetscCalloc1(lpEnd - lpStart, remoteOffsets);
378:   PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart], MPI_REPLACE);
379:   PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart], MPI_REPLACE);
380:   PetscSFDestroy(&embedSF);
381:   PetscLogEventEnd(PETSCSF_RemoteOff, sf, 0, 0, 0);
382:   return 0;
383: }

385: /*@C
386:   PetscSFCreateSectionSF - Create an expanded `PetscSF` of dofs, assuming the input `PetscSF` relates points

388:   Collective

390:   Input Parameters:
391: + sf - The `PetscSF`
392: . rootSection - Data layout of remote points for outgoing data (this is usually the serial section)
393: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
394: - leafSection - Data layout of local points for incoming data  (this is the distributed section)

396:   Output Parameters:
397: - sectionSF - The new `PetscSF`

399:   Level: advanced

401:   Note:
402:   Either rootSection or remoteOffsets can be specified

404: .seealso:  `PetscSF`, `PetscSFCreate()`
405: @*/
406: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
407: {
408:   MPI_Comm           comm;
409:   const PetscInt    *localPoints;
410:   const PetscSFNode *remotePoints;
411:   PetscInt           lpStart, lpEnd;
412:   PetscInt           numRoots, numSectionRoots, numPoints, numIndices = 0;
413:   PetscInt          *localIndices;
414:   PetscSFNode       *remoteIndices;
415:   PetscInt           i, ind;

422:   PetscObjectGetComm((PetscObject)sf, &comm);
423:   PetscSFCreate(comm, sectionSF);
424:   PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
425:   PetscSectionGetStorageSize(rootSection, &numSectionRoots);
426:   PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
427:   if (numRoots < 0) return 0;
428:   PetscLogEventBegin(PETSCSF_SectSF, sf, 0, 0, 0);
429:   for (i = 0; i < numPoints; ++i) {
430:     PetscInt localPoint = localPoints ? localPoints[i] : i;
431:     PetscInt dof;

433:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
434:       PetscSectionGetDof(leafSection, localPoint, &dof);
435:       numIndices += dof < 0 ? 0 : dof;
436:     }
437:   }
438:   PetscMalloc1(numIndices, &localIndices);
439:   PetscMalloc1(numIndices, &remoteIndices);
440:   /* Create new index graph */
441:   for (i = 0, ind = 0; i < numPoints; ++i) {
442:     PetscInt localPoint = localPoints ? localPoints[i] : i;
443:     PetscInt rank       = remotePoints[i].rank;

445:     if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
446:       PetscInt remoteOffset = remoteOffsets[localPoint - lpStart];
447:       PetscInt loff, dof, d;

449:       PetscSectionGetOffset(leafSection, localPoint, &loff);
450:       PetscSectionGetDof(leafSection, localPoint, &dof);
451:       for (d = 0; d < dof; ++d, ++ind) {
452:         localIndices[ind]        = loff + d;
453:         remoteIndices[ind].rank  = rank;
454:         remoteIndices[ind].index = remoteOffset + d;
455:       }
456:     }
457:   }
459:   PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
460:   PetscSFSetUp(*sectionSF);
461:   PetscLogEventEnd(PETSCSF_SectSF, sf, 0, 0, 0);
462:   return 0;
463: }

465: /*@C
466:    PetscSFCreateFromLayouts - Creates a parallel star forest mapping two `PetscLayout` objects

468:    Collective

470:    Input Parameters:
471: +  rmap - `PetscLayout` defining the global root space
472: -  lmap - `PetscLayout` defining the global leaf space

474:    Output Parameter:
475: .  sf - The parallel star forest

477:    Level: intermediate

479: .seealso: `PetscSF`, `PetscSFCreate()`, `PetscLayoutCreate()`, `PetscSFSetGraphLayout()`
480: @*/
481: PetscErrorCode PetscSFCreateFromLayouts(PetscLayout rmap, PetscLayout lmap, PetscSF *sf)
482: {
483:   PetscInt     i, nroots, nleaves = 0;
484:   PetscInt     rN, lst, len;
485:   PetscMPIInt  owner = -1;
486:   PetscSFNode *remote;
487:   MPI_Comm     rcomm = rmap->comm;
488:   MPI_Comm     lcomm = lmap->comm;
489:   PetscMPIInt  flg;

494:   MPI_Comm_compare(rcomm, lcomm, &flg);
496:   PetscSFCreate(rcomm, sf);
497:   PetscLayoutGetLocalSize(rmap, &nroots);
498:   PetscLayoutGetSize(rmap, &rN);
499:   PetscLayoutGetRange(lmap, &lst, &len);
500:   PetscMalloc1(len - lst, &remote);
501:   for (i = lst; i < len && i < rN; i++) {
502:     if (owner < -1 || i >= rmap->range[owner + 1]) PetscLayoutFindOwner(rmap, i, &owner);
503:     remote[nleaves].rank  = owner;
504:     remote[nleaves].index = i - rmap->range[owner];
505:     nleaves++;
506:   }
507:   PetscSFSetGraph(*sf, nroots, nleaves, NULL, PETSC_OWN_POINTER, remote, PETSC_COPY_VALUES);
508:   PetscFree(remote);
509:   return 0;
510: }

512: /* TODO: handle nooffprocentries like MatZeroRowsMapLocal_Private, since this code is the same */
513: PetscErrorCode PetscLayoutMapLocal(PetscLayout map, PetscInt N, const PetscInt idxs[], PetscInt *on, PetscInt **oidxs, PetscInt **ogidxs)
514: {
515:   PetscInt    *owners = map->range;
516:   PetscInt     n      = map->n;
517:   PetscSF      sf;
518:   PetscInt    *lidxs, *work = NULL;
519:   PetscSFNode *ridxs;
520:   PetscMPIInt  rank, p = 0;
521:   PetscInt     r, len  = 0;

523:   if (on) *on = 0; /* squelch -Wmaybe-uninitialized */
524:   /* Create SF where leaves are input idxs and roots are owned idxs */
525:   MPI_Comm_rank(map->comm, &rank);
526:   PetscMalloc1(n, &lidxs);
527:   for (r = 0; r < n; ++r) lidxs[r] = -1;
528:   PetscMalloc1(N, &ridxs);
529:   for (r = 0; r < N; ++r) {
530:     const PetscInt idx = idxs[r];
532:     if (idx < owners[p] || owners[p + 1] <= idx) { /* short-circuit the search if the last p owns this idx too */
533:       PetscLayoutFindOwner(map, idx, &p);
534:     }
535:     ridxs[r].rank  = p;
536:     ridxs[r].index = idxs[r] - owners[p];
537:   }
538:   PetscSFCreate(map->comm, &sf);
539:   PetscSFSetGraph(sf, n, N, NULL, PETSC_OWN_POINTER, ridxs, PETSC_OWN_POINTER);
540:   PetscSFReduceBegin(sf, MPIU_INT, (PetscInt *)idxs, lidxs, MPI_LOR);
541:   PetscSFReduceEnd(sf, MPIU_INT, (PetscInt *)idxs, lidxs, MPI_LOR);
542:   if (ogidxs) { /* communicate global idxs */
543:     PetscInt cum = 0, start, *work2;

545:     PetscMalloc1(n, &work);
546:     PetscCalloc1(N, &work2);
547:     for (r = 0; r < N; ++r)
548:       if (idxs[r] >= 0) cum++;
549:     MPI_Scan(&cum, &start, 1, MPIU_INT, MPI_SUM, map->comm);
550:     start -= cum;
551:     cum = 0;
552:     for (r = 0; r < N; ++r)
553:       if (idxs[r] >= 0) work2[r] = start + cum++;
554:     PetscSFReduceBegin(sf, MPIU_INT, work2, work, MPI_REPLACE);
555:     PetscSFReduceEnd(sf, MPIU_INT, work2, work, MPI_REPLACE);
556:     PetscFree(work2);
557:   }
558:   PetscSFDestroy(&sf);
559:   /* Compress and put in indices */
560:   for (r = 0; r < n; ++r)
561:     if (lidxs[r] >= 0) {
562:       if (work) work[len] = work[r];
563:       lidxs[len++] = r;
564:     }
565:   if (on) *on = len;
566:   if (oidxs) *oidxs = lidxs;
567:   if (ogidxs) *ogidxs = work;
568:   return 0;
569: }

571: /*@
572:   PetscSFCreateByMatchingIndices - Create `PetscSF` by matching root and leaf indices

574:   Collective

576:   Input Parameters:
577: + layout - `PetscLayout` defining the global index space and the rank that brokers each index
578: . numRootIndices - size of rootIndices
579: . rootIndices - `PetscInt` array of global indices of which this process requests ownership
580: . rootLocalIndices - root local index permutation (NULL if no permutation)
581: . rootLocalOffset - offset to be added to root local indices
582: . numLeafIndices - size of leafIndices
583: . leafIndices - `PetscInt` array of global indices with which this process requires data associated
584: . leafLocalIndices - leaf local index permutation (NULL if no permutation)
585: - leafLocalOffset - offset to be added to leaf local indices

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

591:   Level: advanced

593:   Example 1:
594: .vb
595:   rank             : 0            1            2
596:   rootIndices      : [1 0 2]      [3]          [3]
597:   rootLocalOffset  : 100          200          300
598:   layout           : [0 1]        [2]          [3]
599:   leafIndices      : [0]          [2]          [0 3]
600:   leafLocalOffset  : 400          500          600

602: would build the following PetscSF

604:   [0] 400 <- (0,101)
605:   [1] 500 <- (0,102)
606:   [2] 600 <- (0,101)
607:   [2] 601 <- (2,300)
608: .ve

610:   Example 2:
611: .vb
612:   rank             : 0               1               2
613:   rootIndices      : [1 0 2]         [3]             [3]
614:   rootLocalOffset  : 100             200             300
615:   layout           : [0 1]           [2]             [3]
616:   leafIndices      : rootIndices     rootIndices     rootIndices
617:   leafLocalOffset  : rootLocalOffset rootLocalOffset rootLocalOffset

619: would build the following PetscSF

621:   [1] 200 <- (2,300)
622: .ve

624:   Example 3:
625: .vb
626:   No process requests ownership of global index 1, but no process needs it.

628:   rank             : 0            1            2
629:   numRootIndices   : 2            1            1
630:   rootIndices      : [0 2]        [3]          [3]
631:   rootLocalOffset  : 100          200          300
632:   layout           : [0 1]        [2]          [3]
633:   numLeafIndices   : 1            1            2
634:   leafIndices      : [0]          [2]          [0 3]
635:   leafLocalOffset  : 400          500          600

637: would build the following PetscSF

639:   [0] 400 <- (0,100)
640:   [1] 500 <- (0,101)
641:   [2] 600 <- (0,100)
642:   [2] 601 <- (2,300)
643: .ve

645:   Notes:
646:   The layout parameter represents any partitioning of [0, N), where N is the total number of global indices, and its
647:   local size can be set to `PETSC_DECIDE`.

649:   If a global index x lies in the partition owned by process i, each process whose rootIndices contains x requests
650:   ownership of x and sends its own rank and the local index of x to process i.
651:   If multiple processes request ownership of x, the one with the highest rank is to own x.
652:   Process i then broadcasts the ownership information, so that each process whose leafIndices contains x knows the
653:   ownership information of x.
654:   The output sf is constructed by associating each leaf point to a root point in this way.

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

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

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

665:   Developer Note:
666:   Current approach of a process of the highest rank gaining the ownership may cause load imbalance; consider using
667:   hash(rank, root_local_index) as the bid for the ownership determination.

669: .seealso: `PetscSF`, `PetscSFCreate()`
670: @*/
671: 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)
672: {
673:   MPI_Comm     comm = layout->comm;
674:   PetscMPIInt  size, rank;
675:   PetscSF      sf1;
676:   PetscSFNode *owners, *buffer, *iremote;
677:   PetscInt    *ilocal, nleaves, N, n, i;
678: #if defined(PETSC_USE_DEBUG)
679:   PetscInt N1;
680: #endif
681:   PetscBool flag;

691:   MPI_Comm_size(comm, &size);
692:   MPI_Comm_rank(comm, &rank);
693:   PetscLayoutSetUp(layout);
694:   PetscLayoutGetSize(layout, &N);
695:   PetscLayoutGetLocalSize(layout, &n);
696:   flag = (PetscBool)(leafIndices == rootIndices);
697:   MPIU_Allreduce(MPI_IN_PLACE, &flag, 1, MPIU_BOOL, MPI_LAND, comm);
699: #if defined(PETSC_USE_DEBUG)
700:   N1 = PETSC_MIN_INT;
701:   for (i = 0; i < numRootIndices; i++)
702:     if (rootIndices[i] > N1) N1 = rootIndices[i];
703:   MPI_Allreduce(MPI_IN_PLACE, &N1, 1, MPIU_INT, MPI_MAX, comm);
705:   if (!flag) {
706:     N1 = PETSC_MIN_INT;
707:     for (i = 0; i < numLeafIndices; i++)
708:       if (leafIndices[i] > N1) N1 = leafIndices[i];
709:     MPI_Allreduce(MPI_IN_PLACE, &N1, 1, MPIU_INT, MPI_MAX, comm);
711:   }
712: #endif
713:   /* Reduce: owners -> buffer */
714:   PetscMalloc1(n, &buffer);
715:   PetscSFCreate(comm, &sf1);
716:   PetscSFSetFromOptions(sf1);
717:   PetscSFSetGraphLayout(sf1, layout, numRootIndices, NULL, PETSC_OWN_POINTER, rootIndices);
718:   PetscMalloc1(numRootIndices, &owners);
719:   for (i = 0; i < numRootIndices; ++i) {
720:     owners[i].rank  = rank;
721:     owners[i].index = rootLocalOffset + (rootLocalIndices ? rootLocalIndices[i] : i);
722:   }
723:   for (i = 0; i < n; ++i) {
724:     buffer[i].index = -1;
725:     buffer[i].rank  = -1;
726:   }
727:   PetscSFReduceBegin(sf1, MPIU_2INT, owners, buffer, MPI_MAXLOC);
728:   PetscSFReduceEnd(sf1, MPIU_2INT, owners, buffer, MPI_MAXLOC);
729:   /* Bcast: buffer -> owners */
730:   if (!flag) {
731:     /* leafIndices is different from rootIndices */
732:     PetscFree(owners);
733:     PetscSFSetGraphLayout(sf1, layout, numLeafIndices, NULL, PETSC_OWN_POINTER, leafIndices);
734:     PetscMalloc1(numLeafIndices, &owners);
735:   }
736:   PetscSFBcastBegin(sf1, MPIU_2INT, buffer, owners, MPI_REPLACE);
737:   PetscSFBcastEnd(sf1, MPIU_2INT, buffer, owners, MPI_REPLACE);
739:   PetscFree(buffer);
740:   if (sfA) {
741:     *sfA = sf1;
742:   } else PetscSFDestroy(&sf1);
743:   /* Create sf */
744:   if (flag && rootLocalIndices == leafLocalIndices && leafLocalOffset == rootLocalOffset) {
745:     /* leaf space == root space */
746:     for (i = 0, nleaves = 0; i < numLeafIndices; ++i)
747:       if (owners[i].rank != rank) ++nleaves;
748:     PetscMalloc1(nleaves, &ilocal);
749:     PetscMalloc1(nleaves, &iremote);
750:     for (i = 0, nleaves = 0; i < numLeafIndices; ++i) {
751:       if (owners[i].rank != rank) {
752:         ilocal[nleaves]        = leafLocalOffset + i;
753:         iremote[nleaves].rank  = owners[i].rank;
754:         iremote[nleaves].index = owners[i].index;
755:         ++nleaves;
756:       }
757:     }
758:     PetscFree(owners);
759:   } else {
760:     nleaves = numLeafIndices;
761:     PetscMalloc1(nleaves, &ilocal);
762:     for (i = 0; i < nleaves; ++i) ilocal[i] = leafLocalOffset + (leafLocalIndices ? leafLocalIndices[i] : i);
763:     iremote = owners;
764:   }
765:   PetscSFCreate(comm, sf);
766:   PetscSFSetFromOptions(*sf);
767:   PetscSFSetGraph(*sf, rootLocalOffset + numRootIndices, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
768:   return 0;
769: }