Actual source code: sfutils.c
1: #include <petsc/private/sfimpl.h>
2: #include <petsc/private/sectionimpl.h>
4: /*@
5: PetscSFSetGraphLayout - Set a `PetscSF` communication pattern using global indices and a `PetscLayout`
7: Collective
9: Input Parameters:
10: + sf - star forest
11: . layout - `PetscLayout` defining the global space for roots, i.e. which roots are owned by each MPI process
12: . nleaves - number of leaf vertices on the current process, each of these references a root on any MPI process
13: . ilocal - locations of leaves in leafdata buffers, pass `NULL` for contiguous storage, that is the locations are in [0,`nleaves`)
14: . localmode - copy mode for `ilocal`
15: - gremote - root vertices in global numbering corresponding to the leaves
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 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: [](sec_petscsf), `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: PetscFunctionBegin;
39: PetscAssertPointer(layout, 2);
40: if (nleaves > 0 && ilocal) PetscAssertPointer(ilocal, 4);
41: if (nleaves > 0) PetscAssertPointer(gremote, 6);
42: PetscCall(PetscLayoutSetUp(layout));
43: PetscCall(PetscLayoutGetLocalSize(layout, &nroots));
44: PetscCall(PetscLayoutGetRanges(layout, &range));
45: PetscCall(PetscMalloc1(nleaves, &remote));
46: if (nleaves) ls = gremote[0] + 1;
47: for (i = 0; i < nleaves; i++) {
48: const PetscInt idx = gremote[i] - ls;
49: if (idx < 0 || idx >= ln) { /* short-circuit the search */
50: PetscCall(PetscLayoutFindOwnerIndex(layout, gremote[i], &lr, &remote[i].index));
51: remote[i].rank = lr;
52: ls = range[lr];
53: ln = range[lr + 1] - ls;
54: } else {
55: remote[i].rank = lr;
56: remote[i].index = idx;
57: }
58: }
59: PetscCall(PetscSFSetGraph(sf, nroots, nleaves, ilocal, localmode, remote, PETSC_OWN_POINTER));
60: PetscFunctionReturn(PETSC_SUCCESS);
61: }
63: /*@C
64: PetscSFGetGraphLayout - Get the global indices and `PetscLayout` that describe a `PetscSF`
66: Collective
68: Input Parameter:
69: . sf - star forest
71: Output Parameters:
72: + layout - `PetscLayout` defining the global space for roots
73: . nleaves - number of leaf vertices on the current process, each of these references a root on any process
74: . ilocal - locations of leaves in leafdata buffers, or `NULL` for contiguous storage
75: - gremote - root vertices in global numbering corresponding to the leaves
77: Level: intermediate
79: Notes:
80: The outputs are such that passing them as inputs to `PetscSFSetGraphLayout()` would lead to the same star forest.
81: The outputs `layout` and `gremote` are freshly created each time this function is called,
82: so they need to be freed (with `PetscLayoutDestroy()` and `PetscFree()`) by the user.
84: .seealso: [](sec_petscsf), `PetscSF`, `PetscSFSetGraphLayout()`, `PetscSFCreate()`, `PetscSFView()`, `PetscSFSetGraph()`, `PetscSFGetGraph()`
85: @*/
86: PetscErrorCode PetscSFGetGraphLayout(PetscSF sf, PetscLayout *layout, PetscInt *nleaves, const PetscInt *ilocal[], PetscInt *gremote[])
87: {
88: PetscInt nr, nl;
89: const PetscSFNode *ir;
90: PetscLayout lt;
92: PetscFunctionBegin;
94: if (layout) PetscAssertPointer(layout, 2);
95: if (nleaves) PetscAssertPointer(nleaves, 3);
96: if (ilocal) PetscAssertPointer(ilocal, 4);
97: if (gremote) PetscAssertPointer(gremote, 5);
98: PetscCall(PetscSFGetGraph(sf, &nr, &nl, ilocal, &ir));
99: PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)sf), nr, PETSC_DECIDE, 1, <));
100: if (gremote) {
101: PetscInt i;
102: const PetscInt *range;
103: PetscInt *gr;
105: PetscCall(PetscLayoutGetRanges(lt, &range));
106: PetscCall(PetscMalloc1(nl, &gr));
107: for (i = 0; i < nl; i++) gr[i] = range[ir[i].rank] + ir[i].index;
108: *gremote = gr;
109: }
110: if (nleaves) *nleaves = nl;
111: if (layout) *layout = lt;
112: else PetscCall(PetscLayoutDestroy(<));
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: /*@
117: PetscSFSetGraphSection - Sets the `PetscSF` graph (communication pattern) encoding the parallel dof overlap based upon the `PetscSection` describing the data layout.
119: Input Parameters:
120: + sf - The `PetscSF`
121: . localSection - `PetscSection` describing the local data layout
122: - globalSection - `PetscSection` describing the global data layout
124: Level: developer
126: .seealso: [](sec_petscsf), `PetscSF`, `PetscSFSetGraph()`, `PetscSFSetGraphLayout()`
127: @*/
128: PetscErrorCode PetscSFSetGraphSection(PetscSF sf, PetscSection localSection, PetscSection globalSection)
129: {
130: MPI_Comm comm;
131: PetscLayout layout;
132: const PetscInt *ranges;
133: PetscInt *local;
134: PetscSFNode *remote;
135: PetscInt pStart, pEnd, p, nroots, nleaves = 0, l;
136: PetscMPIInt size, rank;
138: PetscFunctionBegin;
143: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
144: PetscCallMPI(MPI_Comm_size(comm, &size));
145: PetscCallMPI(MPI_Comm_rank(comm, &rank));
146: PetscCall(PetscSectionGetChart(globalSection, &pStart, &pEnd));
147: PetscCall(PetscSectionGetConstrainedStorageSize(globalSection, &nroots));
148: PetscCall(PetscLayoutCreate(comm, &layout));
149: PetscCall(PetscLayoutSetBlockSize(layout, 1));
150: PetscCall(PetscLayoutSetLocalSize(layout, nroots));
151: PetscCall(PetscLayoutSetUp(layout));
152: PetscCall(PetscLayoutGetRanges(layout, &ranges));
153: for (p = pStart; p < pEnd; ++p) {
154: PetscInt gdof, gcdof;
156: PetscCall(PetscSectionGetDof(globalSection, p, &gdof));
157: PetscCall(PetscSectionGetConstraintDof(globalSection, p, &gcdof));
158: PetscCheck(gcdof <= (gdof < 0 ? -(gdof + 1) : gdof), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " has %" PetscInt_FMT " constraints > %" PetscInt_FMT " dof", p, gcdof, gdof < 0 ? -(gdof + 1) : gdof);
159: nleaves += gdof < 0 ? -(gdof + 1) - gcdof : gdof - gcdof;
160: }
161: PetscCall(PetscMalloc1(nleaves, &local));
162: PetscCall(PetscMalloc1(nleaves, &remote));
163: for (p = pStart, l = 0; p < pEnd; ++p) {
164: const PetscInt *cind;
165: PetscInt dof, cdof, off, gdof, gcdof, goff, gsize, d, c;
167: PetscCall(PetscSectionGetDof(localSection, p, &dof));
168: PetscCall(PetscSectionGetOffset(localSection, p, &off));
169: PetscCall(PetscSectionGetConstraintDof(localSection, p, &cdof));
170: PetscCall(PetscSectionGetConstraintIndices(localSection, p, &cind));
171: PetscCall(PetscSectionGetDof(globalSection, p, &gdof));
172: PetscCall(PetscSectionGetConstraintDof(globalSection, p, &gcdof));
173: PetscCall(PetscSectionGetOffset(globalSection, p, &goff));
174: if (!gdof) continue; /* Censored point */
175: gsize = gdof < 0 ? -(gdof + 1) - gcdof : gdof - gcdof;
176: if (gsize != dof - cdof) {
177: PetscCheck(gsize == dof, comm, PETSC_ERR_ARG_WRONG, "Global dof %" PetscInt_FMT " for point %" PetscInt_FMT " is neither the constrained size %" PetscInt_FMT ", nor the unconstrained %" PetscInt_FMT, gsize, p, dof - cdof, dof);
178: cdof = 0; /* Ignore constraints */
179: }
180: for (d = 0, c = 0; d < dof; ++d) {
181: if ((c < cdof) && (cind[c] == d)) {
182: ++c;
183: continue;
184: }
185: local[l + d - c] = off + d;
186: }
187: PetscCheck(d - c == gsize, comm, PETSC_ERR_ARG_WRONG, "Point %" PetscInt_FMT ": Global dof %" PetscInt_FMT " != %" PetscInt_FMT " size - number of constraints", p, gsize, d - c);
188: if (gdof < 0) {
189: for (d = 0; d < gsize; ++d, ++l) {
190: PetscInt offset = -(goff + 1) + d, ir;
191: PetscMPIInt r;
193: PetscCall(PetscFindInt(offset, size + 1, ranges, &ir));
194: PetscCall(PetscMPIIntCast(ir, &r));
195: if (r < 0) r = -(r + 2);
196: PetscCheck(!(r < 0) && !(r >= size), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " mapped to invalid process %d (%" PetscInt_FMT ", %" PetscInt_FMT ")", p, r, gdof, goff);
197: remote[l].rank = r;
198: remote[l].index = offset - ranges[r];
199: }
200: } else {
201: for (d = 0; d < gsize; ++d, ++l) {
202: remote[l].rank = rank;
203: remote[l].index = goff + d - ranges[rank];
204: }
205: }
206: }
207: PetscCheck(l == nleaves, comm, PETSC_ERR_PLIB, "Iteration error, l %" PetscInt_FMT " != nleaves %" PetscInt_FMT, l, nleaves);
208: PetscCall(PetscLayoutDestroy(&layout));
209: PetscCall(PetscSFSetGraph(sf, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER));
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: /*@C
214: PetscSFDistributeSection - Create a new `PetscSection` reorganized, moving from the root to the leaves of the `PetscSF`
216: Collective
218: Input Parameters:
219: + sf - The `PetscSF`
220: - rootSection - Section defined on root space
222: Output Parameters:
223: + remoteOffsets - root offsets in leaf storage, or `NULL`, its length will be the size of the chart of `leafSection`
224: - leafSection - Section defined on the leaf space
226: Level: advanced
228: Note:
229: Caller must `PetscFree()` `remoteOffsets` if it was requested
231: To distribute data from the `rootSection` to the `leafSection`, see `PetscSFCreateSectionSF()` or `PetscSectionMigrateData()`.
233: Fortran Note:
234: Use `PetscSFDestroyRemoteOffsets()` when `remoteOffsets` is no longer needed.
236: .seealso: [](sec_petscsf), `PetscSF`, `PetscSFCreate()`, `PetscSFCreateSectionSF()`
237: @*/
238: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt *remoteOffsets[], PetscSection leafSection)
239: {
240: PetscSF embedSF;
241: const PetscInt *indices;
242: IS selected;
243: PetscInt numFields, nroots, rpStart, rpEnd, lpStart = PETSC_INT_MAX, lpEnd = -1, f, c;
244: PetscBool *sub, hasc;
246: PetscFunctionBegin;
249: if (remoteOffsets) PetscAssertPointer(remoteOffsets, 3);
251: PetscCall(PetscLogEventBegin(PETSCSF_DistSect, sf, 0, 0, 0));
252: PetscCall(PetscSectionGetNumFields(rootSection, &numFields));
253: if (numFields) {
254: IS perm;
256: /* PetscSectionSetNumFields() calls PetscSectionReset(), which destroys
257: leafSection->perm. To keep this permutation set by the user, we grab
258: the reference before calling PetscSectionSetNumFields() and set it
259: back after. */
260: PetscCall(PetscSectionGetPermutation(leafSection, &perm));
261: PetscCall(PetscObjectReference((PetscObject)perm));
262: PetscCall(PetscSectionSetNumFields(leafSection, numFields));
263: PetscCall(PetscSectionSetPermutation(leafSection, perm));
264: PetscCall(ISDestroy(&perm));
265: }
266: PetscCall(PetscMalloc1(numFields + 2, &sub));
267: sub[1] = rootSection->bc ? PETSC_TRUE : PETSC_FALSE;
268: for (f = 0; f < numFields; ++f) {
269: PetscSectionSym sym, dsym = NULL;
270: const char *name = NULL;
271: PetscInt numComp = 0;
273: sub[2 + f] = rootSection->field[f]->bc ? PETSC_TRUE : PETSC_FALSE;
274: PetscCall(PetscSectionGetFieldComponents(rootSection, f, &numComp));
275: PetscCall(PetscSectionGetFieldName(rootSection, f, &name));
276: PetscCall(PetscSectionGetFieldSym(rootSection, f, &sym));
277: if (sym) PetscCall(PetscSectionSymDistribute(sym, sf, &dsym));
278: PetscCall(PetscSectionSetFieldComponents(leafSection, f, numComp));
279: PetscCall(PetscSectionSetFieldName(leafSection, f, name));
280: PetscCall(PetscSectionSetFieldSym(leafSection, f, dsym));
281: PetscCall(PetscSectionSymDestroy(&dsym));
282: for (c = 0; c < rootSection->numFieldComponents[f]; ++c) {
283: PetscCall(PetscSectionGetComponentName(rootSection, f, c, &name));
284: PetscCall(PetscSectionSetComponentName(leafSection, f, c, name));
285: }
286: }
287: PetscCall(PetscSectionGetChart(rootSection, &rpStart, &rpEnd));
288: PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
289: rpEnd = PetscMin(rpEnd, nroots);
290: rpEnd = PetscMax(rpStart, rpEnd);
291: /* see if we can avoid creating the embedded SF, since it can cost more than an allreduce */
292: sub[0] = (PetscBool)(nroots != rpEnd - rpStart);
293: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, sub, 2 + numFields, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)sf)));
294: if (sub[0]) {
295: PetscCall(ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected));
296: PetscCall(ISGetIndices(selected, &indices));
297: PetscCall(PetscSFCreateEmbeddedRootSF(sf, rpEnd - rpStart, indices, &embedSF));
298: PetscCall(ISRestoreIndices(selected, &indices));
299: PetscCall(ISDestroy(&selected));
300: } else {
301: PetscCall(PetscObjectReference((PetscObject)sf));
302: embedSF = sf;
303: }
304: PetscCall(PetscSFGetLeafRange(embedSF, &lpStart, &lpEnd));
305: lpEnd++;
307: PetscCall(PetscSectionSetChart(leafSection, lpStart, lpEnd));
309: /* Constrained dof section */
310: hasc = sub[1];
311: for (f = 0; f < numFields; ++f) hasc = (PetscBool)(hasc || sub[2 + f]);
313: /* Could fuse these at the cost of copies and extra allocation */
314: PetscCall(PetscSFBcastBegin(embedSF, MPIU_INT, PetscSafePointerPlusOffset(rootSection->atlasDof, -rpStart), PetscSafePointerPlusOffset(leafSection->atlasDof, -lpStart), MPI_REPLACE));
315: PetscCall(PetscSFBcastEnd(embedSF, MPIU_INT, PetscSafePointerPlusOffset(rootSection->atlasDof, -rpStart), PetscSafePointerPlusOffset(leafSection->atlasDof, -lpStart), MPI_REPLACE));
316: if (sub[1]) {
317: PetscCall(PetscSectionCheckConstraints_Private(rootSection));
318: PetscCall(PetscSectionCheckConstraints_Private(leafSection));
319: PetscCall(PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart], MPI_REPLACE));
320: PetscCall(PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart], MPI_REPLACE));
321: }
322: for (f = 0; f < numFields; ++f) {
323: PetscCall(PetscSFBcastBegin(embedSF, MPIU_INT, PetscSafePointerPlusOffset(rootSection->field[f]->atlasDof, -rpStart), PetscSafePointerPlusOffset(leafSection->field[f]->atlasDof, -lpStart), MPI_REPLACE));
324: PetscCall(PetscSFBcastEnd(embedSF, MPIU_INT, PetscSafePointerPlusOffset(rootSection->field[f]->atlasDof, -rpStart), PetscSafePointerPlusOffset(leafSection->field[f]->atlasDof, -lpStart), MPI_REPLACE));
325: if (sub[2 + f]) {
326: PetscCall(PetscSectionCheckConstraints_Private(rootSection->field[f]));
327: PetscCall(PetscSectionCheckConstraints_Private(leafSection->field[f]));
328: PetscCall(PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart], MPI_REPLACE));
329: PetscCall(PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart], MPI_REPLACE));
330: }
331: }
332: if (remoteOffsets) {
333: PetscCall(PetscMalloc1(lpEnd - lpStart, remoteOffsets));
334: PetscCall(PetscSFBcastBegin(embedSF, MPIU_INT, PetscSafePointerPlusOffset(rootSection->atlasOff, -rpStart), PetscSafePointerPlusOffset(*remoteOffsets, -lpStart), MPI_REPLACE));
335: PetscCall(PetscSFBcastEnd(embedSF, MPIU_INT, PetscSafePointerPlusOffset(rootSection->atlasOff, -rpStart), PetscSafePointerPlusOffset(*remoteOffsets, -lpStart), MPI_REPLACE));
336: }
337: PetscCall(PetscSectionInvalidateMaxDof_Internal(leafSection));
338: PetscCall(PetscSectionSetUp(leafSection));
339: if (hasc) { /* need to communicate bcIndices */
340: PetscSF bcSF;
341: PetscInt *rOffBc;
343: PetscCall(PetscMalloc1(lpEnd - lpStart, &rOffBc));
344: if (sub[1]) {
345: PetscCall(PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart], MPI_REPLACE));
346: PetscCall(PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart], MPI_REPLACE));
347: PetscCall(PetscSFCreateSectionSF(embedSF, rootSection->bc, rOffBc, leafSection->bc, &bcSF));
348: PetscCall(PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices, MPI_REPLACE));
349: PetscCall(PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices, MPI_REPLACE));
350: PetscCall(PetscSFDestroy(&bcSF));
351: }
352: for (f = 0; f < numFields; ++f) {
353: if (sub[2 + f]) {
354: PetscCall(PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart], MPI_REPLACE));
355: PetscCall(PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart], MPI_REPLACE));
356: PetscCall(PetscSFCreateSectionSF(embedSF, rootSection->field[f]->bc, rOffBc, leafSection->field[f]->bc, &bcSF));
357: PetscCall(PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices, MPI_REPLACE));
358: PetscCall(PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices, MPI_REPLACE));
359: PetscCall(PetscSFDestroy(&bcSF));
360: }
361: }
362: PetscCall(PetscFree(rOffBc));
363: }
364: PetscCall(PetscSFDestroy(&embedSF));
365: PetscCall(PetscFree(sub));
366: PetscCall(PetscLogEventEnd(PETSCSF_DistSect, sf, 0, 0, 0));
367: PetscFunctionReturn(PETSC_SUCCESS);
368: }
370: /*@C
371: PetscSFCreateRemoteOffsets - Create offsets for point data on remote processes
373: Collective
375: Input Parameters:
376: + sf - The `PetscSF`
377: . rootSection - Data layout of remote points for outgoing data (this is layout for roots)
378: - leafSection - Data layout of local points for incoming data (this is layout for leaves)
380: Output Parameter:
381: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or `NULL`
383: Level: developer
385: Note:
386: Caller must `PetscFree()` `remoteOffsets` if it was requested
388: Fortran Note:
389: Use `PetscSFDestroyRemoteOffsets()` when `remoteOffsets` is no longer needed.
391: .seealso: [](sec_petscsf), `PetscSF`, `PetscSFCreate()`
392: @*/
393: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt *remoteOffsets[])
394: {
395: PetscSF embedSF;
396: const PetscInt *indices;
397: IS selected;
398: PetscInt numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
400: PetscFunctionBegin;
404: PetscAssertPointer(remoteOffsets, 4);
405: *remoteOffsets = NULL;
406: PetscCall(PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL));
407: if (numRoots < 0) PetscFunctionReturn(PETSC_SUCCESS);
408: PetscCall(PetscLogEventBegin(PETSCSF_RemoteOff, sf, 0, 0, 0));
409: PetscCall(PetscSectionGetChart(rootSection, &rpStart, &rpEnd));
410: PetscCall(PetscSectionGetChart(leafSection, &lpStart, &lpEnd));
411: PetscCall(ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected));
412: PetscCall(ISGetIndices(selected, &indices));
413: PetscCall(PetscSFCreateEmbeddedRootSF(sf, rpEnd - rpStart, indices, &embedSF));
414: PetscCall(ISRestoreIndices(selected, &indices));
415: PetscCall(ISDestroy(&selected));
416: PetscCall(PetscCalloc1(lpEnd - lpStart, remoteOffsets));
417: PetscCall(PetscSFBcastBegin(embedSF, MPIU_INT, PetscSafePointerPlusOffset(rootSection->atlasOff, -rpStart), PetscSafePointerPlusOffset(*remoteOffsets, -lpStart), MPI_REPLACE));
418: PetscCall(PetscSFBcastEnd(embedSF, MPIU_INT, PetscSafePointerPlusOffset(rootSection->atlasOff, -rpStart), PetscSafePointerPlusOffset(*remoteOffsets, -lpStart), MPI_REPLACE));
419: PetscCall(PetscSFDestroy(&embedSF));
420: PetscCall(PetscLogEventEnd(PETSCSF_RemoteOff, sf, 0, 0, 0));
421: PetscFunctionReturn(PETSC_SUCCESS);
422: }
424: /*@
425: PetscSFCreateSectionSF - Create an expanded `PetscSF` of dofs, assuming the input `PetscSF` relates points
427: Collective
429: Input Parameters:
430: + sf - The `PetscSF`
431: . rootSection - Data layout of remote points for outgoing data (this is usually the serial section)
432: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or `NULL`
433: - leafSection - Data layout of local points for incoming data (this is the distributed section)
435: Output Parameter:
436: . sectionSF - The new `PetscSF`
438: Level: advanced
440: Notes:
441: `remoteOffsets` can be `NULL` if `sf` does not reference any points in `leafSection`
443: .seealso: [](sec_petscsf), `PetscSF`, `PetscSFCreate()`, `PetscSFDistributeSection()`
444: @*/
445: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
446: {
447: MPI_Comm comm;
448: const PetscInt *localPoints;
449: const PetscSFNode *remotePoints;
450: PetscInt lpStart, lpEnd;
451: PetscInt numRoots, numSectionRoots, numPoints, numIndices = 0;
452: PetscInt *localIndices;
453: PetscSFNode *remoteIndices;
454: PetscInt i, ind;
456: PetscFunctionBegin;
458: PetscAssertPointer(rootSection, 2);
459: /* Cannot check PetscAssertPointer(remoteOffsets,3) because it can be NULL if sf does not reference any points in leafSection */
460: PetscAssertPointer(leafSection, 4);
461: PetscAssertPointer(sectionSF, 5);
462: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
463: PetscCall(PetscSFCreate(comm, sectionSF));
464: PetscCall(PetscSectionGetChart(leafSection, &lpStart, &lpEnd));
465: PetscCall(PetscSectionGetStorageSize(rootSection, &numSectionRoots));
466: PetscCall(PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints));
467: if (numRoots < 0) PetscFunctionReturn(PETSC_SUCCESS);
468: PetscCall(PetscLogEventBegin(PETSCSF_SectSF, sf, 0, 0, 0));
469: for (i = 0; i < numPoints; ++i) {
470: PetscInt localPoint = localPoints ? localPoints[i] : i;
471: PetscInt dof;
473: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
474: PetscCall(PetscSectionGetDof(leafSection, localPoint, &dof));
475: numIndices += dof < 0 ? 0 : dof;
476: }
477: }
478: PetscCall(PetscMalloc1(numIndices, &localIndices));
479: PetscCall(PetscMalloc1(numIndices, &remoteIndices));
480: /* Create new index graph */
481: for (i = 0, ind = 0; i < numPoints; ++i) {
482: PetscInt localPoint = localPoints ? localPoints[i] : i;
483: PetscInt rank = remotePoints[i].rank;
485: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
486: PetscInt remoteOffset = remoteOffsets[localPoint - lpStart];
487: PetscInt loff, dof, d;
489: PetscCall(PetscSectionGetOffset(leafSection, localPoint, &loff));
490: PetscCall(PetscSectionGetDof(leafSection, localPoint, &dof));
491: for (d = 0; d < dof; ++d, ++ind) {
492: localIndices[ind] = loff + d;
493: remoteIndices[ind].rank = rank;
494: remoteIndices[ind].index = remoteOffset + d;
495: }
496: }
497: }
498: PetscCheck(numIndices == ind, comm, PETSC_ERR_PLIB, "Inconsistency in indices, %" PetscInt_FMT " should be %" PetscInt_FMT, ind, numIndices);
499: PetscCall(PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER));
500: PetscCall(PetscSFSetUp(*sectionSF));
501: PetscCall(PetscLogEventEnd(PETSCSF_SectSF, sf, 0, 0, 0));
502: PetscFunctionReturn(PETSC_SUCCESS);
503: }
505: /*@C
506: PetscSFCreateFromLayouts - Creates a parallel star forest mapping between two `PetscLayout` objects
508: Collective
510: Input Parameters:
511: + rmap - `PetscLayout` defining the global root space
512: - lmap - `PetscLayout` defining the global leaf space
514: Output Parameter:
515: . sf - The parallel star forest
517: Level: intermediate
519: Notes:
520: If the global length of `lmap` differs from the global length of `rmap` then the excess entries are ignored.
522: The resulting `sf` used with `PetscSFBcastBegin()` and `PetscSFBcastEnd()` merely copies the array entries of `rootdata` to
523: `leafdata`; moving them between MPI processes if needed. For example,
524: if rmap is [0, 3, 5) and lmap is [0, 2, 6) and `rootdata` is (1, 2, 3) on MPI rank 0 and (4, 5) on MPI rank 1 then the
525: `leafdata` would become (1, 2) on MPI rank 0 and (3, 4, 5, x) on MPI rank 1.
527: .seealso: [](sec_petscsf), `PetscSF`, `PetscLayout`, `PetscSFCreate()`, `PetscSFSetGraph()`, `PetscLayoutCreate()`, `PetscSFSetGraphLayout()`
528: @*/
529: PetscErrorCode PetscSFCreateFromLayouts(PetscLayout rmap, PetscLayout lmap, PetscSF *sf)
530: {
531: PetscInt i, nroots, nleaves = 0;
532: PetscInt rN, lst, len;
533: PetscMPIInt owner = -1;
534: PetscSFNode *remote;
535: MPI_Comm rcomm = rmap->comm;
536: MPI_Comm lcomm = lmap->comm;
537: PetscMPIInt flg;
539: PetscFunctionBegin;
540: PetscAssertPointer(rmap, 1);
541: PetscAssertPointer(lmap, 2);
542: PetscAssertPointer(sf, 3);
543: PetscCheck(rmap->setupcalled, rcomm, PETSC_ERR_ARG_WRONGSTATE, "Root layout not setup");
544: PetscCheck(lmap->setupcalled, lcomm, PETSC_ERR_ARG_WRONGSTATE, "Leaf layout not setup");
545: PetscCallMPI(MPI_Comm_compare(rcomm, lcomm, &flg));
546: PetscCheck(flg == MPI_CONGRUENT || flg == MPI_IDENT, rcomm, PETSC_ERR_SUP, "cannot map two layouts with non-matching communicators");
547: PetscCall(PetscSFCreate(rcomm, sf));
548: PetscCall(PetscLayoutGetLocalSize(rmap, &nroots));
549: PetscCall(PetscLayoutGetSize(rmap, &rN));
550: PetscCall(PetscLayoutGetRange(lmap, &lst, &len));
551: PetscCall(PetscMalloc1(len - lst, &remote));
552: for (i = lst; i < len && i < rN; i++) {
553: if (owner < -1 || i >= rmap->range[owner + 1]) PetscCall(PetscLayoutFindOwner(rmap, i, &owner));
554: remote[nleaves].rank = owner;
555: remote[nleaves].index = i - rmap->range[owner];
556: nleaves++;
557: }
558: PetscCall(PetscSFSetGraph(*sf, nroots, nleaves, NULL, PETSC_OWN_POINTER, remote, PETSC_COPY_VALUES));
559: PetscCall(PetscFree(remote));
560: PetscFunctionReturn(PETSC_SUCCESS);
561: }
563: /*@
564: PetscLayoutMapLocal - Maps a set of global indices to the subset owned locally by each MPI process according to a `PetscLayout`
566: Collective
568: Input Parameters:
569: + map - the `PetscLayout` describing global ownership
570: . N - the number of input indices
571: - idxs - the global indices; negative entries are ignored
573: Output Parameters:
574: + on - number of indices in the returned local set (may be `NULL`)
575: . oidxs - the local (0-based) indices owned by this MPI process (may be `NULL`); caller must `PetscFree()`
576: - ogidxs - the corresponding global indices in a compact numbering (may be `NULL`); caller must `PetscFree()`
578: Level: developer
580: Developer Note:
581: Uses a `PetscSF` internally to route each input index to its owner and reduce with `MPI_LOR`, producing on each MPI
582: process the sorted, deduplicated list of local indices for which at least one process supplied a matching global index.
584: .seealso: `PetscLayout`, `PetscSF`, `PetscLayoutFindOwner()`, `PetscSFCreateFromLayouts()`
585: @*/
586: PetscErrorCode PetscLayoutMapLocal(PetscLayout map, PetscInt N, const PetscInt idxs[], PetscInt *on, PetscInt *oidxs[], PetscInt *ogidxs[])
587: {
588: PetscInt *owners = map->range;
589: PetscInt n = map->n;
590: PetscSF sf;
591: PetscInt *lidxs, *work = NULL, *ilocal;
592: PetscSFNode *ridxs;
593: PetscMPIInt rank, p = 0;
594: PetscInt r, len = 0, nleaves = 0;
596: PetscFunctionBegin;
597: if (on) *on = 0; /* squelch -Wmaybe-uninitialized */
598: /* Create SF where leaves are input idxs and roots are owned idxs */
599: PetscCallMPI(MPI_Comm_rank(map->comm, &rank));
600: PetscCall(PetscMalloc1(n, &lidxs));
601: for (r = 0; r < n; ++r) lidxs[r] = -1;
602: PetscCall(PetscMalloc1(N, &ridxs));
603: PetscCall(PetscMalloc1(N, &ilocal));
604: for (r = 0; r < N; ++r) {
605: const PetscInt idx = idxs[r];
607: if (idx < 0) continue;
608: PetscCheck(idx < map->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Index %" PetscInt_FMT " out of range [0,%" PetscInt_FMT ")", idx, map->N);
609: if (idx < owners[p] || owners[p + 1] <= idx) { /* short-circuit the search if the last p owns this idx too */
610: PetscCall(PetscLayoutFindOwner(map, idx, &p));
611: }
612: ridxs[nleaves].rank = p;
613: ridxs[nleaves].index = idxs[r] - owners[p];
614: ilocal[nleaves] = r;
615: nleaves++;
616: }
617: PetscCall(PetscSFCreate(map->comm, &sf));
618: PetscCall(PetscSFSetGraph(sf, n, nleaves, ilocal, PETSC_OWN_POINTER, ridxs, PETSC_OWN_POINTER));
619: PetscCall(PetscSFReduceBegin(sf, MPIU_INT, (PetscInt *)idxs, lidxs, MPI_LOR));
620: PetscCall(PetscSFReduceEnd(sf, MPIU_INT, (PetscInt *)idxs, lidxs, MPI_LOR));
621: if (ogidxs) { /* communicate global idxs */
622: PetscInt cum = 0, start, *work2;
624: PetscCall(PetscMalloc1(n, &work));
625: PetscCall(PetscCalloc1(N, &work2));
626: for (r = 0; r < N; ++r)
627: if (idxs[r] >= 0) cum++;
628: PetscCallMPI(MPI_Scan(&cum, &start, 1, MPIU_INT, MPI_SUM, map->comm));
629: start -= cum;
630: cum = 0;
631: for (r = 0; r < N; ++r)
632: if (idxs[r] >= 0) work2[r] = start + cum++;
633: PetscCall(PetscSFReduceBegin(sf, MPIU_INT, work2, work, MPI_REPLACE));
634: PetscCall(PetscSFReduceEnd(sf, MPIU_INT, work2, work, MPI_REPLACE));
635: PetscCall(PetscFree(work2));
636: }
637: PetscCall(PetscSFDestroy(&sf));
638: /* Compress and put in indices */
639: for (r = 0; r < n; ++r)
640: if (lidxs[r] >= 0) {
641: if (work) work[len] = work[r];
642: lidxs[len++] = r;
643: }
644: if (on) *on = len;
645: if (oidxs) *oidxs = lidxs;
646: if (ogidxs) *ogidxs = work;
647: PetscFunctionReturn(PETSC_SUCCESS);
648: }
650: /*@
651: PetscSFCreateByMatchingIndices - Create `PetscSF` by matching root and leaf indices
653: Collective
655: Input Parameters:
656: + layout - `PetscLayout` defining the global index space and the MPI rank that brokers each index
657: . numRootIndices - size of `rootIndices`
658: . rootIndices - array of global indices of which this process requests ownership
659: . rootLocalIndices - root local index permutation (`NULL` if no permutation)
660: . rootLocalOffset - offset to be added to `rootLocalIndices`
661: . numLeafIndices - size of `leafIndices`
662: . leafIndices - array of global indices with which this process requires data associated
663: . leafLocalIndices - leaf local index permutation (`NULL` if no permutation)
664: - leafLocalOffset - offset to be added to `leafLocalIndices`
666: Output Parameters:
667: + sfA - star forest representing the communication pattern from the layout space to the leaf space (`NULL` if not needed)
668: - sf - star forest representing the communication pattern from the root space to the leaf space
670: Level: advanced
672: Example 1:
673: .vb
674: rank : 0 1 2
675: rootIndices : [1 0 2] [3] [3]
676: rootLocalOffset : 100 200 300
677: layout : [0 1] [2] [3]
678: leafIndices : [0] [2] [0 3]
679: leafLocalOffset : 400 500 600
681: would build the following PetscSF
683: [0] 400 <- (0,101)
684: [1] 500 <- (0,102)
685: [2] 600 <- (0,101)
686: [2] 601 <- (2,300)
687: .ve
689: Example 2:
690: .vb
691: rank : 0 1 2
692: rootIndices : [1 0 2] [3] [3]
693: rootLocalOffset : 100 200 300
694: layout : [0 1] [2] [3]
695: leafIndices : rootIndices rootIndices rootIndices
696: leafLocalOffset : rootLocalOffset rootLocalOffset rootLocalOffset
698: would build the following PetscSF
700: [1] 200 <- (2,300)
701: .ve
703: Example 3:
704: .vb
705: No process requests ownership of global index 1, but no process needs it.
707: rank : 0 1 2
708: numRootIndices : 2 1 1
709: rootIndices : [0 2] [3] [3]
710: rootLocalOffset : 100 200 300
711: layout : [0 1] [2] [3]
712: numLeafIndices : 1 1 2
713: leafIndices : [0] [2] [0 3]
714: leafLocalOffset : 400 500 600
716: would build the following PetscSF
718: [0] 400 <- (0,100)
719: [1] 500 <- (0,101)
720: [2] 600 <- (0,100)
721: [2] 601 <- (2,300)
722: .ve
724: Notes:
725: `layout` represents any partitioning of [0, N), where N is the total number of global indices, and its
726: local size can be set to `PETSC_DECIDE`.
728: If a global index x lies in the partition owned by process i, each process whose `rootIndices` contains x requests
729: ownership of x and sends its own rank and the local index of x to process i.
730: If multiple processes request ownership of x, the one with the highest rank is to own x.
731: Process i then broadcasts the ownership information, so that each process whose `leafIndices` contains x knows the
732: ownership information of x.
733: The output `sf` is constructed by associating each leaf point to a root point in this way.
735: Suppose there is point data ordered according to the global indices and partitioned according to the given layout.
736: The optional output `sfA` can be used to push such data to leaf points.
738: All indices in `rootIndices` and `leafIndices` must lie in the layout range. The union (over all processes) of `rootIndices`
739: must cover that of `leafIndices`, but need not cover the entire layout.
741: If (leafIndices, leafLocalIndices, leafLocalOffset) == (rootIndices, rootLocalIndices, rootLocalOffset), the output
742: star forest is almost identity, so will only include non-trivial part of the map.
744: Developer Notes:
745: Current approach of a process of the highest rank gaining the ownership may cause load imbalance; consider using
746: hash(rank, root_local_index) as the bid for the ownership determination.
748: .seealso: [](sec_petscsf), `PetscSF`, `PetscSFCreate()`
749: @*/
750: 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)
751: {
752: MPI_Comm comm = layout->comm;
753: PetscMPIInt rank;
754: PetscSF sf1;
755: PetscSFNode *owners, *buffer, *iremote;
756: PetscInt *ilocal, nleaves, N, n, i;
757: PetscBool areIndicesSame;
759: PetscFunctionBegin;
760: PetscAssertPointer(layout, 1);
761: if (rootIndices) PetscAssertPointer(rootIndices, 3);
762: if (rootLocalIndices) PetscAssertPointer(rootLocalIndices, 4);
763: if (leafIndices) PetscAssertPointer(leafIndices, 7);
764: if (leafLocalIndices) PetscAssertPointer(leafLocalIndices, 8);
765: if (sfA) PetscAssertPointer(sfA, 10);
766: PetscAssertPointer(sf, 11);
767: PetscCheck(numRootIndices >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "numRootIndices (%" PetscInt_FMT ") must be non-negative", numRootIndices);
768: PetscCheck(numLeafIndices >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "numLeafIndices (%" PetscInt_FMT ") must be non-negative", numLeafIndices);
769: PetscCallMPI(MPI_Comm_rank(comm, &rank));
770: PetscCall(PetscLayoutSetUp(layout));
771: PetscCall(PetscLayoutGetSize(layout, &N));
772: PetscCall(PetscLayoutGetLocalSize(layout, &n));
773: areIndicesSame = (PetscBool)(leafIndices == rootIndices);
774: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &areIndicesSame, 1, MPI_C_BOOL, MPI_LAND, comm));
775: PetscCheck(!areIndicesSame || numLeafIndices == numRootIndices, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "leafIndices == rootIndices, but numLeafIndices (%" PetscInt_FMT ") != numRootIndices(%" PetscInt_FMT ")", numLeafIndices, numRootIndices);
776: if (PetscDefined(USE_DEBUG)) {
777: PetscInt N1 = PETSC_INT_MIN;
778: for (i = 0; i < numRootIndices; i++)
779: if (rootIndices[i] > N1) N1 = rootIndices[i];
780: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &N1, 1, MPIU_INT, MPI_MAX, comm));
781: PetscCheck(N1 < N, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Max. root index (%" PetscInt_FMT ") out of layout range [0,%" PetscInt_FMT ")", N1, N);
782: if (!areIndicesSame) {
783: N1 = PETSC_INT_MIN;
784: for (i = 0; i < numLeafIndices; i++)
785: if (leafIndices[i] > N1) N1 = leafIndices[i];
786: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &N1, 1, MPIU_INT, MPI_MAX, comm));
787: PetscCheck(N1 < N, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Max. leaf index (%" PetscInt_FMT ") out of layout range [0,%" PetscInt_FMT ")", N1, N);
788: }
789: }
791: /* Reduce: owners -> buffer */
792: PetscCall(PetscMalloc1(n, &buffer));
793: PetscCall(PetscSFCreate(comm, &sf1));
794: PetscCall(PetscSFSetFromOptions(sf1));
795: PetscCall(PetscSFSetGraphLayout(sf1, layout, numRootIndices, NULL, PETSC_OWN_POINTER, rootIndices));
796: PetscCall(PetscMalloc1(numRootIndices, &owners));
797: for (i = 0; i < numRootIndices; ++i) {
798: owners[i].rank = rank;
799: owners[i].index = rootLocalOffset + (rootLocalIndices ? rootLocalIndices[i] : i);
800: }
801: for (i = 0; i < n; ++i) {
802: buffer[i].index = -1;
803: buffer[i].rank = -1;
804: }
805: PetscCall(PetscSFReduceBegin(sf1, MPIU_SF_NODE, owners, buffer, MPI_MAXLOC));
806: PetscCall(PetscSFReduceEnd(sf1, MPIU_SF_NODE, owners, buffer, MPI_MAXLOC));
807: /* Bcast: buffer -> owners */
808: if (!areIndicesSame) {
809: PetscCall(PetscFree(owners));
810: PetscCall(PetscSFSetGraphLayout(sf1, layout, numLeafIndices, NULL, PETSC_OWN_POINTER, leafIndices));
811: PetscCall(PetscMalloc1(numLeafIndices, &owners));
812: }
813: PetscCall(PetscSFBcastBegin(sf1, MPIU_SF_NODE, buffer, owners, MPI_REPLACE));
814: PetscCall(PetscSFBcastEnd(sf1, MPIU_SF_NODE, buffer, owners, MPI_REPLACE));
815: for (i = 0; i < numLeafIndices; ++i) PetscCheck(owners[i].rank >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Global point %" PetscInt_FMT " was unclaimed", leafIndices[i]);
816: PetscCall(PetscFree(buffer));
817: if (sfA) *sfA = sf1;
818: else PetscCall(PetscSFDestroy(&sf1));
819: /* Create sf */
820: if (areIndicesSame && rootLocalIndices == leafLocalIndices && leafLocalOffset == rootLocalOffset) {
821: /* leaf space == root space */
822: for (i = 0, nleaves = 0; i < numLeafIndices; ++i)
823: if (owners[i].rank != rank) ++nleaves;
824: PetscCall(PetscMalloc1(nleaves, &ilocal));
825: PetscCall(PetscMalloc1(nleaves, &iremote));
826: for (i = 0, nleaves = 0; i < numLeafIndices; ++i) {
827: if (owners[i].rank != rank) {
828: ilocal[nleaves] = leafLocalOffset + i;
829: iremote[nleaves].rank = owners[i].rank;
830: iremote[nleaves].index = owners[i].index;
831: ++nleaves;
832: }
833: }
834: PetscCall(PetscFree(owners));
835: } else {
836: nleaves = numLeafIndices;
837: PetscCall(PetscMalloc1(nleaves, &ilocal));
838: for (i = 0; i < nleaves; ++i) ilocal[i] = leafLocalOffset + (leafLocalIndices ? leafLocalIndices[i] : i);
839: iremote = owners;
840: }
841: PetscCall(PetscSFCreate(comm, sf));
842: PetscCall(PetscSFSetFromOptions(*sf));
843: PetscCall(PetscSFSetGraph(*sf, rootLocalOffset + numRootIndices, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
844: PetscFunctionReturn(PETSC_SUCCESS);
845: }
847: /*@
848: PetscSFMerge - append/merge indices of `sfb` into `sfa`, with preference for `sfb`
850: Collective
852: Input Parameters:
853: + sfa - default `PetscSF`
854: - sfb - additional edges to add/replace edges in `sfa`
856: Output Parameter:
857: . merged - new `PetscSF` with combined edges
859: Level: intermediate
861: .seealso: [](sec_petscsf), `PetscSF`, `PetscSFCompose()`
862: @*/
863: PetscErrorCode PetscSFMerge(PetscSF sfa, PetscSF sfb, PetscSF *merged)
864: {
865: PetscInt maxleaf;
867: PetscFunctionBegin;
870: PetscCheckSameComm(sfa, 1, sfb, 2);
871: PetscAssertPointer(merged, 3);
872: {
873: PetscInt aleaf, bleaf;
874: PetscCall(PetscSFGetLeafRange(sfa, NULL, &aleaf));
875: PetscCall(PetscSFGetLeafRange(sfb, NULL, &bleaf));
876: maxleaf = PetscMax(aleaf, bleaf) + 1; // One more than the last index
877: }
878: PetscInt *clocal, aroots, aleaves, broots, bleaves;
879: PetscSFNode *cremote;
880: const PetscInt *alocal, *blocal;
881: const PetscSFNode *aremote, *bremote;
882: PetscCall(PetscMalloc2(maxleaf, &clocal, maxleaf, &cremote));
883: for (PetscInt i = 0; i < maxleaf; i++) clocal[i] = -1;
884: PetscCall(PetscSFGetGraph(sfa, &aroots, &aleaves, &alocal, &aremote));
885: PetscCall(PetscSFGetGraph(sfb, &broots, &bleaves, &blocal, &bremote));
886: PetscCheck(aroots == broots, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Both sfa and sfb must have the same root space");
887: for (PetscInt i = 0; i < aleaves; i++) {
888: PetscInt a = alocal ? alocal[i] : i;
889: clocal[a] = a;
890: cremote[a] = aremote[i];
891: }
892: for (PetscInt i = 0; i < bleaves; i++) {
893: PetscInt b = blocal ? blocal[i] : i;
894: clocal[b] = b;
895: cremote[b] = bremote[i];
896: }
897: PetscInt nleaves = 0;
898: for (PetscInt i = 0; i < maxleaf; i++) {
899: if (clocal[i] < 0) continue;
900: clocal[nleaves] = clocal[i];
901: cremote[nleaves] = cremote[i];
902: nleaves++;
903: }
904: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)sfa), merged));
905: PetscCall(PetscSFSetGraph(*merged, aroots, nleaves, clocal, PETSC_COPY_VALUES, cremote, PETSC_COPY_VALUES));
906: PetscCall(PetscFree2(clocal, cremote));
907: PetscFunctionReturn(PETSC_SUCCESS);
908: }
910: /*@
911: PetscSFCreateStridedSF - Create an `PetscSF` to communicate interleaved blocks of data
913: Collective
915: Input Parameters:
916: + sf - star forest
917: . bs - stride
918: . ldr - leading dimension of root space
919: - ldl - leading dimension of leaf space
921: Output Parameter:
922: . vsf - the new `PetscSF`
924: Level: intermediate
926: Notes:
927: This can be useful to perform communications on multiple right-hand sides stored in a Fortran-style two dimensional array.
928: For example, the calling sequence
929: .vb
930: c_datatype *roots, *leaves;
931: for i in [0,bs) do
932: PetscSFBcastBegin(sf, mpi_datatype, roots + i*ldr, leaves + i*ldl, op)
933: PetscSFBcastEnd(sf, mpi_datatype, roots + i*ldr, leaves + i*ldl, op)
934: .ve
935: is equivalent to
936: .vb
937: c_datatype *roots, *leaves;
938: PetscSFCreateStridedSF(sf, bs, ldr, ldl, &vsf)
939: PetscSFBcastBegin(vsf, mpi_datatype, roots, leaves, op)
940: PetscSFBcastEnd(vsf, mpi_datatype, roots, leaves, op)
941: .ve
943: Developer Notes:
944: Should this functionality be handled with a new API instead of creating a new object?
946: .seealso: [](sec_petscsf), `PetscSF`, `PetscSFCreate()`, `PetscSFSetGraph()`
947: @*/
948: PetscErrorCode PetscSFCreateStridedSF(PetscSF sf, PetscInt bs, PetscInt ldr, PetscInt ldl, PetscSF *vsf)
949: {
950: PetscSF rankssf;
951: const PetscSFNode *iremote, *sfrremote;
952: PetscSFNode *viremote;
953: const PetscInt *ilocal;
954: PetscInt *vilocal = NULL, *ldrs;
955: PetscInt nranks, nr, nl, vnr, vnl, maxl;
956: PetscMPIInt rank;
957: MPI_Comm comm;
958: PetscSFType sftype;
960: PetscFunctionBegin;
963: PetscAssertPointer(vsf, 5);
964: if (bs == 1) {
965: PetscCall(PetscObjectReference((PetscObject)sf));
966: *vsf = sf;
967: PetscFunctionReturn(PETSC_SUCCESS);
968: }
969: PetscCall(PetscSFSetUp(sf));
970: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
971: PetscCallMPI(MPI_Comm_rank(comm, &rank));
972: PetscCall(PetscSFGetGraph(sf, &nr, &nl, &ilocal, &iremote));
973: PetscCall(PetscSFGetLeafRange(sf, NULL, &maxl));
974: maxl += 1;
975: if (ldl == PETSC_DECIDE) ldl = maxl;
976: if (ldr == PETSC_DECIDE) ldr = nr;
977: ldl /= PetscMax(1, sf->vscat.bs); // SFs created from VecScatterCreate() may have a nonzero block size. If not 0, we need to scale ldl and ldr
978: ldr /= PetscMax(1, sf->vscat.bs);
979: PetscCheck(ldr >= nr, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid leading dimension %" PetscInt_FMT " must be smaller than number of roots %" PetscInt_FMT, ldr, nr);
980: PetscCheck(ldl >= maxl, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid leading dimension %" PetscInt_FMT " must be larger than leaf range %" PetscInt_FMT, ldl, maxl - 1);
981: vnr = nr * bs;
982: vnl = nl * bs;
983: PetscCall(PetscMalloc1(vnl, &viremote));
984: PetscCall(PetscMalloc1(vnl, &vilocal));
986: /* Communicate root leading dimensions to leaf ranks */
987: PetscCall(PetscSFGetRanksSF(sf, &rankssf));
988: PetscCall(PetscSFGetGraph(rankssf, NULL, &nranks, NULL, &sfrremote));
989: PetscCall(PetscMalloc1(nranks, &ldrs));
990: PetscCall(PetscSFBcastBegin(rankssf, MPIU_INT, &ldr, ldrs, MPI_REPLACE));
991: PetscCall(PetscSFBcastEnd(rankssf, MPIU_INT, &ldr, ldrs, MPI_REPLACE));
993: for (PetscInt i = 0, rold = -1, lda = -1; i < nl; i++) {
994: const PetscInt r = iremote[i].rank;
995: const PetscInt ii = iremote[i].index;
997: if (r == rank) lda = ldr;
998: else if (rold != r) {
999: PetscInt j;
1001: for (j = 0; j < nranks; j++)
1002: if (sfrremote[j].rank == r) break;
1003: PetscCheck(j < nranks, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unable to locate neighbor rank %" PetscInt_FMT, r);
1004: lda = ldrs[j];
1005: }
1006: rold = r;
1007: for (PetscInt v = 0; v < bs; v++) {
1008: viremote[v * nl + i].rank = r;
1009: viremote[v * nl + i].index = v * lda + ii;
1010: vilocal[v * nl + i] = v * ldl + (ilocal ? ilocal[i] : i);
1011: }
1012: }
1013: PetscCall(PetscFree(ldrs));
1014: PetscCall(PetscSFCreate(comm, vsf));
1015: if (sf->vscat.bs > 1) {
1016: (*vsf)->vscat.bs = sf->vscat.bs;
1017: PetscCallMPI(MPI_Type_dup(sf->vscat.unit, &(*vsf)->vscat.unit));
1018: (*vsf)->vscat.to_n = bs * sf->vscat.to_n;
1019: (*vsf)->vscat.from_n = bs * sf->vscat.from_n;
1020: }
1021: PetscCall(PetscSFGetType(sf, &sftype));
1022: PetscCall(PetscSFSetType(*vsf, sftype));
1023: PetscCall(PetscSFSetGraph(*vsf, vnr, vnl, vilocal, PETSC_OWN_POINTER, viremote, PETSC_OWN_POINTER));
1024: PetscFunctionReturn(PETSC_SUCCESS);
1025: }