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