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 Arguments:
10: + sf - star forest
11: . layout - PetscLayout defining the global space
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 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 (gdof < 0) {
125: for (d = 0; d < gsize; ++d, ++l) {
126: PetscInt offset = -(goff+1) + d, r;
128: PetscFindInt(offset,size+1,ranges,&r);
129: if (r < 0) r = -(r+2);
130: 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);
131: remote[l].rank = r;
132: remote[l].index = offset - ranges[r];
133: }
134: } else {
135: for (d = 0; d < gsize; ++d, ++l) {
136: remote[l].rank = rank;
137: remote[l].index = goff+d - ranges[rank];
138: }
139: }
140: }
141: if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %D != nleaves %D", l, nleaves);
142: PetscLayoutDestroy(&layout);
143: PetscSFSetGraph(sf, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
144: return(0);
145: }
147: static PetscErrorCode PetscSectionCheckConstraints_Static(PetscSection s)
148: {
152: if (!s->bc) {
153: PetscSectionCreate(PETSC_COMM_SELF, &s->bc);
154: PetscSectionSetChart(s->bc, s->pStart, s->pEnd);
155: }
156: return(0);
157: }
159: /*@C
160: PetscSFDistributeSection - Create a new PetscSection reorganized, moving from the root to the leaves of the SF
162: Collective on sf
164: Input Parameters:
165: + sf - The SF
166: - rootSection - Section defined on root space
168: Output Parameters:
169: + remoteOffsets - root offsets in leaf storage, or NULL
170: - leafSection - Section defined on the leaf space
172: Level: advanced
174: .seealso: PetscSFCreate()
175: @*/
176: PetscErrorCode PetscSFDistributeSection(PetscSF sf, PetscSection rootSection, PetscInt **remoteOffsets, PetscSection leafSection)
177: {
178: PetscSF embedSF;
179: const PetscInt *indices;
180: IS selected;
181: PetscInt numFields, nroots, rpStart, rpEnd, lpStart = PETSC_MAX_INT, lpEnd = -1, f, c;
182: PetscBool *sub, hasc;
186: PetscLogEventBegin(PETSCSF_DistSect,sf,0,0,0);
187: PetscSectionGetNumFields(rootSection, &numFields);
188: if (numFields) {PetscSectionSetNumFields(leafSection, numFields);}
189: PetscMalloc1(numFields+2, &sub);
190: sub[1] = rootSection->bc ? PETSC_TRUE : PETSC_FALSE;
191: for (f = 0; f < numFields; ++f) {
192: PetscSectionSym sym;
193: const char *name = NULL;
194: PetscInt numComp = 0;
196: sub[2 + f] = rootSection->field[f]->bc ? PETSC_TRUE : PETSC_FALSE;
197: PetscSectionGetFieldComponents(rootSection, f, &numComp);
198: PetscSectionGetFieldName(rootSection, f, &name);
199: PetscSectionGetFieldSym(rootSection, f, &sym);
200: PetscSectionSetFieldComponents(leafSection, f, numComp);
201: PetscSectionSetFieldName(leafSection, f, name);
202: PetscSectionSetFieldSym(leafSection, f, sym);
203: for (c = 0; c < rootSection->numFieldComponents[f]; ++c) {
204: PetscSectionGetComponentName(rootSection, f, c, &name);
205: PetscSectionSetComponentName(leafSection, f, c, name);
206: }
207: }
208: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
209: PetscSFGetGraph(sf,&nroots,NULL,NULL,NULL);
210: rpEnd = PetscMin(rpEnd,nroots);
211: rpEnd = PetscMax(rpStart,rpEnd);
212: /* see if we can avoid creating the embedded SF, since it can cost more than an allreduce */
213: sub[0] = (PetscBool)(nroots != rpEnd - rpStart);
214: MPIU_Allreduce(MPI_IN_PLACE, sub, 2+numFields, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)sf));
215: if (sub[0]) {
216: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
217: ISGetIndices(selected, &indices);
218: PetscSFCreateEmbeddedRootSF(sf, rpEnd - rpStart, indices, &embedSF);
219: ISRestoreIndices(selected, &indices);
220: ISDestroy(&selected);
221: } else {
222: PetscObjectReference((PetscObject)sf);
223: embedSF = sf;
224: }
225: PetscSFGetLeafRange(embedSF, &lpStart, &lpEnd);
226: lpEnd++;
228: PetscSectionSetChart(leafSection, lpStart, lpEnd);
230: /* Constrained dof section */
231: hasc = sub[1];
232: for (f = 0; f < numFields; ++f) hasc = (PetscBool)(hasc || sub[2+f]);
234: /* Could fuse these at the cost of copies and extra allocation */
235: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart],MPI_REPLACE);
236: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasDof[-rpStart], &leafSection->atlasDof[-lpStart],MPI_REPLACE);
237: if (sub[1]) {
238: PetscSectionCheckConstraints_Static(rootSection);
239: PetscSectionCheckConstraints_Static(leafSection);
240: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart],MPI_REPLACE);
241: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasDof[-rpStart], &leafSection->bc->atlasDof[-lpStart],MPI_REPLACE);
242: }
243: for (f = 0; f < numFields; ++f) {
244: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart],MPI_REPLACE);
245: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->atlasDof[-rpStart], &leafSection->field[f]->atlasDof[-lpStart],MPI_REPLACE);
246: if (sub[2+f]) {
247: PetscSectionCheckConstraints_Static(rootSection->field[f]);
248: PetscSectionCheckConstraints_Static(leafSection->field[f]);
249: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart],MPI_REPLACE);
250: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasDof[-rpStart], &leafSection->field[f]->bc->atlasDof[-lpStart],MPI_REPLACE);
251: }
252: }
253: if (remoteOffsets) {
254: PetscMalloc1(lpEnd - lpStart, remoteOffsets);
255: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart],MPI_REPLACE);
256: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart],MPI_REPLACE);
257: }
258: PetscSectionSetUp(leafSection);
259: if (hasc) { /* need to communicate bcIndices */
260: PetscSF bcSF;
261: PetscInt *rOffBc;
263: PetscMalloc1(lpEnd - lpStart, &rOffBc);
264: if (sub[1]) {
265: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart],MPI_REPLACE);
266: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->bc->atlasOff[-rpStart], &rOffBc[-lpStart],MPI_REPLACE);
267: PetscSFCreateSectionSF(embedSF, rootSection->bc, rOffBc, leafSection->bc, &bcSF);
268: PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices,MPI_REPLACE);
269: PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->bcIndices, leafSection->bcIndices,MPI_REPLACE);
270: PetscSFDestroy(&bcSF);
271: }
272: for (f = 0; f < numFields; ++f) {
273: if (sub[2+f]) {
274: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart],MPI_REPLACE);
275: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->field[f]->bc->atlasOff[-rpStart], &rOffBc[-lpStart],MPI_REPLACE);
276: PetscSFCreateSectionSF(embedSF, rootSection->field[f]->bc, rOffBc, leafSection->field[f]->bc, &bcSF);
277: PetscSFBcastBegin(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices,MPI_REPLACE);
278: PetscSFBcastEnd(bcSF, MPIU_INT, rootSection->field[f]->bcIndices, leafSection->field[f]->bcIndices,MPI_REPLACE);
279: PetscSFDestroy(&bcSF);
280: }
281: }
282: PetscFree(rOffBc);
283: }
284: PetscSFDestroy(&embedSF);
285: PetscFree(sub);
286: PetscLogEventEnd(PETSCSF_DistSect,sf,0,0,0);
287: return(0);
288: }
290: /*@C
291: PetscSFCreateRemoteOffsets - Create offsets for point data on remote processes
293: Collective on sf
295: Input Parameters:
296: + sf - The SF
297: . rootSection - Data layout of remote points for outgoing data (this is layout for SF roots)
298: - leafSection - Data layout of local points for incoming data (this is layout for SF leaves)
300: Output Parameter:
301: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
303: Level: developer
305: .seealso: PetscSFCreate()
306: @*/
307: PetscErrorCode PetscSFCreateRemoteOffsets(PetscSF sf, PetscSection rootSection, PetscSection leafSection, PetscInt **remoteOffsets)
308: {
309: PetscSF embedSF;
310: const PetscInt *indices;
311: IS selected;
312: PetscInt numRoots, rpStart = 0, rpEnd = 0, lpStart = 0, lpEnd = 0;
313: PetscErrorCode ierr;
316: *remoteOffsets = NULL;
317: PetscSFGetGraph(sf, &numRoots, NULL, NULL, NULL);
318: if (numRoots < 0) return(0);
319: PetscLogEventBegin(PETSCSF_RemoteOff,sf,0,0,0);
320: PetscSectionGetChart(rootSection, &rpStart, &rpEnd);
321: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
322: ISCreateStride(PETSC_COMM_SELF, rpEnd - rpStart, rpStart, 1, &selected);
323: ISGetIndices(selected, &indices);
324: PetscSFCreateEmbeddedRootSF(sf, rpEnd - rpStart, indices, &embedSF);
325: ISRestoreIndices(selected, &indices);
326: ISDestroy(&selected);
327: PetscCalloc1(lpEnd - lpStart, remoteOffsets);
328: PetscSFBcastBegin(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart],MPI_REPLACE);
329: PetscSFBcastEnd(embedSF, MPIU_INT, &rootSection->atlasOff[-rpStart], &(*remoteOffsets)[-lpStart],MPI_REPLACE);
330: PetscSFDestroy(&embedSF);
331: PetscLogEventEnd(PETSCSF_RemoteOff,sf,0,0,0);
332: return(0);
333: }
335: /*@C
336: PetscSFCreateSectionSF - Create an expanded SF of dofs, assuming the input SF relates points
338: Collective on sf
340: Input Parameters:
341: + sf - The SF
342: . rootSection - Data layout of remote points for outgoing data (this is usually the serial section)
343: . remoteOffsets - Offsets for point data on remote processes (these are offsets from the root section), or NULL
344: - leafSection - Data layout of local points for incoming data (this is the distributed section)
346: Output Parameters:
347: - sectionSF - The new SF
349: Note: Either rootSection or remoteOffsets can be specified
351: Level: advanced
353: .seealso: PetscSFCreate()
354: @*/
355: PetscErrorCode PetscSFCreateSectionSF(PetscSF sf, PetscSection rootSection, PetscInt remoteOffsets[], PetscSection leafSection, PetscSF *sectionSF)
356: {
357: MPI_Comm comm;
358: const PetscInt *localPoints;
359: const PetscSFNode *remotePoints;
360: PetscInt lpStart, lpEnd;
361: PetscInt numRoots, numSectionRoots, numPoints, numIndices = 0;
362: PetscInt *localIndices;
363: PetscSFNode *remoteIndices;
364: PetscInt i, ind;
365: PetscErrorCode ierr;
373: PetscObjectGetComm((PetscObject)sf,&comm);
374: PetscSFCreate(comm, sectionSF);
375: PetscSectionGetChart(leafSection, &lpStart, &lpEnd);
376: PetscSectionGetStorageSize(rootSection, &numSectionRoots);
377: PetscSFGetGraph(sf, &numRoots, &numPoints, &localPoints, &remotePoints);
378: if (numRoots < 0) return(0);
379: PetscLogEventBegin(PETSCSF_SectSF,sf,0,0,0);
380: for (i = 0; i < numPoints; ++i) {
381: PetscInt localPoint = localPoints ? localPoints[i] : i;
382: PetscInt dof;
384: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
385: PetscSectionGetDof(leafSection, localPoint, &dof);
386: numIndices += dof;
387: }
388: }
389: PetscMalloc1(numIndices, &localIndices);
390: PetscMalloc1(numIndices, &remoteIndices);
391: /* Create new index graph */
392: for (i = 0, ind = 0; i < numPoints; ++i) {
393: PetscInt localPoint = localPoints ? localPoints[i] : i;
394: PetscInt rank = remotePoints[i].rank;
396: if ((localPoint >= lpStart) && (localPoint < lpEnd)) {
397: PetscInt remoteOffset = remoteOffsets[localPoint-lpStart];
398: PetscInt loff, dof, d;
400: PetscSectionGetOffset(leafSection, localPoint, &loff);
401: PetscSectionGetDof(leafSection, localPoint, &dof);
402: for (d = 0; d < dof; ++d, ++ind) {
403: localIndices[ind] = loff+d;
404: remoteIndices[ind].rank = rank;
405: remoteIndices[ind].index = remoteOffset+d;
406: }
407: }
408: }
409: if (numIndices != ind) SETERRQ2(comm, PETSC_ERR_PLIB, "Inconsistency in indices, %D should be %D", ind, numIndices);
410: PetscSFSetGraph(*sectionSF, numSectionRoots, numIndices, localIndices, PETSC_OWN_POINTER, remoteIndices, PETSC_OWN_POINTER);
411: PetscSFSetUp(*sectionSF);
412: PetscLogEventEnd(PETSCSF_SectSF,sf,0,0,0);
413: return(0);
414: }
416: /*@C
417: PetscSFCreateFromLayouts - Creates a parallel star forest mapping two PetscLayout objects
419: Collective
421: Input Arguments:
422: + rmap - PetscLayout defining the global root space
423: - lmap - PetscLayout defining the global leaf space
425: Output Arguments:
426: . sf - The parallel star forest
428: Level: intermediate
430: .seealso: PetscSFCreate(), PetscLayoutCreate(), PetscSFSetGraphLayout()
431: @*/
432: PetscErrorCode PetscSFCreateFromLayouts(PetscLayout rmap, PetscLayout lmap, PetscSF* sf)
433: {
435: PetscInt i,nroots,nleaves = 0;
436: PetscInt rN, lst, len;
437: PetscMPIInt owner = -1;
438: PetscSFNode *remote;
439: MPI_Comm rcomm = rmap->comm;
440: MPI_Comm lcomm = lmap->comm;
441: PetscMPIInt flg;
445: if (!rmap->setupcalled) SETERRQ(rcomm,PETSC_ERR_ARG_WRONGSTATE,"Root layout not setup");
446: if (!lmap->setupcalled) SETERRQ(lcomm,PETSC_ERR_ARG_WRONGSTATE,"Leaf layout not setup");
447: MPI_Comm_compare(rcomm,lcomm,&flg);
448: if (flg != MPI_CONGRUENT && flg != MPI_IDENT) SETERRQ(rcomm,PETSC_ERR_SUP,"cannot map two layouts with non-matching communicators");
449: PetscSFCreate(rcomm,sf);
450: PetscLayoutGetLocalSize(rmap,&nroots);
451: PetscLayoutGetSize(rmap,&rN);
452: PetscLayoutGetRange(lmap,&lst,&len);
453: PetscMalloc1(len-lst,&remote);
454: for (i = lst; i < len && i < rN; i++) {
455: if (owner < -1 || i >= rmap->range[owner+1]) {
456: PetscLayoutFindOwner(rmap,i,&owner);
457: }
458: remote[nleaves].rank = owner;
459: remote[nleaves].index = i - rmap->range[owner];
460: nleaves++;
461: }
462: PetscSFSetGraph(*sf,nroots,nleaves,NULL,PETSC_OWN_POINTER,remote,PETSC_COPY_VALUES);
463: PetscFree(remote);
464: return(0);
465: }
467: /* TODO: handle nooffprocentries like MatZeroRowsMapLocal_Private, since this code is the same */
468: PetscErrorCode PetscLayoutMapLocal(PetscLayout map,PetscInt N,const PetscInt idxs[], PetscInt *on,PetscInt **oidxs,PetscInt **ogidxs)
469: {
470: PetscInt *owners = map->range;
471: PetscInt n = map->n;
472: PetscSF sf;
473: PetscInt *lidxs,*work = NULL;
474: PetscSFNode *ridxs;
475: PetscMPIInt rank, p = 0;
476: PetscInt r, len = 0;
480: if (on) *on = 0; /* squelch -Wmaybe-uninitialized */
481: /* Create SF where leaves are input idxs and roots are owned idxs */
482: MPI_Comm_rank(map->comm,&rank);
483: PetscMalloc1(n,&lidxs);
484: for (r = 0; r < n; ++r) lidxs[r] = -1;
485: PetscMalloc1(N,&ridxs);
486: for (r = 0; r < N; ++r) {
487: const PetscInt idx = idxs[r];
488: if (idx < 0 || map->N <= idx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Index %D out of range [0,%D)",idx,map->N);
489: if (idx < owners[p] || owners[p+1] <= idx) { /* short-circuit the search if the last p owns this idx too */
490: PetscLayoutFindOwner(map,idx,&p);
491: }
492: ridxs[r].rank = p;
493: ridxs[r].index = idxs[r] - owners[p];
494: }
495: PetscSFCreate(map->comm,&sf);
496: PetscSFSetGraph(sf,n,N,NULL,PETSC_OWN_POINTER,ridxs,PETSC_OWN_POINTER);
497: PetscSFReduceBegin(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);
498: PetscSFReduceEnd(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);
499: if (ogidxs) { /* communicate global idxs */
500: PetscInt cum = 0,start,*work2;
502: PetscMalloc1(n,&work);
503: PetscCalloc1(N,&work2);
504: for (r = 0; r < N; ++r) if (idxs[r] >=0) cum++;
505: MPI_Scan(&cum,&start,1,MPIU_INT,MPI_SUM,map->comm);
506: start -= cum;
507: cum = 0;
508: for (r = 0; r < N; ++r) if (idxs[r] >=0) work2[r] = start+cum++;
509: PetscSFReduceBegin(sf,MPIU_INT,work2,work,MPI_REPLACE);
510: PetscSFReduceEnd(sf,MPIU_INT,work2,work,MPI_REPLACE);
511: PetscFree(work2);
512: }
513: PetscSFDestroy(&sf);
514: /* Compress and put in indices */
515: for (r = 0; r < n; ++r)
516: if (lidxs[r] >= 0) {
517: if (work) work[len] = work[r];
518: lidxs[len++] = r;
519: }
520: if (on) *on = len;
521: if (oidxs) *oidxs = lidxs;
522: if (ogidxs) *ogidxs = work;
523: return(0);
524: }
526: /*@
527: PetscSFCreateByMatchingIndices - Create SF by matching root and leaf indices
529: Collective
531: Input Parameters:
532: + layout - PetscLayout defining the global index space and the rank that brokers each index
533: . numRootIndices - size of rootIndices
534: . rootIndices - PetscInt array of global indices of which this process requests ownership
535: . rootLocalIndices - root local index permutation (NULL if no permutation)
536: . rootLocalOffset - offset to be added to root local indices
537: . numLeafIndices - size of leafIndices
538: . leafIndices - PetscInt array of global indices with which this process requires data associated
539: . leafLocalIndices - leaf local index permutation (NULL if no permutation)
540: - leafLocalOffset - offset to be added to leaf local indices
542: Output Parameter:
543: + sfA - star forest representing the communication pattern from the layout space to the leaf space (NULL if not needed)
544: - sf - star forest representing the communication pattern from the root space to the leaf space
546: Example 1:
547: $
548: $ rank : 0 1 2
549: $ rootIndices : [1 0 2] [3] [3]
550: $ rootLocalOffset : 100 200 300
551: $ layout : [0 1] [2] [3]
552: $ leafIndices : [0] [2] [0 3]
553: $ leafLocalOffset : 400 500 600
554: $
555: would build the following SF:
556: $
557: $ [0] 400 <- (0,101)
558: $ [1] 500 <- (0,102)
559: $ [2] 600 <- (0,101)
560: $ [2] 601 <- (2,300)
561: $
562: Example 2:
563: $
564: $ rank : 0 1 2
565: $ rootIndices : [1 0 2] [3] [3]
566: $ rootLocalOffset : 100 200 300
567: $ layout : [0 1] [2] [3]
568: $ leafIndices : rootIndices rootIndices rootIndices
569: $ leafLocalOffset : rootLocalOffset rootLocalOffset rootLocalOffset
570: $
571: would build the following SF:
572: $
573: $ [1] 200 <- (2,300)
574: $
575: Example 3:
576: $
577: $ No process requests ownership of global index 1, but no process needs it.
578: $
579: $ rank : 0 1 2
580: $ numRootIndices : 2 1 1
581: $ rootIndices : [0 2] [3] [3]
582: $ rootLocalOffset : 100 200 300
583: $ layout : [0 1] [2] [3]
584: $ numLeafIndices : 1 1 2
585: $ leafIndices : [0] [2] [0 3]
586: $ leafLocalOffset : 400 500 600
587: $
588: would build the following SF:
589: $
590: $ [0] 400 <- (0,100)
591: $ [1] 500 <- (0,101)
592: $ [2] 600 <- (0,100)
593: $ [2] 601 <- (2,300)
594: $
596: Notes:
597: The layout parameter represents any partitioning of [0, N), where N is the total number of global indices, and its
598: local size can be set to PETSC_DECIDE.
599: If a global index x lies in the partition owned by process i, each process whose rootIndices contains x requests
600: ownership of x and sends its own rank and the local index of x to process i.
601: If multiple processes request ownership of x, the one with the highest rank is to own x.
602: Process i then broadcasts the ownership information, so that each process whose leafIndices contains x knows the
603: ownership information of x.
604: The output sf is constructed by associating each leaf point to a root point in this way.
606: Suppose there is point data ordered according to the global indices and partitioned according to the given layout.
607: The optional output PetscSF object sfA can be used to push such data to leaf points.
609: All indices in rootIndices and leafIndices must lie in the layout range. The union (over all processes) of rootIndices
610: must cover that of leafIndices, but need not cover the entire layout.
612: If (leafIndices, leafLocalIndices, leafLocalOffset) == (rootIndices, rootLocalIndices, rootLocalOffset), the output
613: star forest is almost identity, so will only include non-trivial part of the map.
615: Developer Notes:
616: Current approach of a process of the highest rank gaining the ownership may cause load imbalance; consider using
617: hash(rank, root_local_index) as the bid for the ownership determination.
619: Level: advanced
621: .seealso: PetscSFCreate()
622: @*/
623: 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)
624: {
625: MPI_Comm comm = layout->comm;
626: PetscMPIInt size, rank;
627: PetscSF sf1;
628: PetscSFNode *owners, *buffer, *iremote;
629: PetscInt *ilocal, nleaves, N, n, i;
630: #if defined(PETSC_USE_DEBUG)
631: PetscInt N1;
632: #endif
633: PetscBool flag;
634: PetscErrorCode ierr;
643: if (numRootIndices < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "numRootIndices (%D) must be non-negative", numRootIndices);
644: if (numLeafIndices < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "numLeafIndices (%D) must be non-negative", numLeafIndices);
645: MPI_Comm_size(comm, &size);
646: MPI_Comm_rank(comm, &rank);
647: PetscLayoutSetUp(layout);
648: PetscLayoutGetSize(layout, &N);
649: PetscLayoutGetLocalSize(layout, &n);
650: flag = (PetscBool)(leafIndices == rootIndices);
651: MPIU_Allreduce(MPI_IN_PLACE, &flag, 1, MPIU_BOOL, MPI_LAND, comm);
652: if (flag && numLeafIndices != numRootIndices) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "leafIndices == rootIndices, but numLeafIndices (%D) != numRootIndices(%D)", numLeafIndices, numRootIndices);
653: #if defined(PETSC_USE_DEBUG)
654: N1 = PETSC_MIN_INT;
655: for (i = 0; i < numRootIndices; i++) if (rootIndices[i] > N1) N1 = rootIndices[i];
656: MPI_Allreduce(MPI_IN_PLACE, &N1, 1, MPIU_INT, MPI_MAX, comm);
657: if (N1 >= N) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Max. root index (%D) out of layout range [0,%D)", N1, N);
658: if (!flag) {
659: N1 = PETSC_MIN_INT;
660: for (i = 0; i < numLeafIndices; i++) if (leafIndices[i] > N1) N1 = leafIndices[i];
661: MPI_Allreduce(MPI_IN_PLACE, &N1, 1, MPIU_INT, MPI_MAX, comm);
662: if (N1 >= N) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Max. leaf index (%D) out of layout range [0,%D)", N1, N);
663: }
664: #endif
665: /* Reduce: owners -> buffer */
666: PetscMalloc1(n, &buffer);
667: PetscSFCreate(comm, &sf1);
668: PetscSFSetFromOptions(sf1);
669: PetscSFSetGraphLayout(sf1, layout, numRootIndices, NULL, PETSC_OWN_POINTER, rootIndices);
670: PetscMalloc1(numRootIndices, &owners);
671: for (i = 0; i < numRootIndices; ++i) {
672: owners[i].rank = rank;
673: owners[i].index = rootLocalOffset + (rootLocalIndices ? rootLocalIndices[i] : i);
674: }
675: for (i = 0; i < n; ++i) {
676: buffer[i].index = -1;
677: buffer[i].rank = -1;
678: }
679: PetscSFReduceBegin(sf1, MPIU_2INT, owners, buffer, MPI_MAXLOC);
680: PetscSFReduceEnd(sf1, MPIU_2INT, owners, buffer, MPI_MAXLOC);
681: /* Bcast: buffer -> owners */
682: if (!flag) {
683: /* leafIndices is different from rootIndices */
684: PetscFree(owners);
685: PetscSFSetGraphLayout(sf1, layout, numLeafIndices, NULL, PETSC_OWN_POINTER, leafIndices);
686: PetscMalloc1(numLeafIndices, &owners);
687: }
688: PetscSFBcastBegin(sf1, MPIU_2INT, buffer, owners, MPI_REPLACE);
689: PetscSFBcastEnd(sf1, MPIU_2INT, buffer, owners, MPI_REPLACE);
690: 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]);
691: PetscFree(buffer);
692: if (sfA) {*sfA = sf1;}
693: else {PetscSFDestroy(&sf1);}
694: /* Create sf */
695: if (flag && rootLocalIndices == leafLocalIndices && leafLocalOffset == rootLocalOffset) {
696: /* leaf space == root space */
697: for (i = 0, nleaves = 0; i < numLeafIndices; ++i) if (owners[i].rank != rank) ++nleaves;
698: PetscMalloc1(nleaves, &ilocal);
699: PetscMalloc1(nleaves, &iremote);
700: for (i = 0, nleaves = 0; i < numLeafIndices; ++i) {
701: if (owners[i].rank != rank) {
702: ilocal[nleaves] = leafLocalOffset + i;
703: iremote[nleaves].rank = owners[i].rank;
704: iremote[nleaves].index = owners[i].index;
705: ++nleaves;
706: }
707: }
708: PetscFree(owners);
709: } else {
710: nleaves = numLeafIndices;
711: PetscMalloc1(nleaves, &ilocal);
712: for (i = 0; i < nleaves; ++i) {ilocal[i] = leafLocalOffset + (leafLocalIndices ? leafLocalIndices[i] : i);}
713: iremote = owners;
714: }
715: PetscSFCreate(comm, sf);
716: PetscSFSetFromOptions(*sf);
717: PetscSFSetGraph(*sf, rootLocalOffset + numRootIndices, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
718: return(0);
719: }