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, <);
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(<);
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: }