Actual source code: sfutils.c
1: #include <petsc/private/sfimpl.h>
2: #include <petsc/private/sectionimpl.h>
4: /*@C
5: PetscSFSetGraphLayout - Set a parallel star forest via global indices and a PetscLayout
7: Collective
9: Input Parameters:
10: + sf - star forest
11: . layout - PetscLayout defining the global space for roots
12: . nleaves - number of leaf vertices on the current process, each of these references a root on any process
13: . ilocal - locations of leaves in leafdata buffers, pass NULL for contiguous storage
14: . localmode - copy mode for ilocal
15: - iremote - remote locations (in global indices) of root vertices for each leaf on the current process
17: Level: intermediate
19: Developers Note: Local indices which are the identity permutation in the range [0,nleaves) are discarded as they
20: encode contiguous storage. In such case, if localmode is PETSC_OWN_POINTER, the memory is deallocated as it is not
21: needed
23: .seealso: PetscSFCreate(), PetscSFView(), PetscSFSetGraph(), PetscSFGetGraph()
24: @*/
25: PetscErrorCode PetscSFSetGraphLayout(PetscSF sf,PetscLayout layout,PetscInt nleaves,const PetscInt *ilocal,PetscCopyMode localmode,const PetscInt *iremote)
26: {
28: const PetscInt *range;
29: PetscInt i, nroots, ls = -1, ln = -1;
30: PetscMPIInt lr = -1;
31: PetscSFNode *remote;
34: PetscLayoutGetLocalSize(layout,&nroots);
35: PetscLayoutGetRanges(layout,&range);
36: PetscMalloc1(nleaves,&remote);
37: if (nleaves) { ls = iremote[0] + 1; }
38: for (i=0; i<nleaves; i++) {
39: const PetscInt idx = iremote[i] - ls;
40: if (idx < 0 || idx >= ln) { /* short-circuit the search */
41: PetscLayoutFindOwnerIndex(layout,iremote[i],&lr,&remote[i].index);
42: remote[i].rank = lr;
43: ls = range[lr];
44: ln = range[lr+1] - ls;
45: } else {
46: remote[i].rank = lr;
47: remote[i].index = idx;
48: }
49: }
50: PetscSFSetGraph(sf,nroots,nleaves,ilocal,localmode,remote,PETSC_OWN_POINTER);
51: return(0);
52: }
54: /*@
55: PetscSFSetGraphSection - Sets the PetscSF graph encoding the parallel dof overlap based upon the PetscSections
56: describing the data layout.
58: Input Parameters:
59: + sf - The SF
60: . localSection - PetscSection describing the local data layout
61: - globalSection - PetscSection describing the global data layout
63: Level: developer
65: .seealso: PetscSFSetGraph(), PetscSFSetGraphLayout()
66: @*/
67: PetscErrorCode PetscSFSetGraphSection(PetscSF sf, PetscSection localSection, PetscSection globalSection)
68: {
69: MPI_Comm comm;
70: PetscLayout layout;
71: const PetscInt *ranges;
72: PetscInt *local;
73: PetscSFNode *remote;
74: PetscInt pStart, pEnd, p, nroots, nleaves = 0, l;
75: PetscMPIInt size, rank;
83: PetscObjectGetComm((PetscObject)sf,&comm);
84: MPI_Comm_size(comm, &size);
85: MPI_Comm_rank(comm, &rank);
86: PetscSectionGetChart(globalSection, &pStart, &pEnd);
87: PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
88: PetscLayoutCreate(comm, &layout);
89: PetscLayoutSetBlockSize(layout, 1);
90: PetscLayoutSetLocalSize(layout, nroots);
91: PetscLayoutSetUp(layout);
92: PetscLayoutGetRanges(layout, &ranges);
93: for (p = pStart; p < pEnd; ++p) {
94: PetscInt gdof, gcdof;
96: PetscSectionGetDof(globalSection, p, &gdof);
97: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
98: if (gcdof > (gdof < 0 ? -(gdof+1) : gdof)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %D has %D constraints > %D dof", p, gcdof, (gdof < 0 ? -(gdof+1) : gdof));
99: nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
100: }
101: PetscMalloc1(nleaves, &local);
102: PetscMalloc1(nleaves, &remote);
103: for (p = pStart, l = 0; p < pEnd; ++p) {
104: const PetscInt *cind;
105: PetscInt dof, cdof, off, gdof, gcdof, goff, gsize, d, c;
107: PetscSectionGetDof(localSection, p, &dof);
108: PetscSectionGetOffset(localSection, p, &off);
109: PetscSectionGetConstraintDof(localSection, p, &cdof);
110: PetscSectionGetConstraintIndices(localSection, p, &cind);
111: PetscSectionGetDof(globalSection, p, &gdof);
112: PetscSectionGetConstraintDof(globalSection, p, &gcdof);
113: PetscSectionGetOffset(globalSection, p, &goff);
114: if (!gdof) continue; /* Censored point */
115: gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
116: if (gsize != dof-cdof) {
117: if (gsize != dof) SETERRQ4(comm, PETSC_ERR_ARG_WRONG, "Global dof %D for point %D is neither the constrained size %D, nor the unconstrained %D", gsize, p, dof-cdof, dof);
118: cdof = 0; /* Ignore constraints */
119: }
120: for (d = 0, c = 0; d < dof; ++d) {
121: if ((c < cdof) && (cind[c] == d)) {++c; continue;}
122: local[l+d-c] = off+d;
123: }
124: if (d-c != gsize) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Point %D: Global dof %D != %D size - number of constraints", p, gsize, d-c);
125: if (gdof < 0) {
126: for (d = 0; d < gsize; ++d, ++l) {
127: PetscInt offset = -(goff+1) + d, r;
129: PetscFindInt(offset,size+1,ranges,&r);
130: if (r < 0) r = -(r+2);
131: if ((r < 0) || (r >= size)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %D mapped to invalid process %D (%D, %D)", p, r, gdof, goff);
132: remote[l].rank = r;
133: remote[l].index = offset - ranges[r];
134: }
135: } else {
136: for (d = 0; d < gsize; ++d, ++l) {
137: remote[l].rank = rank;
138: remote[l].index = goff+d - ranges[rank];
139: }
140: }
141: }
142: if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %D != nleaves %D", l, nleaves);
143: PetscLayoutDestroy(&layout);
144: PetscSFSetGraph(sf, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
145: return(0);
146: }
148: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
149: {
153: if (!s->bc) {
154: PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
155: PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
156: }
157: return(0);
158: }
160: /*@C
161: PetscSFDistributeSection - Create a new PetscSection reorganized, moving from the root to the leaves of the SF
163: Collective on sf
165: Input Parameters:
166: + sf - The SF
167: - rootSection - Section defined on root space
169: Output Parameters:
170: + remoteOffsets - root offsets in leaf storage, or NULL
171: - leafSection - Section defined on the leaf space
173: Level: advanced
175: .seealso: PetscSFCreate()
176: @*/
177: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
178: {
179: PetscSF embedSF;
180: const PetscInt *indices;
181: IS selected;
182: PetscInt numFields, nroots, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, f, c;
183: PetscBool *sub, hasc;
187: PetscLogEventBegin(PETSCSF_DistSect,sf,0,0,0);
188: PetscSectionGetNumFields(rootSection, &numFields);
189: if (numFields) {
190: IS perm;
192: /* PetscSectionSetNumFields() calls PetscSectionReset(), which destroys
193: leafSection->perm. To keep this permutation set by the user, we grab
194: the reference before calling PetscSectionSetNumFields() and set it
195: back after. */
196: PetscSectionGetPermutation(leafSection, &perm);
197: PetscObjectReference((PetscObject)perm);
198: PetscSectionSetNumFields(leafSection, numFields);
199: PetscSectionSetPermutation(leafSection, perm);
200: ISDestroy(&perm);
201: }
202: PetscMalloc1(numFields+2, &sub);
203: sub[1] = rootSection->bc ? PETSC_TRUE : PETSC_FALSE;
204: for (f = 0; f < numFields; ++f) {
205: PetscSectionSym sym;
206: const char *name = NULL;
207: PetscInt numComp = 0;
209: sub[2 + f] = rootSection->field[f]->bc ? PETSC_TRUE : PETSC_FALSE;
210: PetscSectionGetFieldComponents(rootSection, f, &numComp);
211: PetscSectionGetFieldName(rootSection, f, &name);
212: PetscSectionGetFieldSym(rootSection, f, &sym);
213: PetscSectionSetFieldComponents(leafSection, f, numComp);
214: PetscSectionSetFieldName(leafSection, f, name);
215: PetscSectionSetFieldSym(leafSection, f, sym);
216: for (c = 0; c < rootSection->numFieldComponents[f]; ++c) {
217: PetscSectionGetComponentName(rootSection, f, c, &name);
218: PetscSectionSetComponentName(leafSection, f, c, name);
219: }
220: }
221: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
222: PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);
223: rpEnd = PetscMin(rpEnd,nroots);
224: rpEnd = PetscMax(rpStart,rpEnd);
225: /* see if we can avoid creating the embedded SF, since it can cost more than an allreduce */
226: sub[0] = (PetscBool)(nroots != rpEnd - rpStart);
227: MPIU_Allreduce(MPI_IN_PLACE, sub, 2+numFields, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)sf));
228: if (sub[0]) {
229: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
230: ISGetIndices(selected, &indices);
231: PetscSFCreateEmbeddedRootSF(sf, rpEnd - rpStart, indices, &embedSF);
232: ISRestoreIndices(selected, &indices);
233: ISDestroy(&selected);
234: } else {
235: PetscObjectReference((PetscObject)sf);
236: embedSF = sf;
237: }
238: PetscSFGetLeafRange(embedSF, &lpStart, &lpEnd);
239: lpEnd++;
241: PetscSectionSetChart(leafSection, lpStart, lpEnd);
243: /* Constrained dof section */
244: hasc = sub[1];
245: for (f = 0; f < numFields; ++f) hasc = (PetscBool)(hasc || sub[2+f]);
247: /* Could fuse these at the cost of copies and extra allocation */
248: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart],MPI_REPLACE);
249: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart],MPI_REPLACE);
250: if (sub[1]) {
251: PetscSectionCheckConstraints_Static(rootSection);
252: PetscSectionCheckConstraints_Static(leafSection);
253: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart],MPI_REPLACE);
254: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart],MPI_REPLACE);
255: }
256: for (f = 0; f < numFields; ++f) {
257: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart],MPI_REPLACE);
258: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart],MPI_REPLACE);
259: if (sub[2+f]) {
260: PetscSectionCheckConstraints_Static(rootSection->field[f]);
261: PetscSectionCheckConstraints_Static(leafSection->field[f]);
262: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart],MPI_REPLACE);
263: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart],MPI_REPLACE);
264: }
265: }
266: if (remoteOffsets) {
267: PetscMalloc1(lpEnd - lpStart, remoteOffsets);
268: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart],MPI_REPLACE);
269: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart],MPI_REPLACE);
270: }
271: PetscSectionSetUp(leafSection);
272: if (hasc) { /* need to communicate bcIndices */
273: PetscSF bcSF;
274: PetscInt *rOffBc;
276: PetscMalloc1(lpEnd - lpStart, &rOffBc);
277: if (sub[1]) {
278: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart],MPI_REPLACE);
279: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart],MPI_REPLACE);
280: PetscSFCreateSectionSF(embedSF, rootSection->bc, rOffBc, leafSection->bc, &bcSF);
281: PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices,MPI_REPLACE);
282: PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices,MPI_REPLACE);
283: PetscSFDestroy(&bcSF);
284: }
285: for (f = 0; f < numFields; ++f) {
286: if (sub[2+f]) {
287: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart],MPI_REPLACE);
288: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart],MPI_REPLACE);
289: PetscSFCreateSectionSF(embedSF, rootSection->field[f]->bc, rOffBc, leafSection->field[f]->bc, &bcSF);
290: PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices,MPI_REPLACE);
291: PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices,MPI_REPLACE);
292: PetscSFDestroy(&bcSF);
293: }
294: }
295: PetscFree(rOffBc);
296: }
297: PetscSFDestroy(&embedSF);
298: PetscFree(sub);
299: PetscLogEventEnd(PETSCSF_DistSect,sf,0,0,0);
300: return(0);
301: }
303: /*@C
304: PetscSFCreateRemoteOffsets - Create offsets for point data on remote processes
306: Collective on sf
308: Input Parameters:
309: + sf - The SF
310: . rootSection - Data layout of remote points for outgoing data (this is layout for SF roots)
311: - leafSection - Data layout of local points for incoming data (this is layout for SF leaves)
313: Output Parameter:
314: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
316: Level: developer
318: .seealso: PetscSFCreate()
319: @*/
320: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
321: {
322: PetscSF embedSF;
323: const PetscInt *indices;
324: IS selected;
325: PetscInt numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
326: PetscErrorCode ierr;
329: *remoteOffsets = NULL;
330: PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
331: if (numRoots < 0) return(0);
332: PetscLogEventBegin(PETSCSF_RemoteOff,sf,0,0,0);
333: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
334: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
335: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
336: ISGetIndices(selected, &indices);
337: PetscSFCreateEmbeddedRootSF(sf, rpEnd - rpStart, indices, &embedSF);
338: ISRestoreIndices(selected, &indices);
339: ISDestroy(&selected);
340: PetscCalloc1(lpEnd - lpStart, remoteOffsets);
341: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart],MPI_REPLACE);
342: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart],MPI_REPLACE);
343: PetscSFDestroy(&embedSF);
344: PetscLogEventEnd(PETSCSF_RemoteOff,sf,0,0,0);
345: return(0);
346: }
348: /*@C
349: PetscSFCreateSectionSF - Create an expanded SF of dofs, assuming the input SF relates points
351: Collective on sf
353: Input Parameters:
354: + sf - The SF
355: . rootSection - Data layout of remote points for outgoing data (this is usually the serial section)
356: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
357: - leafSection - Data layout of local points for incoming data (this is the distributed section)
359: Output Parameters:
360: - sectionSF - The new SF
362: Note: Either rootSection or remoteOffsets can be specified
364: Level: advanced
366: .seealso: PetscSFCreate()
367: @*/
368: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
369: {
370: MPI_Comm comm;
371: const PetscInt *localPoints;
372: const PetscSFNode *remotePoints;
373: PetscInt lpStart, lpEnd;
374: PetscInt numRoots, numSectionRoots, numPoints, numIndices = 0;
375: PetscInt *localIndices;
376: PetscSFNode *remoteIndices;
377: PetscInt i, ind;
378: PetscErrorCode ierr;
386: PetscObjectGetComm((PetscObject)sf,&comm);
387: PetscSFCreate(comm, sectionSF);
388: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
389: PetscSectionGetStorageSize(rootSection, &numSectionRoots);
390: PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
391: if (numRoots < 0) return(0);
392: PetscLogEventBegin(PETSCSF_SectSF,sf,0,0,0);
393: for (i = 0; i < numPoints; ++i) {
394: PetscInt localPoint = localPoints ? localPoints[i] : i;
395: PetscInt dof;
397: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
398: PetscSectionGetDof(leafSection, localPoint, &dof);
399: numIndices += dof;
400: }
401: }
402: PetscMalloc1(numIndices, &localIndices);
403: PetscMalloc1(numIndices, &remoteIndices);
404: /* Create new index graph */
405: for (i = 0, ind = 0; i < numPoints; ++i) {
406: PetscInt localPoint = localPoints ? localPoints[i] : i;
407: PetscInt rank = remotePoints[i].rank;
409: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
410: PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
411: PetscInt loff, dof, d;
413: PetscSectionGetOffset(leafSection, localPoint, &loff);
414: PetscSectionGetDof(leafSection, localPoint, &dof);
415: for (d = 0; d < dof; ++d, ++ind) {
416: localIndices[ind] = loff+d;
417: remoteIndices[ind].rank = rank;
418: remoteIndices[ind].index = remoteOffset+d;
419: }
420: }
421: }
422: if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %D should be %D", ind, numIndices);
423: PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
424: PetscSFSetUp(*sectionSF);
425: PetscLogEventEnd(PETSCSF_SectSF,sf,0,0,0);
426: return(0);
427: }
429: /*@C
430: PetscSFCreateFromLayouts - Creates a parallel star forest mapping two PetscLayout objects
432: Collective
434: Input Parameters:
435: + rmap - PetscLayout defining the global root space
436: - lmap - PetscLayout defining the global leaf space
438: Output Parameter:
439: . sf - The parallel star forest
441: Level: intermediate
443: .seealso: PetscSFCreate(), PetscLayoutCreate(), PetscSFSetGraphLayout()
444: @*/
445: PetscErrorCode PetscSFCreateFromLayouts(PetscLayout rmap, PetscLayout lmap, PetscSF* sf)
446: {
448: PetscInt i,nroots,nleaves = 0;
449: PetscInt rN, lst, len;
450: PetscMPIInt owner = -1;
451: PetscSFNode *remote;
452: MPI_Comm rcomm = rmap->comm;
453: MPI_Comm lcomm = lmap->comm;
454: PetscMPIInt flg;
458: if (!rmap->setupcalled) SETERRQ(rcomm,PETSC_ERR_ARG_WRONGSTATE,"Root layout not setup");
459: if (!lmap->setupcalled) SETERRQ(lcomm,PETSC_ERR_ARG_WRONGSTATE,"Leaf layout not setup");
460: MPI_Comm_compare(rcomm,lcomm,&flg);
461: if (flg != MPI_CONGRUENT && flg != MPI_IDENT) SETERRQ(rcomm,PETSC_ERR_SUP,"cannot map two layouts with non-matching communicators");
462: PetscSFCreate(rcomm,sf);
463: PetscLayoutGetLocalSize(rmap,&nroots);
464: PetscLayoutGetSize(rmap,&rN);
465: PetscLayoutGetRange(lmap,&lst,&len);
466: PetscMalloc1(len-lst,&remote);
467: for (i = lst; i < len && i < rN; i++) {
468: if (owner < -1 || i >= rmap->range[owner+1]) {
469: PetscLayoutFindOwner(rmap,i,&owner);
470: }
471: remote[nleaves].rank = owner;
472: remote[nleaves].index = i - rmap->range[owner];
473: nleaves++;
474: }
475: PetscSFSetGraph(*sf,nroots,nleaves,NULL,PETSC_OWN_POINTER,remote,PETSC_COPY_VALUES);
476: PetscFree(remote);
477: return(0);
478: }
480: /* TODO: handle nooffprocentries like MatZeroRowsMapLocal_Private, since this code is the same */
481: PetscErrorCode PetscLayoutMapLocal(PetscLayout map,PetscInt N,const PetscInt idxs[], PetscInt *on,PetscInt **oidxs,PetscInt **ogidxs)
482: {
483: PetscInt *owners = map->range;
484: PetscInt n = map->n;
485: PetscSF sf;
486: PetscInt *lidxs,*work = NULL;
487: PetscSFNode *ridxs;
488: PetscMPIInt rank, p = 0;
489: PetscInt r, len = 0;
493: if (on) *on = 0; /* squelch -Wmaybe-uninitialized */
494: /* Create SF where leaves are input idxs and roots are owned idxs */
495: MPI_Comm_rank(map->comm,&rank);
496: PetscMalloc1(n,&lidxs);
497: for (r = 0; r < n; ++r) lidxs[r] = -1;
498: PetscMalloc1(N,&ridxs);
499: for (r = 0; r < N; ++r) {
500: const PetscInt idx = idxs[r];
501: if (idx < 0 || map->N <= idx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Index %D out of range [0,%D)",idx,map->N);
502: if (idx < owners[p] || owners[p+1] <= idx) { /* short-circuit the search if the last p owns this idx too */
503: PetscLayoutFindOwner(map,idx,&p);
504: }
505: ridxs[r].rank = p;
506: ridxs[r].index = idxs[r] - owners[p];
507: }
508: PetscSFCreate(map->comm,&sf);
509: PetscSFSetGraph(sf,n,N,NULL,PETSC_OWN_POINTER,ridxs,PETSC_OWN_POINTER);
510: PetscSFReduceBegin(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);
511: PetscSFReduceEnd(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);
512: if (ogidxs) { /* communicate global idxs */
513: PetscInt cum = 0,start,*work2;
515: PetscMalloc1(n,&work);
516: PetscCalloc1(N,&work2);
517: for (r = 0; r < N; ++r) if (idxs[r] >=0) cum++;
518: MPI_Scan(&cum,&start,1,MPIU_INT,MPI_SUM,map->comm);
519: start -= cum;
520: cum = 0;
521: for (r = 0; r < N; ++r) if (idxs[r] >=0) work2[r] = start+cum++;
522: PetscSFReduceBegin(sf,MPIU_INT,work2,work,MPI_REPLACE);
523: PetscSFReduceEnd(sf,MPIU_INT,work2,work,MPI_REPLACE);
524: PetscFree(work2);
525: }
526: PetscSFDestroy(&sf);
527: /* Compress and put in indices */
528: for (r = 0; r < n; ++r)
529: if (lidxs[r] >= 0) {
530: if (work) work[len] = work[r];
531: lidxs[len++] = r;
532: }
533: if (on) *on = len;
534: if (oidxs) *oidxs = lidxs;
535: if (ogidxs) *ogidxs = work;
536: return(0);
537: }
539: /*@
540: PetscSFCreateByMatchingIndices - Create SF by matching root and leaf indices
542: Collective
544: Input Parameters:
545: + layout - PetscLayout defining the global index space and the rank that brokers each index
546: . numRootIndices - size of rootIndices
547: . rootIndices - PetscInt array of global indices of which this process requests ownership
548: . rootLocalIndices - root local index permutation (NULL if no permutation)
549: . rootLocalOffset - offset to be added to root local indices
550: . numLeafIndices - size of leafIndices
551: . leafIndices - PetscInt array of global indices with which this process requires data associated
552: . leafLocalIndices - leaf local index permutation (NULL if no permutation)
553: - leafLocalOffset - offset to be added to leaf local indices
555: Output Parameters:
556: + sfA - star forest representing the communication pattern from the layout space to the leaf space (NULL if not needed)
557: - sf - star forest representing the communication pattern from the root space to the leaf space
559: Example 1:
560: $
561: $ rank : 0 1 2
562: $ rootIndices : [1 0 2] [3] [3]
563: $ rootLocalOffset : 100 200 300
564: $ layout : [0 1] [2] [3]
565: $ leafIndices : [0] [2] [0 3]
566: $ leafLocalOffset : 400 500 600
567: $
568: would build the following SF:
569: $
570: $ [0] 400 <- (0,101)
571: $ [1] 500 <- (0,102)
572: $ [2] 600 <- (0,101)
573: $ [2] 601 <- (2,300)
574: $
575: Example 2:
576: $
577: $ rank : 0 1 2
578: $ rootIndices : [1 0 2] [3] [3]
579: $ rootLocalOffset : 100 200 300
580: $ layout : [0 1] [2] [3]
581: $ leafIndices : rootIndices rootIndices rootIndices
582: $ leafLocalOffset : rootLocalOffset rootLocalOffset rootLocalOffset
583: $
584: would build the following SF:
585: $
586: $ [1] 200 <- (2,300)
587: $
588: Example 3:
589: $
590: $ No process requests ownership of global index 1, but no process needs it.
591: $
592: $ rank : 0 1 2
593: $ numRootIndices : 2 1 1
594: $ rootIndices : [0 2] [3] [3]
595: $ rootLocalOffset : 100 200 300
596: $ layout : [0 1] [2] [3]
597: $ numLeafIndices : 1 1 2
598: $ leafIndices : [0] [2] [0 3]
599: $ leafLocalOffset : 400 500 600
600: $
601: would build the following SF:
602: $
603: $ [0] 400 <- (0,100)
604: $ [1] 500 <- (0,101)
605: $ [2] 600 <- (0,100)
606: $ [2] 601 <- (2,300)
607: $
609: Notes:
610: The layout parameter represents any partitioning of [0, N), where N is the total number of global indices, and its
611: local size can be set to PETSC_DECIDE.
612: If a global index x lies in the partition owned by process i, each process whose rootIndices contains x requests
613: ownership of x and sends its own rank and the local index of x to process i.
614: If multiple processes request ownership of x, the one with the highest rank is to own x.
615: Process i then broadcasts the ownership information, so that each process whose leafIndices contains x knows the
616: ownership information of x.
617: The output sf is constructed by associating each leaf point to a root point in this way.
619: Suppose there is point data ordered according to the global indices and partitioned according to the given layout.
620: The optional output PetscSF object sfA can be used to push such data to leaf points.
622: All indices in rootIndices and leafIndices must lie in the layout range. The union (over all processes) of rootIndices
623: must cover that of leafIndices, but need not cover the entire layout.
625: If (leafIndices, leafLocalIndices, leafLocalOffset) == (rootIndices, rootLocalIndices, rootLocalOffset), the output
626: star forest is almost identity, so will only include non-trivial part of the map.
628: Developer Notes:
629: Current approach of a process of the highest rank gaining the ownership may cause load imbalance; consider using
630: hash(rank, root_local_index) as the bid for the ownership determination.
632: Level: advanced
634: .seealso: PetscSFCreate()
635: @*/
636: 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)
637: {
638: MPI_Comm comm = layout->comm;
639: PetscMPIInt size, rank;
640: PetscSF sf1;
641: PetscSFNode *owners, *buffer, *iremote;
642: PetscInt *ilocal, nleaves, N, n, i;
643: #if defined(PETSC_USE_DEBUG)
644: PetscInt N1;
645: #endif
646: PetscBool flag;
647: PetscErrorCode ierr;
656: if (numRootIndices < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "numRootIndices (%D) must be non-negative", numRootIndices);
657: if (numLeafIndices < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "numLeafIndices (%D) must be non-negative", numLeafIndices);
658: MPI_Comm_size(comm, &size);
659: MPI_Comm_rank(comm, &rank);
660: PetscLayoutSetUp(layout);
661: PetscLayoutGetSize(layout, &N);
662: PetscLayoutGetLocalSize(layout, &n);
663: flag = (PetscBool)(leafIndices == rootIndices);
664: MPIU_Allreduce(MPI_IN_PLACE, &flag, 1, MPIU_BOOL, MPI_LAND, comm);
665: if (flag && numLeafIndices != numRootIndices) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "leafIndices == rootIndices, but numLeafIndices (%D) != numRootIndices(%D)", numLeafIndices, numRootIndices);
666: #if defined(PETSC_USE_DEBUG)
667: N1 = PETSC_MIN_INT;
668: for (i = 0; i < numRootIndices; i++) if (rootIndices[i] > N1) N1 = rootIndices[i];
669: MPI_Allreduce(MPI_IN_PLACE, &N1, 1, MPIU_INT, MPI_MAX, comm);
670: if (N1 >= N) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Max. root index (%D) out of layout range [0,%D)", N1, N);
671: if (!flag) {
672: N1 = PETSC_MIN_INT;
673: for (i = 0; i < numLeafIndices; i++) if (leafIndices[i] > N1) N1 = leafIndices[i];
674: MPI_Allreduce(MPI_IN_PLACE, &N1, 1, MPIU_INT, MPI_MAX, comm);
675: if (N1 >= N) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Max. leaf index (%D) out of layout range [0,%D)", N1, N);
676: }
677: #endif
678: /* Reduce: owners -> buffer */
679: PetscMalloc1(n, &buffer);
680: PetscSFCreate(comm, &sf1);
681: PetscSFSetFromOptions(sf1);
682: PetscSFSetGraphLayout(sf1, layout, numRootIndices, NULL, PETSC_OWN_POINTER, rootIndices);
683: PetscMalloc1(numRootIndices, &owners);
684: for (i = 0; i < numRootIndices; ++i) {
685: owners[i].rank = rank;
686: owners[i].index = rootLocalOffset + (rootLocalIndices ? rootLocalIndices[i] : i);
687: }
688: for (i = 0; i < n; ++i) {
689: buffer[i].index = -1;
690: buffer[i].rank = -1;
691: }
692: PetscSFReduceBegin(sf1, MPIU_2INT, owners, buffer, MPI_MAXLOC);
693: PetscSFReduceEnd(sf1, MPIU_2INT, owners, buffer, MPI_MAXLOC);
694: /* Bcast: buffer -> owners */
695: if (!flag) {
696: /* leafIndices is different from rootIndices */
697: PetscFree(owners);
698: PetscSFSetGraphLayout(sf1, layout, numLeafIndices, NULL, PETSC_OWN_POINTER, leafIndices);
699: PetscMalloc1(numLeafIndices, &owners);
700: }
701: PetscSFBcastBegin(sf1, MPIU_2INT, buffer, owners, MPI_REPLACE);
702: PetscSFBcastEnd(sf1, MPIU_2INT, buffer, owners, MPI_REPLACE);
703: for (i = 0; i < numLeafIndices; ++i) if (owners[i].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Global point %D was unclaimed", leafIndices[i]);
704: PetscFree(buffer);
705: if (sfA) {*sfA = sf1;}
706: else {PetscSFDestroy(&sf1);}
707: /* Create sf */
708: if (flag && rootLocalIndices == leafLocalIndices && leafLocalOffset == rootLocalOffset) {
709: /* leaf space == root space */
710: for (i = 0, nleaves = 0; i < numLeafIndices; ++i) if (owners[i].rank != rank) ++nleaves;
711: PetscMalloc1(nleaves, &ilocal);
712: PetscMalloc1(nleaves, &iremote);
713: for (i = 0, nleaves = 0; i < numLeafIndices; ++i) {
714: if (owners[i].rank != rank) {
715: ilocal[nleaves] = leafLocalOffset + i;
716: iremote[nleaves].rank = owners[i].rank;
717: iremote[nleaves].index = owners[i].index;
718: ++nleaves;
719: }
720: }
721: PetscFree(owners);
722: } else {
723: nleaves = numLeafIndices;
724: PetscMalloc1(nleaves, &ilocal);
725: for (i = 0; i < nleaves; ++i) {ilocal[i] = leafLocalOffset + (leafLocalIndices ? leafLocalIndices[i] : i);}
726: iremote = owners;
727: }
728: PetscSFCreate(comm, sf);
729: PetscSFSetFromOptions(*sf);
730: PetscSFSetGraph(*sf, rootLocalOffset + numRootIndices, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
731: return(0);
732: }