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