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 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: `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;
38: PetscCall(PetscLayoutSetUp(layout));
39: PetscCall(PetscLayoutGetLocalSize(layout, &nroots));
40: PetscCall(PetscLayoutGetRanges(layout, &range));
41: PetscCall(PetscMalloc1(nleaves, &remote));
42: if (nleaves) ls = gremote[0] + 1;
43: for (i = 0; i < nleaves; i++) {
44: const PetscInt idx = gremote[i] - ls;
45: if (idx < 0 || idx >= ln) { /* short-circuit the search */
46: PetscCall(PetscLayoutFindOwnerIndex(layout, gremote[i], &lr, &remote[i].index));
47: remote[i].rank = lr;
48: ls = range[lr];
49: ln = range[lr + 1] - ls;
50: } else {
51: remote[i].rank = lr;
52: remote[i].index = idx;
53: }
54: }
55: PetscCall(PetscSFSetGraph(sf, nroots, nleaves, ilocal, localmode, remote, PETSC_OWN_POINTER));
56: PetscFunctionReturn(PETSC_SUCCESS);
57: }
59: /*@C
60: PetscSFGetGraphLayout - Get the global indices and `PetscLayout` that describe this star forest
62: Collective
64: Input Parameter:
65: . sf - star forest
67: Output Parameters:
68: + layout - `PetscLayout` defining the global space for roots
69: . nleaves - number of leaf vertices on the current process, each of these references a root on any process
70: . ilocal - locations of leaves in leafdata buffers, or NULL for contiguous storage
71: - gremote - root vertices in global numbering corresponding to leaves in ilocal
73: Level: intermediate
75: Notes:
76: The outputs are such that passing them as inputs to `PetscSFSetGraphLayout()` would lead to the same star forest.
77: The outputs layout and gremote are freshly created each time this function is called,
78: so they need to be freed by user and cannot be qualified as const.
80: .seealso: `PetscSF`, `PetscSFSetGraphLayout()`, `PetscSFCreate()`, `PetscSFView()`, `PetscSFSetGraph()`, `PetscSFGetGraph()`
81: @*/
82: PetscErrorCode PetscSFGetGraphLayout(PetscSF sf, PetscLayout *layout, PetscInt *nleaves, const PetscInt *ilocal[], PetscInt *gremote[])
83: {
84: PetscInt nr, nl;
85: const PetscSFNode *ir;
86: PetscLayout lt;
88: PetscFunctionBegin;
89: PetscCall(PetscSFGetGraph(sf, &nr, &nl, ilocal, &ir));
90: PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)sf), nr, PETSC_DECIDE, 1, <));
91: if (gremote) {
92: PetscInt i;
93: const PetscInt *range;
94: PetscInt *gr;
96: PetscCall(PetscLayoutGetRanges(lt, &range));
97: PetscCall(PetscMalloc1(nl, &gr));
98: for (i = 0; i < nl; i++) gr[i] = range[ir[i].rank] + ir[i].index;
99: *gremote = gr;
100: }
101: if (nleaves) *nleaves = nl;
102: if (layout) *layout = lt;
103: else PetscCall(PetscLayoutDestroy(<));
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
107: /*@
108: PetscSFSetGraphSection - Sets the `PetscSF` graph encoding the parallel dof overlap based upon the `PetscSection` describing the data layout.
110: Input Parameters:
111: + sf - The `PetscSF`
112: . localSection - `PetscSection` describing the local data layout
113: - globalSection - `PetscSection` describing the global data layout
115: Level: developer
117: .seealso: `PetscSF`, `PetscSFSetGraph()`, `PetscSFSetGraphLayout()`
118: @*/
119: PetscErrorCode PetscSFSetGraphSection(PetscSF sf, PetscSection localSection, PetscSection globalSection)
120: {
121: MPI_Comm comm;
122: PetscLayout layout;
123: const PetscInt *ranges;
124: PetscInt *local;
125: PetscSFNode *remote;
126: PetscInt pStart, pEnd, p, nroots, nleaves = 0, l;
127: PetscMPIInt size, rank;
129: PetscFunctionBegin;
134: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
135: PetscCallMPI(MPI_Comm_size(comm, &size));
136: PetscCallMPI(MPI_Comm_rank(comm, &rank));
137: PetscCall(PetscSectionGetChart(globalSection, &pStart, &pEnd));
138: PetscCall(PetscSectionGetConstrainedStorageSize(globalSection, &nroots));
139: PetscCall(PetscLayoutCreate(comm, &layout));
140: PetscCall(PetscLayoutSetBlockSize(layout, 1));
141: PetscCall(PetscLayoutSetLocalSize(layout, nroots));
142: PetscCall(PetscLayoutSetUp(layout));
143: PetscCall(PetscLayoutGetRanges(layout, &ranges));
144: for (p = pStart; p < pEnd; ++p) {
145: PetscInt gdof, gcdof;
147: PetscCall(PetscSectionGetDof(globalSection, p, &gdof));
148: PetscCall(PetscSectionGetConstraintDof(globalSection, p, &gcdof));
149: 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));
150: nleaves += gdof < 0 ? -(gdof + 1) - gcdof : gdof - gcdof;
151: }
152: PetscCall(PetscMalloc1(nleaves, &local));
153: PetscCall(PetscMalloc1(nleaves, &remote));
154: for (p = pStart, l = 0; p < pEnd; ++p) {
155: const PetscInt *cind;
156: PetscInt dof, cdof, off, gdof, gcdof, goff, gsize, d, c;
158: PetscCall(PetscSectionGetDof(localSection, p, &dof));
159: PetscCall(PetscSectionGetOffset(localSection, p, &off));
160: PetscCall(PetscSectionGetConstraintDof(localSection, p, &cdof));
161: PetscCall(PetscSectionGetConstraintIndices(localSection, p, &cind));
162: PetscCall(PetscSectionGetDof(globalSection, p, &gdof));
163: PetscCall(PetscSectionGetConstraintDof(globalSection, p, &gcdof));
164: PetscCall(PetscSectionGetOffset(globalSection, p, &goff));
165: if (!gdof) continue; /* Censored point */
166: gsize = gdof < 0 ? -(gdof + 1) - gcdof : gdof - gcdof;
167: if (gsize != dof - cdof) {
168: 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);
169: cdof = 0; /* Ignore constraints */
170: }
171: for (d = 0, c = 0; d < dof; ++d) {
172: if ((c < cdof) && (cind[c] == d)) {
173: ++c;
174: continue;
175: }
176: local[l + d - c] = off + d;
177: }
178: 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);
179: if (gdof < 0) {
180: for (d = 0; d < gsize; ++d, ++l) {
181: PetscInt offset = -(goff + 1) + d, r;
183: PetscCall(PetscFindInt(offset, size + 1, ranges, &r));
184: if (r < 0) r = -(r + 2);
185: PetscCheck(!(r < 0) && !(r >= size), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " mapped to invalid process %" PetscInt_FMT " (%" PetscInt_FMT ", %" PetscInt_FMT ")", p, r, gdof, goff);
186: remote[l].rank = r;
187: remote[l].index = offset - ranges[r];
188: }
189: } else {
190: for (d = 0; d < gsize; ++d, ++l) {
191: remote[l].rank = rank;
192: remote[l].index = goff + d - ranges[rank];
193: }
194: }
195: }
196: PetscCheck(l == nleaves, comm, PETSC_ERR_PLIB, "Iteration error, l %" PetscInt_FMT " != nleaves %" PetscInt_FMT, l, nleaves);
197: PetscCall(PetscLayoutDestroy(&layout));
198: PetscCall(PetscSFSetGraph(sf, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER));
199: PetscFunctionReturn(PETSC_SUCCESS);
200: }
202: /*@C
203: PetscSFDistributeSection - Create a new `PetscSection` reorganized, moving from the root to the leaves of the `PetscSF`
205: Collective
207: Input Parameters:
208: + sf - The `PetscSF`
209: - rootSection - Section defined on root space
211: Output Parameters:
212: + remoteOffsets - root offsets in leaf storage, or NULL
213: - leafSection - Section defined on the leaf space
215: Level: advanced
217: Fortran Notes:
218: In Fortran, use PetscSFDistributeSectionF90()
220: .seealso: `PetscSF`, `PetscSFCreate()`
221: @*/
222: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
223: {
224: PetscSF embedSF;
225: const PetscInt *indices;
226: IS selected;
227: PetscInt numFields, nroots, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, f, c;
228: PetscBool *sub, hasc;
230: PetscFunctionBegin;
231: PetscCall(PetscLogEventBegin(PETSCSF_DistSect, sf, 0, 0, 0));
232: PetscCall(PetscSectionGetNumFields(rootSection, &numFields));
233: if (numFields) {
234: IS perm;
236: /* PetscSectionSetNumFields() calls PetscSectionReset(), which destroys
237: leafSection->perm. To keep this permutation set by the user, we grab
238: the reference before calling PetscSectionSetNumFields() and set it
239: back after. */
240: PetscCall(PetscSectionGetPermutation(leafSection, &perm));
241: PetscCall(PetscObjectReference((PetscObject)perm));
242: PetscCall(PetscSectionSetNumFields(leafSection, numFields));
243: PetscCall(PetscSectionSetPermutation(leafSection, perm));
244: PetscCall(ISDestroy(&perm));
245: }
246: PetscCall(PetscMalloc1(numFields + 2, &sub));
247: sub[1] = rootSection->bc ? PETSC_TRUE : PETSC_FALSE;
248: for (f = 0; f < numFields; ++f) {
249: PetscSectionSym sym, dsym = NULL;
250: const char *name = NULL;
251: PetscInt numComp = 0;
253: sub[2 + f] = rootSection->field[f]->bc ? PETSC_TRUE : PETSC_FALSE;
254: PetscCall(PetscSectionGetFieldComponents(rootSection, f, &numComp));
255: PetscCall(PetscSectionGetFieldName(rootSection, f, &name));
256: PetscCall(PetscSectionGetFieldSym(rootSection, f, &sym));
257: if (sym) PetscCall(PetscSectionSymDistribute(sym, sf, &dsym));
258: PetscCall(PetscSectionSetFieldComponents(leafSection, f, numComp));
259: PetscCall(PetscSectionSetFieldName(leafSection, f, name));
260: PetscCall(PetscSectionSetFieldSym(leafSection, f, dsym));
261: PetscCall(PetscSectionSymDestroy(&dsym));
262: for (c = 0; c < rootSection->numFieldComponents[f]; ++c) {
263: PetscCall(PetscSectionGetComponentName(rootSection, f, c, &name));
264: PetscCall(PetscSectionSetComponentName(leafSection, f, c, name));
265: }
266: }
267: PetscCall(PetscSectionGetChart(rootSection, &rpStart, &rpEnd));
268: PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL));
269: rpEnd = PetscMin(rpEnd, nroots);
270: rpEnd = PetscMax(rpStart, rpEnd);
271: /* see if we can avoid creating the embedded SF, since it can cost more than an allreduce */
272: sub[0] = (PetscBool)(nroots != rpEnd - rpStart);
273: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, sub, 2 + numFields, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)sf)));
274: if (sub[0]) {
275: PetscCall(ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected));
276: PetscCall(ISGetIndices(selected, &indices));
277: PetscCall(PetscSFCreateEmbeddedRootSF(sf, rpEnd - rpStart, indices, &embedSF));
278: PetscCall(ISRestoreIndices(selected, &indices));
279: PetscCall(ISDestroy(&selected));
280: } else {
281: PetscCall(PetscObjectReference((PetscObject)sf));
282: embedSF = sf;
283: }
284: PetscCall(PetscSFGetLeafRange(embedSF, &lpStart, &lpEnd));
285: lpEnd++;
287: PetscCall(PetscSectionSetChart(leafSection, lpStart, lpEnd));
289: /* Constrained dof section */
290: hasc = sub[1];
291: for (f = 0; f < numFields; ++f) hasc = (PetscBool)(hasc || sub[2 + f]);
293: /* Could fuse these at the cost of copies and extra allocation */
294: PetscCall(PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart], MPI_REPLACE));
295: PetscCall(PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart], MPI_REPLACE));
296: if (sub[1]) {
297: PetscCall(PetscSectionCheckConstraints_Private(rootSection));
298: PetscCall(PetscSectionCheckConstraints_Private(leafSection));
299: PetscCall(PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart], MPI_REPLACE));
300: PetscCall(PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart], MPI_REPLACE));
301: }
302: for (f = 0; f < numFields; ++f) {
303: PetscCall(PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart], MPI_REPLACE));
304: PetscCall(PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart], MPI_REPLACE));
305: if (sub[2 + f]) {
306: PetscCall(PetscSectionCheckConstraints_Private(rootSection->field[f]));
307: PetscCall(PetscSectionCheckConstraints_Private(leafSection->field[f]));
308: PetscCall(PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart], MPI_REPLACE));
309: PetscCall(PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart], MPI_REPLACE));
310: }
311: }
312: if (remoteOffsets) {
313: PetscCall(PetscMalloc1(lpEnd - lpStart, remoteOffsets));
314: PetscCall(PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart], MPI_REPLACE));
315: PetscCall(PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart], MPI_REPLACE));
316: }
317: PetscCall(PetscSectionInvalidateMaxDof_Internal(leafSection));
318: PetscCall(PetscSectionSetUp(leafSection));
319: if (hasc) { /* need to communicate bcIndices */
320: PetscSF bcSF;
321: PetscInt *rOffBc;
323: PetscCall(PetscMalloc1(lpEnd - lpStart, &rOffBc));
324: if (sub[1]) {
325: PetscCall(PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart], MPI_REPLACE));
326: PetscCall(PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart], MPI_REPLACE));
327: PetscCall(PetscSFCreateSectionSF(embedSF, rootSection->bc, rOffBc, leafSection->bc, &bcSF));
328: PetscCall(PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices, MPI_REPLACE));
329: PetscCall(PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices, MPI_REPLACE));
330: PetscCall(PetscSFDestroy(&bcSF));
331: }
332: for (f = 0; f < numFields; ++f) {
333: if (sub[2 + f]) {
334: PetscCall(PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart], MPI_REPLACE));
335: PetscCall(PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart], MPI_REPLACE));
336: PetscCall(PetscSFCreateSectionSF(embedSF, rootSection->field[f]->bc, rOffBc, leafSection->field[f]->bc, &bcSF));
337: PetscCall(PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices, MPI_REPLACE));
338: PetscCall(PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices, MPI_REPLACE));
339: PetscCall(PetscSFDestroy(&bcSF));
340: }
341: }
342: PetscCall(PetscFree(rOffBc));
343: }
344: PetscCall(PetscSFDestroy(&embedSF));
345: PetscCall(PetscFree(sub));
346: PetscCall(PetscLogEventEnd(PETSCSF_DistSect, sf, 0, 0, 0));
347: PetscFunctionReturn(PETSC_SUCCESS);
348: }
350: /*@C
351: PetscSFCreateRemoteOffsets - Create offsets for point data on remote processes
353: Collective
355: Input Parameters:
356: + sf - The `PetscSF`
357: . rootSection - Data layout of remote points for outgoing data (this is layout for SF roots)
358: - leafSection - Data layout of local points for incoming data (this is layout for SF leaves)
360: Output Parameter:
361: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
363: Level: developer
365: Fortran Notes:
366: In Fortran, use PetscSFCreateRemoteOffsetsF90()
368: .seealso: `PetscSF`, `PetscSFCreate()`
369: @*/
370: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
371: {
372: PetscSF embedSF;
373: const PetscInt *indices;
374: IS selected;
375: PetscInt numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
377: PetscFunctionBegin;
378: *remoteOffsets = NULL;
379: PetscCall(PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL));
380: if (numRoots < 0) PetscFunctionReturn(PETSC_SUCCESS);
381: PetscCall(PetscLogEventBegin(PETSCSF_RemoteOff, sf, 0, 0, 0));
382: PetscCall(PetscSectionGetChart(rootSection, &rpStart, &rpEnd));
383: PetscCall(PetscSectionGetChart(leafSection, &lpStart, &lpEnd));
384: PetscCall(ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected));
385: PetscCall(ISGetIndices(selected, &indices));
386: PetscCall(PetscSFCreateEmbeddedRootSF(sf, rpEnd - rpStart, indices, &embedSF));
387: PetscCall(ISRestoreIndices(selected, &indices));
388: PetscCall(ISDestroy(&selected));
389: PetscCall(PetscCalloc1(lpEnd - lpStart, remoteOffsets));
390: PetscCall(PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart], MPI_REPLACE));
391: PetscCall(PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart], MPI_REPLACE));
392: PetscCall(PetscSFDestroy(&embedSF));
393: PetscCall(PetscLogEventEnd(PETSCSF_RemoteOff, sf, 0, 0, 0));
394: PetscFunctionReturn(PETSC_SUCCESS);
395: }
397: /*@C
398: PetscSFCreateSectionSF - Create an expanded `PetscSF` of dofs, assuming the input `PetscSF` relates points
400: Collective
402: Input Parameters:
403: + sf - The `PetscSF`
404: . rootSection - Data layout of remote points for outgoing data (this is usually the serial section)
405: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
406: - leafSection - Data layout of local points for incoming data (this is the distributed section)
408: Output Parameter:
409: . sectionSF - The new `PetscSF`
411: Level: advanced
413: Notes:
414: Either rootSection or remoteOffsets can be specified
416: Fortran Notes:
417: In Fortran, use PetscSFCreateSectionSFF90()
419: .seealso: `PetscSF`, `PetscSFCreate()`
420: @*/
421: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
422: {
423: MPI_Comm comm;
424: const PetscInt *localPoints;
425: const PetscSFNode *remotePoints;
426: PetscInt lpStart, lpEnd;
427: PetscInt numRoots, numSectionRoots, numPoints, numIndices = 0;
428: PetscInt *localIndices;
429: PetscSFNode *remoteIndices;
430: PetscInt i, ind;
432: PetscFunctionBegin;
434: PetscAssertPointer(rootSection, 2);
435: /* Cannot check PetscAssertPointer(remoteOffsets,3) because it can be NULL if sf does not reference any points in leafSection */
436: PetscAssertPointer(leafSection, 4);
437: PetscAssertPointer(sectionSF, 5);
438: PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
439: PetscCall(PetscSFCreate(comm, sectionSF));
440: PetscCall(PetscSectionGetChart(leafSection, &lpStart, &lpEnd));
441: PetscCall(PetscSectionGetStorageSize(rootSection, &numSectionRoots));
442: PetscCall(PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints));
443: if (numRoots < 0) PetscFunctionReturn(PETSC_SUCCESS);
444: PetscCall(PetscLogEventBegin(PETSCSF_SectSF, sf, 0, 0, 0));
445: for (i = 0; i < numPoints; ++i) {
446: PetscInt localPoint = localPoints ? localPoints[i] : i;
447: PetscInt dof;
449: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
450: PetscCall(PetscSectionGetDof(leafSection, localPoint, &dof));
451: numIndices += dof < 0 ? 0 : dof;
452: }
453: }
454: PetscCall(PetscMalloc1(numIndices, &localIndices));
455: PetscCall(PetscMalloc1(numIndices, &remoteIndices));
456: /* Create new index graph */
457: for (i = 0, ind = 0; i < numPoints; ++i) {
458: PetscInt localPoint = localPoints ? localPoints[i] : i;
459: PetscInt rank = remotePoints[i].rank;
461: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
462: PetscInt remoteOffset = remoteOffsets[localPoint - lpStart];
463: PetscInt loff, dof, d;
465: PetscCall(PetscSectionGetOffset(leafSection, localPoint, &loff));
466: PetscCall(PetscSectionGetDof(leafSection, localPoint, &dof));
467: for (d = 0; d < dof; ++d, ++ind) {
468: localIndices[ind] = loff + d;
469: remoteIndices[ind].rank = rank;
470: remoteIndices[ind].index = remoteOffset + d;
471: }
472: }
473: }
474: PetscCheck(numIndices == ind, comm, PETSC_ERR_PLIB, "Inconsistency in indices, %" PetscInt_FMT " should be %" PetscInt_FMT, ind, numIndices);
475: PetscCall(PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER));
476: PetscCall(PetscSFSetUp(*sectionSF));
477: PetscCall(PetscLogEventEnd(PETSCSF_SectSF, sf, 0, 0, 0));
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: /*@C
482: PetscSFCreateFromLayouts - Creates a parallel star forest mapping two `PetscLayout` objects
484: Collective
486: Input Parameters:
487: + rmap - `PetscLayout` defining the global root space
488: - lmap - `PetscLayout` defining the global leaf space
490: Output Parameter:
491: . sf - The parallel star forest
493: Level: intermediate
495: .seealso: `PetscSF`, `PetscSFCreate()`, `PetscLayoutCreate()`, `PetscSFSetGraphLayout()`
496: @*/
497: PetscErrorCode PetscSFCreateFromLayouts(PetscLayout rmap, PetscLayout lmap, PetscSF *sf)
498: {
499: PetscInt i, nroots, nleaves = 0;
500: PetscInt rN, lst, len;
501: PetscMPIInt owner = -1;
502: PetscSFNode *remote;
503: MPI_Comm rcomm = rmap->comm;
504: MPI_Comm lcomm = lmap->comm;
505: PetscMPIInt flg;
507: PetscFunctionBegin;
508: PetscAssertPointer(sf, 3);
509: PetscCheck(rmap->setupcalled, rcomm, PETSC_ERR_ARG_WRONGSTATE, "Root layout not setup");
510: PetscCheck(lmap->setupcalled, lcomm, PETSC_ERR_ARG_WRONGSTATE, "Leaf layout not setup");
511: PetscCallMPI(MPI_Comm_compare(rcomm, lcomm, &flg));
512: PetscCheck(flg == MPI_CONGRUENT || flg == MPI_IDENT, rcomm, PETSC_ERR_SUP, "cannot map two layouts with non-matching communicators");
513: PetscCall(PetscSFCreate(rcomm, sf));
514: PetscCall(PetscLayoutGetLocalSize(rmap, &nroots));
515: PetscCall(PetscLayoutGetSize(rmap, &rN));
516: PetscCall(PetscLayoutGetRange(lmap, &lst, &len));
517: PetscCall(PetscMalloc1(len - lst, &remote));
518: for (i = lst; i < len && i < rN; i++) {
519: if (owner < -1 || i >= rmap->range[owner + 1]) PetscCall(PetscLayoutFindOwner(rmap, i, &owner));
520: remote[nleaves].rank = owner;
521: remote[nleaves].index = i - rmap->range[owner];
522: nleaves++;
523: }
524: PetscCall(PetscSFSetGraph(*sf, nroots, nleaves, NULL, PETSC_OWN_POINTER, remote, PETSC_COPY_VALUES));
525: PetscCall(PetscFree(remote));
526: PetscFunctionReturn(PETSC_SUCCESS);
527: }
529: /* TODO: handle nooffprocentries like MatZeroRowsMapLocal_Private, since this code is the same */
530: PetscErrorCode PetscLayoutMapLocal(PetscLayout map, PetscInt N, const PetscInt idxs[], PetscInt *on, PetscInt **oidxs, PetscInt **ogidxs)
531: {
532: PetscInt *owners = map->range;
533: PetscInt n = map->n;
534: PetscSF sf;
535: PetscInt *lidxs, *work = NULL;
536: PetscSFNode *ridxs;
537: PetscMPIInt rank, p = 0;
538: PetscInt r, len = 0;
540: PetscFunctionBegin;
541: if (on) *on = 0; /* squelch -Wmaybe-uninitialized */
542: /* Create SF where leaves are input idxs and roots are owned idxs */
543: PetscCallMPI(MPI_Comm_rank(map->comm, &rank));
544: PetscCall(PetscMalloc1(n, &lidxs));
545: for (r = 0; r < n; ++r) lidxs[r] = -1;
546: PetscCall(PetscMalloc1(N, &ridxs));
547: for (r = 0; r < N; ++r) {
548: const PetscInt idx = idxs[r];
549: PetscCheck(idx >= 0 && idx < map->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Index %" PetscInt_FMT " out of range [0,%" PetscInt_FMT ")", idx, map->N);
550: if (idx < owners[p] || owners[p + 1] <= idx) { /* short-circuit the search if the last p owns this idx too */
551: PetscCall(PetscLayoutFindOwner(map, idx, &p));
552: }
553: ridxs[r].rank = p;
554: ridxs[r].index = idxs[r] - owners[p];
555: }
556: PetscCall(PetscSFCreate(map->comm, &sf));
557: PetscCall(PetscSFSetGraph(sf, n, N, NULL, PETSC_OWN_POINTER, ridxs, PETSC_OWN_POINTER));
558: PetscCall(PetscSFReduceBegin(sf, MPIU_INT, (PetscInt *)idxs, lidxs, MPI_LOR));
559: PetscCall(PetscSFReduceEnd(sf, MPIU_INT, (PetscInt *)idxs, lidxs, MPI_LOR));
560: if (ogidxs) { /* communicate global idxs */
561: PetscInt cum = 0, start, *work2;
563: PetscCall(PetscMalloc1(n, &work));
564: PetscCall(PetscCalloc1(N, &work2));
565: for (r = 0; r < N; ++r)
566: if (idxs[r] >= 0) cum++;
567: PetscCallMPI(MPI_Scan(&cum, &start, 1, MPIU_INT, MPI_SUM, map->comm));
568: start -= cum;
569: cum = 0;
570: for (r = 0; r < N; ++r)
571: if (idxs[r] >= 0) work2[r] = start + cum++;
572: PetscCall(PetscSFReduceBegin(sf, MPIU_INT, work2, work, MPI_REPLACE));
573: PetscCall(PetscSFReduceEnd(sf, MPIU_INT, work2, work, MPI_REPLACE));
574: PetscCall(PetscFree(work2));
575: }
576: PetscCall(PetscSFDestroy(&sf));
577: /* Compress and put in indices */
578: for (r = 0; r < n; ++r)
579: if (lidxs[r] >= 0) {
580: if (work) work[len] = work[r];
581: lidxs[len++] = r;
582: }
583: if (on) *on = len;
584: if (oidxs) *oidxs = lidxs;
585: if (ogidxs) *ogidxs = work;
586: PetscFunctionReturn(PETSC_SUCCESS);
587: }
589: /*@
590: PetscSFCreateByMatchingIndices - Create `PetscSF` by matching root and leaf indices
592: Collective
594: Input Parameters:
595: + layout - `PetscLayout` defining the global index space and the rank that brokers each index
596: . numRootIndices - size of rootIndices
597: . rootIndices - `PetscInt` array of global indices of which this process requests ownership
598: . rootLocalIndices - root local index permutation (NULL if no permutation)
599: . rootLocalOffset - offset to be added to root local indices
600: . numLeafIndices - size of leafIndices
601: . leafIndices - `PetscInt` array of global indices with which this process requires data associated
602: . leafLocalIndices - leaf local index permutation (NULL if no permutation)
603: - leafLocalOffset - offset to be added to leaf local indices
605: Output Parameters:
606: + sfA - star forest representing the communication pattern from the layout space to the leaf space (NULL if not needed)
607: - sf - star forest representing the communication pattern from the root space to the leaf space
609: Level: advanced
611: Example 1:
612: .vb
613: rank : 0 1 2
614: rootIndices : [1 0 2] [3] [3]
615: rootLocalOffset : 100 200 300
616: layout : [0 1] [2] [3]
617: leafIndices : [0] [2] [0 3]
618: leafLocalOffset : 400 500 600
620: would build the following PetscSF
622: [0] 400 <- (0,101)
623: [1] 500 <- (0,102)
624: [2] 600 <- (0,101)
625: [2] 601 <- (2,300)
626: .ve
628: Example 2:
629: .vb
630: rank : 0 1 2
631: rootIndices : [1 0 2] [3] [3]
632: rootLocalOffset : 100 200 300
633: layout : [0 1] [2] [3]
634: leafIndices : rootIndices rootIndices rootIndices
635: leafLocalOffset : rootLocalOffset rootLocalOffset rootLocalOffset
637: would build the following PetscSF
639: [1] 200 <- (2,300)
640: .ve
642: Example 3:
643: .vb
644: No process requests ownership of global index 1, but no process needs it.
646: rank : 0 1 2
647: numRootIndices : 2 1 1
648: rootIndices : [0 2] [3] [3]
649: rootLocalOffset : 100 200 300
650: layout : [0 1] [2] [3]
651: numLeafIndices : 1 1 2
652: leafIndices : [0] [2] [0 3]
653: leafLocalOffset : 400 500 600
655: would build the following PetscSF
657: [0] 400 <- (0,100)
658: [1] 500 <- (0,101)
659: [2] 600 <- (0,100)
660: [2] 601 <- (2,300)
661: .ve
663: Notes:
664: The layout parameter represents any partitioning of [0, N), where N is the total number of global indices, and its
665: local size can be set to `PETSC_DECIDE`.
667: If a global index x lies in the partition owned by process i, each process whose rootIndices contains x requests
668: ownership of x and sends its own rank and the local index of x to process i.
669: If multiple processes request ownership of x, the one with the highest rank is to own x.
670: Process i then broadcasts the ownership information, so that each process whose leafIndices contains x knows the
671: ownership information of x.
672: The output sf is constructed by associating each leaf point to a root point in this way.
674: Suppose there is point data ordered according to the global indices and partitioned according to the given layout.
675: The optional output `PetscSF` object sfA can be used to push such data to leaf points.
677: All indices in rootIndices and leafIndices must lie in the layout range. The union (over all processes) of rootIndices
678: must cover that of leafIndices, but need not cover the entire layout.
680: If (leafIndices, leafLocalIndices, leafLocalOffset) == (rootIndices, rootLocalIndices, rootLocalOffset), the output
681: star forest is almost identity, so will only include non-trivial part of the map.
683: Developer Notes:
684: Current approach of a process of the highest rank gaining the ownership may cause load imbalance; consider using
685: hash(rank, root_local_index) as the bid for the ownership determination.
687: .seealso: `PetscSF`, `PetscSFCreate()`
688: @*/
689: 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)
690: {
691: MPI_Comm comm = layout->comm;
692: PetscMPIInt size, rank;
693: PetscSF sf1;
694: PetscSFNode *owners, *buffer, *iremote;
695: PetscInt *ilocal, nleaves, N, n, i;
696: #if defined(PETSC_USE_DEBUG)
697: PetscInt N1;
698: #endif
699: PetscBool flag;
701: PetscFunctionBegin;
702: if (rootIndices) PetscAssertPointer(rootIndices, 3);
703: if (rootLocalIndices) PetscAssertPointer(rootLocalIndices, 4);
704: if (leafIndices) PetscAssertPointer(leafIndices, 7);
705: if (leafLocalIndices) PetscAssertPointer(leafLocalIndices, 8);
706: if (sfA) PetscAssertPointer(sfA, 10);
707: PetscAssertPointer(sf, 11);
708: PetscCheck(numRootIndices >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "numRootIndices (%" PetscInt_FMT ") must be non-negative", numRootIndices);
709: PetscCheck(numLeafIndices >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "numLeafIndices (%" PetscInt_FMT ") must be non-negative", numLeafIndices);
710: PetscCallMPI(MPI_Comm_size(comm, &size));
711: PetscCallMPI(MPI_Comm_rank(comm, &rank));
712: PetscCall(PetscLayoutSetUp(layout));
713: PetscCall(PetscLayoutGetSize(layout, &N));
714: PetscCall(PetscLayoutGetLocalSize(layout, &n));
715: flag = (PetscBool)(leafIndices == rootIndices);
716: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &flag, 1, MPIU_BOOL, MPI_LAND, comm));
717: PetscCheck(!flag || numLeafIndices == numRootIndices, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "leafIndices == rootIndices, but numLeafIndices (%" PetscInt_FMT ") != numRootIndices(%" PetscInt_FMT ")", numLeafIndices, numRootIndices);
718: #if defined(PETSC_USE_DEBUG)
719: N1 = PETSC_MIN_INT;
720: for (i = 0; i < numRootIndices; i++)
721: if (rootIndices[i] > N1) N1 = rootIndices[i];
722: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &N1, 1, MPIU_INT, MPI_MAX, comm));
723: PetscCheck(N1 < N, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Max. root index (%" PetscInt_FMT ") out of layout range [0,%" PetscInt_FMT ")", N1, N);
724: if (!flag) {
725: N1 = PETSC_MIN_INT;
726: for (i = 0; i < numLeafIndices; i++)
727: if (leafIndices[i] > N1) N1 = leafIndices[i];
728: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &N1, 1, MPIU_INT, MPI_MAX, comm));
729: PetscCheck(N1 < N, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Max. leaf index (%" PetscInt_FMT ") out of layout range [0,%" PetscInt_FMT ")", N1, N);
730: }
731: #endif
732: /* Reduce: owners -> buffer */
733: PetscCall(PetscMalloc1(n, &buffer));
734: PetscCall(PetscSFCreate(comm, &sf1));
735: PetscCall(PetscSFSetFromOptions(sf1));
736: PetscCall(PetscSFSetGraphLayout(sf1, layout, numRootIndices, NULL, PETSC_OWN_POINTER, rootIndices));
737: PetscCall(PetscMalloc1(numRootIndices, &owners));
738: for (i = 0; i < numRootIndices; ++i) {
739: owners[i].rank = rank;
740: owners[i].index = rootLocalOffset + (rootLocalIndices ? rootLocalIndices[i] : i);
741: }
742: for (i = 0; i < n; ++i) {
743: buffer[i].index = -1;
744: buffer[i].rank = -1;
745: }
746: PetscCall(PetscSFReduceBegin(sf1, MPIU_2INT, owners, buffer, MPI_MAXLOC));
747: PetscCall(PetscSFReduceEnd(sf1, MPIU_2INT, owners, buffer, MPI_MAXLOC));
748: /* Bcast: buffer -> owners */
749: if (!flag) {
750: /* leafIndices is different from rootIndices */
751: PetscCall(PetscFree(owners));
752: PetscCall(PetscSFSetGraphLayout(sf1, layout, numLeafIndices, NULL, PETSC_OWN_POINTER, leafIndices));
753: PetscCall(PetscMalloc1(numLeafIndices, &owners));
754: }
755: PetscCall(PetscSFBcastBegin(sf1, MPIU_2INT, buffer, owners, MPI_REPLACE));
756: PetscCall(PetscSFBcastEnd(sf1, MPIU_2INT, buffer, owners, MPI_REPLACE));
757: 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]);
758: PetscCall(PetscFree(buffer));
759: if (sfA) {
760: *sfA = sf1;
761: } else PetscCall(PetscSFDestroy(&sf1));
762: /* Create sf */
763: if (flag && rootLocalIndices == leafLocalIndices && leafLocalOffset == rootLocalOffset) {
764: /* leaf space == root space */
765: for (i = 0, nleaves = 0; i < numLeafIndices; ++i)
766: if (owners[i].rank != rank) ++nleaves;
767: PetscCall(PetscMalloc1(nleaves, &ilocal));
768: PetscCall(PetscMalloc1(nleaves, &iremote));
769: for (i = 0, nleaves = 0; i < numLeafIndices; ++i) {
770: if (owners[i].rank != rank) {
771: ilocal[nleaves] = leafLocalOffset + i;
772: iremote[nleaves].rank = owners[i].rank;
773: iremote[nleaves].index = owners[i].index;
774: ++nleaves;
775: }
776: }
777: PetscCall(PetscFree(owners));
778: } else {
779: nleaves = numLeafIndices;
780: PetscCall(PetscMalloc1(nleaves, &ilocal));
781: for (i = 0; i < nleaves; ++i) ilocal[i] = leafLocalOffset + (leafLocalIndices ? leafLocalIndices[i] : i);
782: iremote = owners;
783: }
784: PetscCall(PetscSFCreate(comm, sf));
785: PetscCall(PetscSFSetFromOptions(*sf));
786: PetscCall(PetscSFSetGraph(*sf, rootLocalOffset + numRootIndices, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
787: PetscFunctionReturn(PETSC_SUCCESS);
788: }
790: /*@
791: PetscSFMerge - append/merge indices of `sfb` into `sfa`, with preference for `sfb`
793: Collective
795: Input Parameters:
796: + sfa - default `PetscSF`
797: - sfb - additional edges to add/replace edges in sfa
799: Output Parameter:
800: . merged - new `PetscSF` with combined edges
802: Level: intermediate
804: .seealso: `PetscSFCompose()`
805: @*/
806: PetscErrorCode PetscSFMerge(PetscSF sfa, PetscSF sfb, PetscSF *merged)
807: {
808: PetscInt maxleaf;
810: PetscFunctionBegin;
813: PetscCheckSameComm(sfa, 1, sfb, 2);
814: PetscAssertPointer(merged, 3);
815: {
816: PetscInt aleaf, bleaf;
817: PetscCall(PetscSFGetLeafRange(sfa, NULL, &aleaf));
818: PetscCall(PetscSFGetLeafRange(sfb, NULL, &bleaf));
819: maxleaf = PetscMax(aleaf, bleaf) + 1; // One more than the last index
820: }
821: PetscInt *clocal, aroots, aleaves, broots, bleaves;
822: PetscSFNode *cremote;
823: const PetscInt *alocal, *blocal;
824: const PetscSFNode *aremote, *bremote;
825: PetscCall(PetscMalloc2(maxleaf, &clocal, maxleaf, &cremote));
826: for (PetscInt i = 0; i < maxleaf; i++) clocal[i] = -1;
827: PetscCall(PetscSFGetGraph(sfa, &aroots, &aleaves, &alocal, &aremote));
828: PetscCall(PetscSFGetGraph(sfb, &broots, &bleaves, &blocal, &bremote));
829: PetscCheck(aroots == broots, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Both sfa and sfb must have the same root space");
830: for (PetscInt i = 0; i < aleaves; i++) {
831: PetscInt a = alocal ? alocal[i] : i;
832: clocal[a] = a;
833: cremote[a] = aremote[i];
834: }
835: for (PetscInt i = 0; i < bleaves; i++) {
836: PetscInt b = blocal ? blocal[i] : i;
837: clocal[b] = b;
838: cremote[b] = bremote[i];
839: }
840: PetscInt nleaves = 0;
841: for (PetscInt i = 0; i < maxleaf; i++) {
842: if (clocal[i] < 0) continue;
843: clocal[nleaves] = clocal[i];
844: cremote[nleaves] = cremote[i];
845: nleaves++;
846: }
847: PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)sfa), merged));
848: PetscCall(PetscSFSetGraph(*merged, aroots, nleaves, clocal, PETSC_COPY_VALUES, cremote, PETSC_COPY_VALUES));
849: PetscCall(PetscFree2(clocal, cremote));
850: PetscFunctionReturn(PETSC_SUCCESS);
851: }