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: }