Actual source code: sf.c

petsc-3.9.0 2018-04-07
Report Typos and Errors
  1:  #include <petsc/private/sfimpl.h>
  2:  #include <petscctable.h>

  4: #if defined(PETSC_USE_DEBUG)
  5: #  define PetscSFCheckGraphSet(sf,arg) do {                          \
  6:     if (PetscUnlikely(!(sf)->graphset))                              \
  7:       SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call PetscSFSetGraph() on argument %D \"%s\" before %s()",(arg),#sf,PETSC_FUNCTION_NAME); \
  8:   } while (0)
  9: #else
 10: #  define PetscSFCheckGraphSet(sf,arg) do {} while (0)
 11: #endif

 13: const char *const PetscSFDuplicateOptions[] = {"CONFONLY","RANKS","GRAPH","PetscSFDuplicateOption","PETSCSF_DUPLICATE_",0};

 15: /*@
 16:    PetscSFCreate - create a star forest communication context

 18:    Not Collective

 20:    Input Arguments:
 21: .  comm - communicator on which the star forest will operate

 23:    Output Arguments:
 24: .  sf - new star forest context

 26:    Level: intermediate

 28: .seealso: PetscSFSetGraph(), PetscSFDestroy()
 29: @*/
 30: PetscErrorCode PetscSFCreate(MPI_Comm comm,PetscSF *sf)
 31: {
 33:   PetscSF        b;

 37:   PetscSFInitializePackage();

 39:   PetscHeaderCreate(b,PETSCSF_CLASSID,"PetscSF","Star Forest","PetscSF",comm,PetscSFDestroy,PetscSFView);

 41:   b->nroots    = -1;
 42:   b->nleaves   = -1;
 43:   b->minleaf   = PETSC_MAX_INT;
 44:   b->maxleaf   = PETSC_MIN_INT;
 45:   b->nranks    = -1;
 46:   b->rankorder = PETSC_TRUE;
 47:   b->ingroup   = MPI_GROUP_NULL;
 48:   b->outgroup  = MPI_GROUP_NULL;
 49:   b->graphset  = PETSC_FALSE;

 51:   *sf = b;
 52:   return(0);
 53: }

 55: /*@
 56:    PetscSFReset - Reset a star forest so that different sizes or neighbors can be used

 58:    Collective

 60:    Input Arguments:
 61: .  sf - star forest

 63:    Level: advanced

 65: .seealso: PetscSFCreate(), PetscSFSetGraph(), PetscSFDestroy()
 66: @*/
 67: PetscErrorCode PetscSFReset(PetscSF sf)
 68: {

 73:   if (sf->ops->Reset) {(*sf->ops->Reset)(sf);}
 74:   sf->nroots   = -1;
 75:   sf->nleaves  = -1;
 76:   sf->minleaf  = PETSC_MAX_INT;
 77:   sf->maxleaf  = PETSC_MIN_INT;
 78:   sf->mine     = NULL;
 79:   sf->remote   = NULL;
 80:   sf->graphset = PETSC_FALSE;
 81:   PetscFree(sf->mine_alloc);
 82:   PetscFree(sf->remote_alloc);
 83:   sf->nranks = -1;
 84:   PetscFree4(sf->ranks,sf->roffset,sf->rmine,sf->rremote);
 85:   sf->degreeknown = PETSC_FALSE;
 86:   PetscFree(sf->degree);
 87:   if (sf->ingroup  != MPI_GROUP_NULL) {MPI_Group_free(&sf->ingroup);}
 88:   if (sf->outgroup != MPI_GROUP_NULL) {MPI_Group_free(&sf->outgroup);}
 89:   PetscSFDestroy(&sf->multi);
 90:   sf->setupcalled = PETSC_FALSE;
 91:   return(0);
 92: }

 94: /*@C
 95:    PetscSFSetType - Set the PetscSF communication implementation

 97:    Collective on PetscSF

 99:    Input Parameters:
100: +  sf - the PetscSF context
101: -  type - a known method

103:    Options Database Key:
104: .  -sf_type <type> - Sets the method; use -help for a list
105:    of available methods (for instance, window, pt2pt, neighbor)

107:    Notes:
108:    See "include/petscsf.h" for available methods (for instance)
109: +    PETSCSFWINDOW - MPI-2/3 one-sided
110: -    PETSCSFBASIC - basic implementation using MPI-1 two-sided

112:   Level: intermediate

114: .keywords: PetscSF, set, type

116: .seealso: PetscSFType, PetscSFCreate()
117: @*/
118: PetscErrorCode PetscSFSetType(PetscSF sf,PetscSFType type)
119: {
120:   PetscErrorCode ierr,(*r)(PetscSF);
121:   PetscBool      match;


127:   PetscObjectTypeCompare((PetscObject)sf,type,&match);
128:   if (match) return(0);

130:   PetscFunctionListFind(PetscSFList,type,&r);
131:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested PetscSF type %s",type);
132:   /* Destroy the previous PetscSF implementation context */
133:   if (sf->ops->Destroy) {(*(sf)->ops->Destroy)(sf);}
134:   PetscMemzero(sf->ops,sizeof(*sf->ops));
135:   PetscObjectChangeTypeName((PetscObject)sf,type);
136:   (*r)(sf);
137:   return(0);
138: }

140: /*@C
141:   PetscSFGetType - Get the PetscSF communication implementation

143:   Not Collective

145:   Input Parameter:
146: . sf  - the PetscSF context

148:   Output Parameter:
149: . type - the PetscSF type name

151:   Level: intermediate

153: .keywords: PetscSF, get, type
154: .seealso: PetscSFSetType(), PetscSFCreate()
155: @*/
156: PetscErrorCode PetscSFGetType(PetscSF sf, PetscSFType *type)
157: {
161:   *type = ((PetscObject)sf)->type_name;
162:   return(0);
163: }

165: /*@
166:    PetscSFDestroy - destroy star forest

168:    Collective

170:    Input Arguments:
171: .  sf - address of star forest

173:    Level: intermediate

175: .seealso: PetscSFCreate(), PetscSFReset()
176: @*/
177: PetscErrorCode PetscSFDestroy(PetscSF *sf)
178: {

182:   if (!*sf) return(0);
184:   if (--((PetscObject)(*sf))->refct > 0) {*sf = NULL; return(0);}
185:   PetscSFReset(*sf);
186:   if ((*sf)->ops->Destroy) {(*(*sf)->ops->Destroy)(*sf);}
187:   PetscHeaderDestroy(sf);
188:   return(0);
189: }

191: /*@
192:    PetscSFSetUp - set up communication structures

194:    Collective

196:    Input Arguments:
197: .  sf - star forest communication object

199:    Level: beginner

201: .seealso: PetscSFSetFromOptions(), PetscSFSetType()
202: @*/
203: PetscErrorCode PetscSFSetUp(PetscSF sf)
204: {

209:   PetscSFCheckGraphSet(sf,1);
210:   if (sf->setupcalled) return(0);
211:   if (!((PetscObject)sf)->type_name) {PetscSFSetType(sf,PETSCSFBASIC);}
212:   PetscLogEventBegin(PETSCSF_SetUp,sf,0,0,0);
213:   if (sf->ops->SetUp) {(*sf->ops->SetUp)(sf);}
214:   PetscLogEventEnd(PETSCSF_SetUp,sf,0,0,0);
215:   sf->setupcalled = PETSC_TRUE;
216:   return(0);
217: }

219: /*@
220:    PetscSFSetFromOptions - set PetscSF options using the options database

222:    Logically Collective

224:    Input Arguments:
225: .  sf - star forest

227:    Options Database Keys:
228: +  -sf_type - implementation type, see PetscSFSetType()
229: -  -sf_rank_order - sort composite points for gathers and scatters in rank order, gathers are non-deterministic otherwise

231:    Level: intermediate

233: .keywords: KSP, set, from, options, database

235: .seealso: PetscSFWindowSetSyncType()
236: @*/
237: PetscErrorCode PetscSFSetFromOptions(PetscSF sf)
238: {
239:   PetscSFType    deft;
240:   char           type[256];
242:   PetscBool      flg;

246:   PetscObjectOptionsBegin((PetscObject)sf);
247:   deft = ((PetscObject)sf)->type_name ? ((PetscObject)sf)->type_name : PETSCSFBASIC;
248:   PetscOptionsFList("-sf_type","PetscSF implementation type","PetscSFSetType",PetscSFList,deft,type,sizeof(type),&flg);
249:   PetscSFSetType(sf,flg ? type : deft);
250:   PetscOptionsBool("-sf_rank_order","sort composite points for gathers and scatters in rank order, gathers are non-deterministic otherwise","PetscSFSetRankOrder",sf->rankorder,&sf->rankorder,NULL);
251:   if (sf->ops->SetFromOptions) {(*sf->ops->SetFromOptions)(PetscOptionsObject,sf);}
252:   PetscOptionsEnd();
253:   return(0);
254: }

256: /*@
257:    PetscSFSetRankOrder - sort multi-points for gathers and scatters by rank order

259:    Logically Collective

261:    Input Arguments:
262: +  sf - star forest
263: -  flg - PETSC_TRUE to sort, PETSC_FALSE to skip sorting (lower setup cost, but non-deterministic)

265:    Level: advanced

267: .seealso: PetscSFGatherBegin(), PetscSFScatterBegin()
268: @*/
269: PetscErrorCode PetscSFSetRankOrder(PetscSF sf,PetscBool flg)
270: {

275:   if (sf->multi) SETERRQ(PetscObjectComm((PetscObject)sf),PETSC_ERR_ARG_WRONGSTATE,"Rank ordering must be set before first call to PetscSFGatherBegin() or PetscSFScatterBegin()");
276:   sf->rankorder = flg;
277:   return(0);
278: }

280: /*@
281:    PetscSFSetGraph - Set a parallel star forest

283:    Collective

285:    Input Arguments:
286: +  sf - star forest
287: .  nroots - number of root vertices on the current process (these are possible targets for other process to attach leaves)
288: .  nleaves - number of leaf vertices on the current process, each of these references a root on any process
289: .  ilocal - locations of leaves in leafdata buffers, pass NULL for contiguous storage
290: .  localmode - copy mode for ilocal
291: .  iremote - remote locations of root vertices for each leaf on the current process
292: -  remotemode - copy mode for iremote

294:    Level: intermediate

296:    Notes: In Fortran you must use PETSC_COPY_VALUES for localmode and remotemode

298:    Developers Note: Local indices which are the identity permutation in the range [0,nleaves) are discarded as they
299:    encode contiguous storage. In such case, if localmode is PETSC_OWN_POINTER, the memory is deallocated as it is not
300:    needed

302: .seealso: PetscSFCreate(), PetscSFView(), PetscSFGetGraph()
303: @*/
304: PetscErrorCode PetscSFSetGraph(PetscSF sf,PetscInt nroots,PetscInt nleaves,const PetscInt *ilocal,PetscCopyMode localmode,const PetscSFNode *iremote,PetscCopyMode remotemode)
305: {

312:   if (nroots  < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nroots %D, cannot be negative",nroots);
313:   if (nleaves < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nleaves %D, cannot be negative",nleaves);

315:   PetscSFReset(sf);
316:   PetscLogEventBegin(PETSCSF_SetGraph,sf,0,0,0);

318:   sf->nroots  = nroots;
319:   sf->nleaves = nleaves;

321:   if (nleaves && ilocal) {
322:     PetscInt i;
323:     PetscInt minleaf = PETSC_MAX_INT;
324:     PetscInt maxleaf = PETSC_MIN_INT;
325:     int      contiguous = 1;
326:     for (i=0; i<nleaves; i++) {
327:       minleaf = PetscMin(minleaf,ilocal[i]);
328:       maxleaf = PetscMax(maxleaf,ilocal[i]);
329:       contiguous &= (ilocal[i] == i);
330:     }
331:     sf->minleaf = minleaf;
332:     sf->maxleaf = maxleaf;
333:     if (contiguous) {
334:       if (localmode == PETSC_OWN_POINTER) {
335:         PetscFree(ilocal);
336:       }
337:       ilocal = NULL;
338:     }
339:   } else {
340:     sf->minleaf = 0;
341:     sf->maxleaf = nleaves - 1;
342:   }

344:   if (ilocal) {
345:     switch (localmode) {
346:     case PETSC_COPY_VALUES:
347:       PetscMalloc1(nleaves,&sf->mine_alloc);
348:       PetscMemcpy(sf->mine_alloc,ilocal,nleaves*sizeof(*ilocal));
349:       sf->mine = sf->mine_alloc;
350:       break;
351:     case PETSC_OWN_POINTER:
352:       sf->mine_alloc = (PetscInt*)ilocal;
353:       sf->mine       = sf->mine_alloc;
354:       break;
355:     case PETSC_USE_POINTER:
356:       sf->mine_alloc = NULL;
357:       sf->mine       = (PetscInt*)ilocal;
358:       break;
359:     default: SETERRQ(PetscObjectComm((PetscObject)sf),PETSC_ERR_ARG_OUTOFRANGE,"Unknown localmode");
360:     }
361:   }

363:   switch (remotemode) {
364:   case PETSC_COPY_VALUES:
365:     PetscMalloc1(nleaves,&sf->remote_alloc);
366:     PetscMemcpy(sf->remote_alloc,iremote,nleaves*sizeof(*iremote));
367:     sf->remote = sf->remote_alloc;
368:     break;
369:   case PETSC_OWN_POINTER:
370:     sf->remote_alloc = (PetscSFNode*)iremote;
371:     sf->remote       = sf->remote_alloc;
372:     break;
373:   case PETSC_USE_POINTER:
374:     sf->remote_alloc = NULL;
375:     sf->remote       = (PetscSFNode*)iremote;
376:     break;
377:   default: SETERRQ(PetscObjectComm((PetscObject)sf),PETSC_ERR_ARG_OUTOFRANGE,"Unknown remotemode");
378:   }

380:   PetscLogEventEnd(PETSCSF_SetGraph,sf,0,0,0);
381:   sf->graphset = PETSC_TRUE;
382:   return(0);
383: }

385: /*@
386:    PetscSFCreateInverseSF - given a PetscSF in which all vertices have degree 1, creates the inverse map

388:    Collective

390:    Input Arguments:
391: .  sf - star forest to invert

393:    Output Arguments:
394: .  isf - inverse of sf

396:    Level: advanced

398:    Notes:
399:    All roots must have degree 1.

401:    The local space may be a permutation, but cannot be sparse.

403: .seealso: PetscSFSetGraph()
404: @*/
405: PetscErrorCode PetscSFCreateInverseSF(PetscSF sf,PetscSF *isf)
406: {
408:   PetscMPIInt    rank;
409:   PetscInt       i,nroots,nleaves,maxlocal,count,*newilocal;
410:   const PetscInt *ilocal;
411:   PetscSFNode    *roots,*leaves;

415:   PetscSFCheckGraphSet(sf,1);

418:   PetscSFGetGraph(sf,&nroots,&nleaves,&ilocal,NULL);
419:   maxlocal = sf->maxleaf+1; /* TODO: We should use PetscSFGetLeafRange() */

421:   MPI_Comm_rank(PetscObjectComm((PetscObject)sf),&rank);
422:   PetscMalloc2(nroots,&roots,maxlocal,&leaves);
423:   for (i=0; i<maxlocal; i++) {
424:     leaves[i].rank  = rank;
425:     leaves[i].index = i;
426:   }
427:   for (i=0; i <nroots; i++) {
428:     roots[i].rank  = -1;
429:     roots[i].index = -1;
430:   }
431:   PetscSFReduceBegin(sf,MPIU_2INT,leaves,roots,MPIU_REPLACE);
432:   PetscSFReduceEnd(sf,MPIU_2INT,leaves,roots,MPIU_REPLACE);

434:   /* Check whether our leaves are sparse */
435:   for (i=0,count=0; i<nroots; i++) if (roots[i].rank >= 0) count++;
436:   if (count == nroots) newilocal = NULL;
437:   else {                        /* Index for sparse leaves and compact "roots" array (which is to become our leaves). */
438:     PetscMalloc1(count,&newilocal);
439:     for (i=0,count=0; i<nroots; i++) {
440:       if (roots[i].rank >= 0) {
441:         newilocal[count]   = i;
442:         roots[count].rank  = roots[i].rank;
443:         roots[count].index = roots[i].index;
444:         count++;
445:       }
446:     }
447:   }

449:   PetscSFDuplicate(sf,PETSCSF_DUPLICATE_CONFONLY,isf);
450:   PetscSFSetGraph(*isf,maxlocal,count,newilocal,PETSC_OWN_POINTER,roots,PETSC_COPY_VALUES);
451:   PetscFree2(roots,leaves);
452:   return(0);
453: }

455: /*@
456:    PetscSFDuplicate - duplicate a PetscSF, optionally preserving rank connectivity and graph

458:    Collective

460:    Input Arguments:
461: +  sf - communication object to duplicate
462: -  opt - PETSCSF_DUPLICATE_CONFONLY, PETSCSF_DUPLICATE_RANKS, or PETSCSF_DUPLICATE_GRAPH (see PetscSFDuplicateOption)

464:    Output Arguments:
465: .  newsf - new communication object

467:    Level: beginner

469: .seealso: PetscSFCreate(), PetscSFSetType(), PetscSFSetGraph()
470: @*/
471: PetscErrorCode PetscSFDuplicate(PetscSF sf,PetscSFDuplicateOption opt,PetscSF *newsf)
472: {
473:   PetscSFType    type;

480:   PetscSFCreate(PetscObjectComm((PetscObject)sf),newsf);
481:   PetscSFGetType(sf,&type);
482:   if (type) {PetscSFSetType(*newsf,type);}
483:   if (opt == PETSCSF_DUPLICATE_GRAPH) {
484:     PetscInt          nroots,nleaves;
485:     const PetscInt    *ilocal;
486:     const PetscSFNode *iremote;
487:     PetscSFCheckGraphSet(sf,1);
488:     PetscSFGetGraph(sf,&nroots,&nleaves,&ilocal,&iremote);
489:     PetscSFSetGraph(*newsf,nroots,nleaves,ilocal,PETSC_COPY_VALUES,iremote,PETSC_COPY_VALUES);
490:   }
491:   if (sf->ops->Duplicate) {(*sf->ops->Duplicate)(sf,opt,*newsf);}
492:   return(0);
493: }

495: /*@C
496:    PetscSFGetGraph - Get the graph specifying a parallel star forest

498:    Not Collective

500:    Input Arguments:
501: .  sf - star forest

503:    Output Arguments:
504: +  nroots - number of root vertices on the current process (these are possible targets for other process to attach leaves)
505: .  nleaves - number of leaf vertices on the current process, each of these references a root on any process
506: .  ilocal - locations of leaves in leafdata buffers
507: -  iremote - remote locations of root vertices for each leaf on the current process

509:    Notes:
510:    We are not currently requiring that the graph is set, thus returning nroots=-1 if it has not been set yet

512:    Level: intermediate

514: .seealso: PetscSFCreate(), PetscSFView(), PetscSFSetGraph()
515: @*/
516: PetscErrorCode PetscSFGetGraph(PetscSF sf,PetscInt *nroots,PetscInt *nleaves,const PetscInt **ilocal,const PetscSFNode **iremote)
517: {

521:   if (nroots) *nroots = sf->nroots;
522:   if (nleaves) *nleaves = sf->nleaves;
523:   if (ilocal) *ilocal = sf->mine;
524:   if (iremote) *iremote = sf->remote;
525:   return(0);
526: }

528: /*@
529:    PetscSFGetLeafRange - Get the active leaf ranges

531:    Not Collective

533:    Input Arguments:
534: .  sf - star forest

536:    Output Arguments:
537: +  minleaf - minimum active leaf on this process
538: -  maxleaf - maximum active leaf on this process

540:    Level: developer

542: .seealso: PetscSFCreate(), PetscSFView(), PetscSFSetGraph(), PetscSFGetGraph()
543: @*/
544: PetscErrorCode PetscSFGetLeafRange(PetscSF sf,PetscInt *minleaf,PetscInt *maxleaf)
545: {

549:   PetscSFCheckGraphSet(sf,1);
550:   if (minleaf) *minleaf = sf->minleaf;
551:   if (maxleaf) *maxleaf = sf->maxleaf;
552:   return(0);
553: }

555: /*@C
556:    PetscSFView - view a star forest

558:    Collective

560:    Input Arguments:
561: +  sf - star forest
562: -  viewer - viewer to display graph, for example PETSC_VIEWER_STDOUT_WORLD

564:    Level: beginner

566: .seealso: PetscSFCreate(), PetscSFSetGraph()
567: @*/
568: PetscErrorCode PetscSFView(PetscSF sf,PetscViewer viewer)
569: {
570:   PetscErrorCode    ierr;
571:   PetscBool         iascii;
572:   PetscViewerFormat format;

576:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sf),&viewer);}
579:   if (sf->graphset) {PetscSFSetUp(sf);}
580:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
581:   if (iascii) {
582:     PetscMPIInt rank;
583:     PetscInt    ii,i,j;

585:     PetscObjectPrintClassNamePrefixType((PetscObject)sf,viewer);
586:     PetscViewerASCIIPushTab(viewer);
587:     if (sf->ops->View) {(*sf->ops->View)(sf,viewer);}
588:     if (!sf->graphset) {
589:       PetscViewerASCIIPrintf(viewer,"PetscSFSetGraph() has not been called yet\n");
590:       PetscViewerASCIIPopTab(viewer);
591:       return(0);
592:     }
593:     MPI_Comm_rank(PetscObjectComm((PetscObject)sf),&rank);
594:     PetscViewerASCIIPushSynchronized(viewer);
595:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of roots=%D, leaves=%D, remote ranks=%D\n",rank,sf->nroots,sf->nleaves,sf->nranks);
596:     for (i=0; i<sf->nleaves; i++) {
597:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D <- (%D,%D)\n",rank,sf->mine ? sf->mine[i] : i,sf->remote[i].rank,sf->remote[i].index);
598:     }
599:     PetscViewerFlush(viewer);
600:     PetscViewerGetFormat(viewer,&format);
601:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
602:       PetscMPIInt *tmpranks,*perm;
603:       PetscMalloc2(sf->nranks,&tmpranks,sf->nranks,&perm);
604:       PetscMemcpy(tmpranks,sf->ranks,sf->nranks*sizeof(tmpranks[0]));
605:       for (i=0; i<sf->nranks; i++) perm[i] = i;
606:       PetscSortMPIIntWithArray(sf->nranks,tmpranks,perm);
607:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Roots referenced by my leaves, by rank\n",rank);
608:       for (ii=0; ii<sf->nranks; ii++) {
609:         i = perm[ii];
610:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %d: %D edges\n",rank,sf->ranks[i],sf->roffset[i+1]-sf->roffset[i]);
611:         for (j=sf->roffset[i]; j<sf->roffset[i+1]; j++) {
612:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d]    %D <- %D\n",rank,sf->rmine[j],sf->rremote[j]);
613:         }
614:       }
615:       PetscFree2(tmpranks,perm);
616:     }
617:     PetscViewerFlush(viewer);
618:     PetscViewerASCIIPopSynchronized(viewer);
619:     PetscViewerASCIIPopTab(viewer);
620:   }
621:   return(0);
622: }

624: /*@C
625:    PetscSFGetRanks - Get ranks and number of vertices referenced by leaves on this process

627:    Not Collective

629:    Input Arguments:
630: .  sf - star forest

632:    Output Arguments:
633: +  nranks - number of ranks referenced by local part
634: .  ranks - array of ranks
635: .  roffset - offset in rmine/rremote for each rank (length nranks+1)
636: .  rmine - concatenated array holding local indices referencing each remote rank
637: -  rremote - concatenated array holding remote indices referenced for each remote rank

639:    Level: developer

641: .seealso: PetscSFSetGraph()
642: @*/
643: PetscErrorCode PetscSFGetRanks(PetscSF sf,PetscInt *nranks,const PetscMPIInt **ranks,const PetscInt **roffset,const PetscInt **rmine,const PetscInt **rremote)
644: {

648:   if (!sf->setupcalled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call PetscSFSetUp() before obtaining ranks");
649:   if (nranks)  *nranks  = sf->nranks;
650:   if (ranks)   *ranks   = sf->ranks;
651:   if (roffset) *roffset = sf->roffset;
652:   if (rmine)   *rmine   = sf->rmine;
653:   if (rremote) *rremote = sf->rremote;
654:   return(0);
655: }

657: static PetscBool InList(PetscMPIInt needle,PetscMPIInt n,const PetscMPIInt *list) {
658:   PetscInt i;
659:   for (i=0; i<n; i++) {
660:     if (needle == list[i]) return PETSC_TRUE;
661:   }
662:   return PETSC_FALSE;
663: }

665: /*@C
666:    PetscSFSetUpRanks - Set up data structures associated with ranks; this is for internal use by PetscSF implementations.

668:    Collective

670:    Input Arguments:
671: +  sf - PetscSF to set up; PetscSFSetGraph() must have been called
672: -  dgroup - MPI_Group of ranks to be distinguished (e.g., for self or shared memory exchange)

674:    Level: developer

676: .seealso: PetscSFGetRanks()
677: @*/
678: PetscErrorCode PetscSFSetUpRanks(PetscSF sf,MPI_Group dgroup)
679: {
680:   PetscErrorCode     ierr;
681:   PetscTable         table;
682:   PetscTablePosition pos;
683:   PetscMPIInt        size,groupsize,*groupranks;
684:   PetscInt           i,*rcount,*ranks;

688:   PetscSFCheckGraphSet(sf,1);
689:   MPI_Comm_size(PetscObjectComm((PetscObject)sf),&size);
690:   PetscTableCreate(10,size,&table);
691:   for (i=0; i<sf->nleaves; i++) {
692:     /* Log 1-based rank */
693:     PetscTableAdd(table,sf->remote[i].rank+1,1,ADD_VALUES);
694:   }
695:   PetscTableGetCount(table,&sf->nranks);
696:   PetscMalloc4(sf->nranks,&sf->ranks,sf->nranks+1,&sf->roffset,sf->nleaves,&sf->rmine,sf->nleaves,&sf->rremote);
697:   PetscMalloc2(sf->nranks,&rcount,sf->nranks,&ranks);
698:   PetscTableGetHeadPosition(table,&pos);
699:   for (i=0; i<sf->nranks; i++) {
700:     PetscTableGetNext(table,&pos,&ranks[i],&rcount[i]);
701:     ranks[i]--;             /* Convert back to 0-based */
702:   }
703:   PetscTableDestroy(&table);

705:   /* We expect that dgroup is reliably "small" while nranks could be large */
706:   {
707:     MPI_Group group = MPI_GROUP_NULL;
708:     PetscMPIInt *dgroupranks;
709:     MPI_Comm_group(PetscObjectComm((PetscObject)sf),&group);
710:     MPI_Group_size(dgroup,&groupsize);
711:     PetscMalloc1(groupsize,&dgroupranks);
712:     PetscMalloc1(groupsize,&groupranks);
713:     for (i=0; i<groupsize; i++) dgroupranks[i] = i;
714:     if (groupsize) {MPI_Group_translate_ranks(dgroup,groupsize,dgroupranks,group,groupranks);}
715:     MPI_Group_free(&group);
716:     PetscFree(dgroupranks);
717:   }

719:   /* Partition ranks[] into distinguished (first sf->ndranks) followed by non-distinguished */
720:   for (sf->ndranks=0,i=sf->nranks; sf->ndranks<i; ) {
721:     for (i--; sf->ndranks<i; i--) { /* Scan i backward looking for distinguished rank */
722:       if (InList(ranks[i],groupsize,groupranks)) break;
723:     }
724:     for ( ; sf->ndranks<=i; sf->ndranks++) { /* Scan sf->ndranks forward looking for non-distinguished rank */
725:       if (!InList(ranks[sf->ndranks],groupsize,groupranks)) break;
726:     }
727:     if (sf->ndranks < i) {                         /* Swap ranks[sf->ndranks] with ranks[i] */
728:       PetscInt    tmprank,tmpcount;
729:       tmprank = ranks[i];
730:       tmpcount = rcount[i];
731:       ranks[i] = ranks[sf->ndranks];
732:       rcount[i] = rcount[sf->ndranks];
733:       ranks[sf->ndranks] = tmprank;
734:       rcount[sf->ndranks] = tmpcount;
735:       sf->ndranks++;
736:     }
737:   }
738:   PetscFree(groupranks);
739:   PetscSortIntWithArray(sf->ndranks,ranks,rcount);
740:   PetscSortIntWithArray(sf->nranks-sf->ndranks,ranks+sf->ndranks,rcount+sf->ndranks);
741:   sf->roffset[0] = 0;
742:   for (i=0; i<sf->nranks; i++) {
743:     PetscMPIIntCast(ranks[i],sf->ranks+i);
744:     sf->roffset[i+1] = sf->roffset[i] + rcount[i];
745:     rcount[i]        = 0;
746:   }
747:   for (i=0; i<sf->nleaves; i++) {
748:     PetscInt irank;
749:     /* Search for index of iremote[i].rank in sf->ranks */
750:     PetscFindMPIInt(sf->remote[i].rank,sf->ndranks,sf->ranks,&irank);
751:     if (irank < 0) {
752:       PetscFindMPIInt(sf->remote[i].rank,sf->nranks-sf->ndranks,sf->ranks+sf->ndranks,&irank);
753:       if (irank >= 0) irank += sf->ndranks;
754:     }
755:     if (irank < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Could not find rank %D in array",sf->remote[i].rank);
756:     sf->rmine[sf->roffset[irank] + rcount[irank]]   = sf->mine ? sf->mine[i] : i;
757:     sf->rremote[sf->roffset[irank] + rcount[irank]] = sf->remote[i].index;
758:     rcount[irank]++;
759:   }
760:   PetscFree2(rcount,ranks);
761:   return(0);
762: }

764: /*@C
765:    PetscSFGetGroups - gets incoming and outgoing process groups

767:    Collective

769:    Input Argument:
770: .  sf - star forest

772:    Output Arguments:
773: +  incoming - group of origin processes for incoming edges (leaves that reference my roots)
774: -  outgoing - group of destination processes for outgoing edges (roots that I reference)

776:    Level: developer

778: .seealso: PetscSFGetWindow(), PetscSFRestoreWindow()
779: @*/
780: PetscErrorCode PetscSFGetGroups(PetscSF sf,MPI_Group *incoming,MPI_Group *outgoing)
781: {
783:   MPI_Group      group = MPI_GROUP_NULL;

786:   if (!sf->setupcalled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call PetscSFSetUp() before obtaining groups");
787:   if (sf->ingroup == MPI_GROUP_NULL) {
788:     PetscInt       i;
789:     const PetscInt *indegree;
790:     PetscMPIInt    rank,*outranks,*inranks;
791:     PetscSFNode    *remote;
792:     PetscSF        bgcount;

794:     /* Compute the number of incoming ranks */
795:     PetscMalloc1(sf->nranks,&remote);
796:     for (i=0; i<sf->nranks; i++) {
797:       remote[i].rank  = sf->ranks[i];
798:       remote[i].index = 0;
799:     }
800:     PetscSFDuplicate(sf,PETSCSF_DUPLICATE_CONFONLY,&bgcount);
801:     PetscSFSetGraph(bgcount,1,sf->nranks,NULL,PETSC_COPY_VALUES,remote,PETSC_OWN_POINTER);
802:     PetscSFComputeDegreeBegin(bgcount,&indegree);
803:     PetscSFComputeDegreeEnd(bgcount,&indegree);

805:     /* Enumerate the incoming ranks */
806:     PetscMalloc2(indegree[0],&inranks,sf->nranks,&outranks);
807:     MPI_Comm_rank(PetscObjectComm((PetscObject)sf),&rank);
808:     for (i=0; i<sf->nranks; i++) outranks[i] = rank;
809:     PetscSFGatherBegin(bgcount,MPI_INT,outranks,inranks);
810:     PetscSFGatherEnd(bgcount,MPI_INT,outranks,inranks);
811:     MPI_Comm_group(PetscObjectComm((PetscObject)sf),&group);
812:     MPI_Group_incl(group,indegree[0],inranks,&sf->ingroup);
813:     MPI_Group_free(&group);
814:     PetscFree2(inranks,outranks);
815:     PetscSFDestroy(&bgcount);
816:   }
817:   *incoming = sf->ingroup;

819:   if (sf->outgroup == MPI_GROUP_NULL) {
820:     MPI_Comm_group(PetscObjectComm((PetscObject)sf),&group);
821:     MPI_Group_incl(group,sf->nranks,sf->ranks,&sf->outgroup);
822:     MPI_Group_free(&group);
823:   }
824:   *outgoing = sf->outgroup;
825:   return(0);
826: }

828: /*@
829:    PetscSFGetMultiSF - gets the inner SF implemeting gathers and scatters

831:    Collective

833:    Input Argument:
834: .  sf - star forest that may contain roots with 0 or with more than 1 vertex

836:    Output Arguments:
837: .  multi - star forest with split roots, such that each root has degree exactly 1

839:    Level: developer

841:    Notes:

843:    In most cases, users should use PetscSFGatherBegin() and PetscSFScatterBegin() instead of manipulating multi
844:    directly. Since multi satisfies the stronger condition that each entry in the global space has exactly one incoming
845:    edge, it is a candidate for future optimization that might involve its removal.

847: .seealso: PetscSFSetGraph(), PetscSFGatherBegin(), PetscSFScatterBegin()
848: @*/
849: PetscErrorCode PetscSFGetMultiSF(PetscSF sf,PetscSF *multi)
850: {

856:   if (sf->nroots < 0) {         /* Graph has not been set yet; why do we need this? */
857:     PetscSFDuplicate(sf,PETSCSF_DUPLICATE_RANKS,&sf->multi);
858:     *multi = sf->multi;
859:     return(0);
860:   }
861:   if (!sf->multi) {
862:     const PetscInt *indegree;
863:     PetscInt       i,*inoffset,*outones,*outoffset,maxlocal;
864:     PetscSFNode    *remote;
865:     maxlocal = sf->maxleaf+1; /* TODO: We should use PetscSFGetLeafRange() */
866:     PetscSFComputeDegreeBegin(sf,&indegree);
867:     PetscSFComputeDegreeEnd(sf,&indegree);
868:     PetscMalloc3(sf->nroots+1,&inoffset,maxlocal,&outones,maxlocal,&outoffset);
869:     inoffset[0] = 0;
870:     for (i=0; i<sf->nroots; i++) inoffset[i+1] = inoffset[i] + indegree[i];
871:     for (i=0; i<maxlocal; i++) outones[i] = 1;
872:     PetscSFFetchAndOpBegin(sf,MPIU_INT,inoffset,outones,outoffset,MPI_SUM);
873:     PetscSFFetchAndOpEnd(sf,MPIU_INT,inoffset,outones,outoffset,MPI_SUM);
874:     for (i=0; i<sf->nroots; i++) inoffset[i] -= indegree[i]; /* Undo the increment */
875: #if 0
876: #if defined(PETSC_USE_DEBUG)                                 /* Check that the expected number of increments occurred */
877:     for (i=0; i<sf->nroots; i++) {
878:       if (inoffset[i] + indegree[i] != inoffset[i+1]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Incorrect result after PetscSFFetchAndOp");
879:     }
880: #endif
881: #endif
882:     PetscMalloc1(sf->nleaves,&remote);
883:     for (i=0; i<sf->nleaves; i++) {
884:       remote[i].rank  = sf->remote[i].rank;
885:       remote[i].index = outoffset[sf->mine ? sf->mine[i] : i];
886:     }
887:     PetscSFDuplicate(sf,PETSCSF_DUPLICATE_RANKS,&sf->multi);
888:     PetscSFSetGraph(sf->multi,inoffset[sf->nroots],sf->nleaves,sf->mine,PETSC_COPY_VALUES,remote,PETSC_OWN_POINTER);
889:     if (sf->rankorder) {        /* Sort the ranks */
890:       PetscMPIInt rank;
891:       PetscInt    *inranks,*newoffset,*outranks,*newoutoffset,*tmpoffset,maxdegree;
892:       PetscSFNode *newremote;
893:       MPI_Comm_rank(PetscObjectComm((PetscObject)sf),&rank);
894:       for (i=0,maxdegree=0; i<sf->nroots; i++) maxdegree = PetscMax(maxdegree,indegree[i]);
895:       PetscMalloc5(sf->multi->nroots,&inranks,sf->multi->nroots,&newoffset,maxlocal,&outranks,maxlocal,&newoutoffset,maxdegree,&tmpoffset);
896:       for (i=0; i<maxlocal; i++) outranks[i] = rank;
897:       PetscSFReduceBegin(sf->multi,MPIU_INT,outranks,inranks,MPIU_REPLACE);
898:       PetscSFReduceEnd(sf->multi,MPIU_INT,outranks,inranks,MPIU_REPLACE);
899:       /* Sort the incoming ranks at each vertex, build the inverse map */
900:       for (i=0; i<sf->nroots; i++) {
901:         PetscInt j;
902:         for (j=0; j<indegree[i]; j++) tmpoffset[j] = j;
903:         PetscSortIntWithArray(indegree[i],inranks+inoffset[i],tmpoffset);
904:         for (j=0; j<indegree[i]; j++) newoffset[inoffset[i] + tmpoffset[j]] = inoffset[i] + j;
905:       }
906:       PetscSFBcastBegin(sf->multi,MPIU_INT,newoffset,newoutoffset);
907:       PetscSFBcastEnd(sf->multi,MPIU_INT,newoffset,newoutoffset);
908:       PetscMalloc1(sf->nleaves,&newremote);
909:       for (i=0; i<sf->nleaves; i++) {
910:         newremote[i].rank  = sf->remote[i].rank;
911:         newremote[i].index = newoutoffset[sf->mine ? sf->mine[i] : i];
912:       }
913:       PetscSFSetGraph(sf->multi,inoffset[sf->nroots],sf->nleaves,sf->mine,PETSC_COPY_VALUES,newremote,PETSC_OWN_POINTER);
914:       PetscFree5(inranks,newoffset,outranks,newoutoffset,tmpoffset);
915:     }
916:     PetscFree3(inoffset,outones,outoffset);
917:   }
918:   *multi = sf->multi;
919:   return(0);
920: }

922: /*@C
923:    PetscSFCreateEmbeddedSF - removes edges from all but the selected roots, does not remap indices

925:    Collective

927:    Input Arguments:
928: +  sf - original star forest
929: .  nroots - number of roots to select on this process
930: -  selected - selected roots on this process

932:    Output Arguments:
933: .  newsf - new star forest

935:    Level: advanced

937:    Note:
938:    To use the new PetscSF, it may be necessary to know the indices of the leaves that are still participating. This can
939:    be done by calling PetscSFGetGraph().

941: .seealso: PetscSFSetGraph(), PetscSFGetGraph()
942: @*/
943: PetscErrorCode PetscSFCreateEmbeddedSF(PetscSF sf,PetscInt nroots,const PetscInt *selected,PetscSF *newsf)
944: {
945:   PetscInt      *rootdata, *leafdata, *ilocal;
946:   PetscSFNode   *iremote;
947:   PetscInt       leafsize = 0, nleaves = 0, n, i;

952:   PetscSFCheckGraphSet(sf,1);
955:   if (sf->mine) for (i = 0; i < sf->nleaves; ++i) {leafsize = PetscMax(leafsize, sf->mine[i]+1);}
956:   else leafsize = sf->nleaves;
957:   PetscCalloc2(sf->nroots,&rootdata,leafsize,&leafdata);
958:   for (i=0; i<nroots; ++i) rootdata[selected[i]] = 1;
959:   PetscSFBcastBegin(sf,MPIU_INT,rootdata,leafdata);
960:   PetscSFBcastEnd(sf,MPIU_INT,rootdata,leafdata);
961:   for (i = 0; i < leafsize; ++i) nleaves += leafdata[i];
962:   PetscMalloc1(nleaves,&ilocal);
963:   PetscMalloc1(nleaves,&iremote);
964:   for (i = 0, n = 0; i < sf->nleaves; ++i) {
965:     const PetscInt lidx = sf->mine ? sf->mine[i] : i;

967:     if (leafdata[lidx]) {
968:       ilocal[n]        = lidx;
969:       iremote[n].rank  = sf->remote[i].rank;
970:       iremote[n].index = sf->remote[i].index;
971:       ++n;
972:     }
973:   }
974:   if (n != nleaves) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "There is a size mismatch in the SF embedding, %d != %d", n, nleaves);
975:   PetscSFDuplicate(sf,PETSCSF_DUPLICATE_RANKS,newsf);
976:   PetscSFSetGraph(*newsf,sf->nroots,nleaves,ilocal,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
977:   PetscFree2(rootdata,leafdata);
978:   return(0);
979: }

981: /*@C
982:   PetscSFCreateEmbeddedLeafSF - removes edges from all but the selected leaves, does not remap indices

984:   Collective

986:   Input Arguments:
987: + sf - original star forest
988: . nleaves - number of leaves to select on this process
989: - selected - selected leaves on this process

991:   Output Arguments:
992: .  newsf - new star forest

994:   Level: advanced

996: .seealso: PetscSFCreateEmbeddedSF(), PetscSFSetGraph(), PetscSFGetGraph()
997: @*/
998: PetscErrorCode PetscSFCreateEmbeddedLeafSF(PetscSF sf, PetscInt nleaves, const PetscInt *selected, PetscSF *newsf)
999: {
1000:   PetscSFNode   *iremote;
1001:   PetscInt      *ilocal;
1002:   PetscInt       i;

1007:   PetscSFCheckGraphSet(sf, 1);
1010:   PetscMalloc1(nleaves, &ilocal);
1011:   PetscMalloc1(nleaves, &iremote);
1012:   for (i = 0; i < nleaves; ++i) {
1013:     const PetscInt l = selected[i];

1015:     ilocal[i]        = sf->mine ? sf->mine[l] : l;
1016:     iremote[i].rank  = sf->remote[l].rank;
1017:     iremote[i].index = sf->remote[l].index;
1018:   }
1019:   PetscSFDuplicate(sf, PETSCSF_DUPLICATE_RANKS, newsf);
1020:   PetscSFSetGraph(*newsf, sf->nroots, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
1021:   return(0);
1022: }

1024: /*@C
1025:    PetscSFBcastBegin - begin pointwise broadcast to be concluded with call to PetscSFBcastEnd()

1027:    Collective on PetscSF

1029:    Input Arguments:
1030: +  sf - star forest on which to communicate
1031: .  unit - data type associated with each node
1032: -  rootdata - buffer to broadcast

1034:    Output Arguments:
1035: .  leafdata - buffer to update with values from each leaf's respective root

1037:    Level: intermediate

1039: .seealso: PetscSFCreate(), PetscSFSetGraph(), PetscSFView(), PetscSFBcastEnd(), PetscSFReduceBegin()
1040: @*/
1041: PetscErrorCode PetscSFBcastBegin(PetscSF sf,MPI_Datatype unit,const void *rootdata,void *leafdata)
1042: {

1047:   PetscSFSetUp(sf);
1048:   PetscLogEventBegin(PETSCSF_BcastBegin,sf,0,0,0);
1049:   (*sf->ops->BcastBegin)(sf,unit,rootdata,leafdata);
1050:   PetscLogEventEnd(PETSCSF_BcastBegin,sf,0,0,0);
1051:   return(0);
1052: }

1054: /*@C
1055:    PetscSFBcastEnd - end a broadcast operation started with PetscSFBcastBegin()

1057:    Collective

1059:    Input Arguments:
1060: +  sf - star forest
1061: .  unit - data type
1062: -  rootdata - buffer to broadcast

1064:    Output Arguments:
1065: .  leafdata - buffer to update with values from each leaf's respective root

1067:    Level: intermediate

1069: .seealso: PetscSFSetGraph(), PetscSFReduceEnd()
1070: @*/
1071: PetscErrorCode PetscSFBcastEnd(PetscSF sf,MPI_Datatype unit,const void *rootdata,void *leafdata)
1072: {

1077:   PetscSFSetUp(sf);
1078:   PetscLogEventBegin(PETSCSF_BcastEnd,sf,0,0,0);
1079:   (*sf->ops->BcastEnd)(sf,unit,rootdata,leafdata);
1080:   PetscLogEventEnd(PETSCSF_BcastEnd,sf,0,0,0);
1081:   return(0);
1082: }

1084: /*@C
1085:    PetscSFReduceBegin - begin reduction of leafdata into rootdata, to be completed with call to PetscSFReduceEnd()

1087:    Collective

1089:    Input Arguments:
1090: +  sf - star forest
1091: .  unit - data type
1092: .  leafdata - values to reduce
1093: -  op - reduction operation

1095:    Output Arguments:
1096: .  rootdata - result of reduction of values from all leaves of each root

1098:    Level: intermediate

1100: .seealso: PetscSFBcastBegin()
1101: @*/
1102: PetscErrorCode PetscSFReduceBegin(PetscSF sf,MPI_Datatype unit,const void *leafdata,void *rootdata,MPI_Op op)
1103: {

1108:   PetscSFSetUp(sf);
1109:   PetscLogEventBegin(PETSCSF_ReduceBegin,sf,0,0,0);
1110:   (sf->ops->ReduceBegin)(sf,unit,leafdata,rootdata,op);
1111:   PetscLogEventEnd(PETSCSF_ReduceBegin,sf,0,0,0);
1112:   return(0);
1113: }

1115: /*@C
1116:    PetscSFReduceEnd - end a reduction operation started with PetscSFReduceBegin()

1118:    Collective

1120:    Input Arguments:
1121: +  sf - star forest
1122: .  unit - data type
1123: .  leafdata - values to reduce
1124: -  op - reduction operation

1126:    Output Arguments:
1127: .  rootdata - result of reduction of values from all leaves of each root

1129:    Level: intermediate

1131: .seealso: PetscSFSetGraph(), PetscSFBcastEnd()
1132: @*/
1133: PetscErrorCode PetscSFReduceEnd(PetscSF sf,MPI_Datatype unit,const void *leafdata,void *rootdata,MPI_Op op)
1134: {

1139:   PetscSFSetUp(sf);
1140:   PetscLogEventBegin(PETSCSF_ReduceEnd,sf,0,0,0);
1141:   (*sf->ops->ReduceEnd)(sf,unit,leafdata,rootdata,op);
1142:   PetscLogEventEnd(PETSCSF_ReduceEnd,sf,0,0,0);
1143:   return(0);
1144: }

1146: /*@C
1147:    PetscSFComputeDegreeBegin - begin computation of degree for each root vertex, to be completed with PetscSFComputeDegreeEnd()

1149:    Collective

1151:    Input Arguments:
1152: .  sf - star forest

1154:    Output Arguments:
1155: .  degree - degree of each root vertex

1157:    Level: advanced

1159: .seealso: PetscSFGatherBegin()
1160: @*/
1161: PetscErrorCode PetscSFComputeDegreeBegin(PetscSF sf,const PetscInt **degree)
1162: {

1167:   PetscSFCheckGraphSet(sf,1);
1169:   if (!sf->degreeknown) {
1170:     PetscInt i, nroots = sf->nroots, maxlocal = sf->maxleaf+1;  /* TODO: We should use PetscSFGetLeafRange() */
1171:     if (sf->degree) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Calls to PetscSFComputeDegreeBegin() cannot be nested.");
1172:     PetscMalloc1(nroots,&sf->degree);
1173:     PetscMalloc1(PetscMax(maxlocal,1),&sf->degreetmp); /* allocate at least one entry, see check in PetscSFComputeDegreeEnd() */
1174:     for (i=0; i<nroots; i++) sf->degree[i] = 0;
1175:     for (i=0; i<maxlocal; i++) sf->degreetmp[i] = 1;
1176:     PetscSFReduceBegin(sf,MPIU_INT,sf->degreetmp,sf->degree,MPI_SUM);
1177:   }
1178:   *degree = NULL;
1179:   return(0);
1180: }

1182: /*@C
1183:    PetscSFComputeDegreeEnd - complete computation of degree for each root vertex, started with PetscSFComputeDegreeBegin()

1185:    Collective

1187:    Input Arguments:
1188: .  sf - star forest

1190:    Output Arguments:
1191: .  degree - degree of each root vertex

1193:    Level: developer

1195: .seealso:
1196: @*/
1197: PetscErrorCode PetscSFComputeDegreeEnd(PetscSF sf,const PetscInt **degree)
1198: {

1203:   PetscSFCheckGraphSet(sf,1);
1205:   if (!sf->degreeknown) {
1206:     if (!sf->degreetmp) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PetscSFComputeDegreeBegin() before PetscSFComputeDegreeEnd()");
1207:     PetscSFReduceEnd(sf,MPIU_INT,sf->degreetmp,sf->degree,MPI_SUM);
1208:     PetscFree(sf->degreetmp);
1209:     sf->degreeknown = PETSC_TRUE;
1210:   }
1211:   *degree = sf->degree;
1212:   return(0);
1213: }

1215: /*@C
1216:    PetscSFFetchAndOpBegin - begin operation that fetches values from root and updates atomically by applying operation using my leaf value, to be completed with PetscSFFetchAndOpEnd()

1218:    Collective

1220:    Input Arguments:
1221: +  sf - star forest
1222: .  unit - data type
1223: .  leafdata - leaf values to use in reduction
1224: -  op - operation to use for reduction

1226:    Output Arguments:
1227: +  rootdata - root values to be updated, input state is seen by first process to perform an update
1228: -  leafupdate - state at each leaf's respective root immediately prior to my atomic update

1230:    Level: advanced

1232:    Note:
1233:    The update is only atomic at the granularity provided by the hardware. Different roots referenced by the same process
1234:    might be updated in a different order. Furthermore, if a composite type is used for the unit datatype, atomicity is
1235:    not guaranteed across the whole vertex. Therefore, this function is mostly only used with primitive types such as
1236:    integers.

1238: .seealso: PetscSFComputeDegreeBegin(), PetscSFReduceBegin(), PetscSFSetGraph()
1239: @*/
1240: PetscErrorCode PetscSFFetchAndOpBegin(PetscSF sf,MPI_Datatype unit,void *rootdata,const void *leafdata,void *leafupdate,MPI_Op op)
1241: {

1246:   PetscSFSetUp(sf);
1247:   PetscLogEventBegin(PETSCSF_FetchAndOpBegin,sf,0,0,0);
1248:   (*sf->ops->FetchAndOpBegin)(sf,unit,rootdata,leafdata,leafupdate,op);
1249:   PetscLogEventEnd(PETSCSF_FetchAndOpBegin,sf,0,0,0);
1250:   return(0);
1251: }

1253: /*@C
1254:    PetscSFFetchAndOpEnd - end operation started in matching call to PetscSFFetchAndOpBegin() to fetch values from roots and update atomically by applying operation using my leaf value

1256:    Collective

1258:    Input Arguments:
1259: +  sf - star forest
1260: .  unit - data type
1261: .  leafdata - leaf values to use in reduction
1262: -  op - operation to use for reduction

1264:    Output Arguments:
1265: +  rootdata - root values to be updated, input state is seen by first process to perform an update
1266: -  leafupdate - state at each leaf's respective root immediately prior to my atomic update

1268:    Level: advanced

1270: .seealso: PetscSFComputeDegreeEnd(), PetscSFReduceEnd(), PetscSFSetGraph()
1271: @*/
1272: PetscErrorCode PetscSFFetchAndOpEnd(PetscSF sf,MPI_Datatype unit,void *rootdata,const void *leafdata,void *leafupdate,MPI_Op op)
1273: {

1278:   PetscSFSetUp(sf);
1279:   PetscLogEventBegin(PETSCSF_FetchAndOpEnd,sf,0,0,0);
1280:   (*sf->ops->FetchAndOpEnd)(sf,unit,rootdata,leafdata,leafupdate,op);
1281:   PetscLogEventEnd(PETSCSF_FetchAndOpEnd,sf,0,0,0);
1282:   return(0);
1283: }

1285: /*@C
1286:    PetscSFGatherBegin - begin pointwise gather of all leaves into multi-roots, to be completed with PetscSFGatherEnd()

1288:    Collective

1290:    Input Arguments:
1291: +  sf - star forest
1292: .  unit - data type
1293: -  leafdata - leaf data to gather to roots

1295:    Output Argument:
1296: .  multirootdata - root buffer to gather into, amount of space per root is equal to its degree

1298:    Level: intermediate

1300: .seealso: PetscSFComputeDegreeBegin(), PetscSFScatterBegin()
1301: @*/
1302: PetscErrorCode PetscSFGatherBegin(PetscSF sf,MPI_Datatype unit,const void *leafdata,void *multirootdata)
1303: {
1305:   PetscSF        multi;

1309:   PetscSFSetUp(sf);
1310:   PetscSFGetMultiSF(sf,&multi);
1311:   PetscSFReduceBegin(multi,unit,leafdata,multirootdata,MPIU_REPLACE);
1312:   return(0);
1313: }

1315: /*@C
1316:    PetscSFGatherEnd - ends pointwise gather operation that was started with PetscSFGatherBegin()

1318:    Collective

1320:    Input Arguments:
1321: +  sf - star forest
1322: .  unit - data type
1323: -  leafdata - leaf data to gather to roots

1325:    Output Argument:
1326: .  multirootdata - root buffer to gather into, amount of space per root is equal to its degree

1328:    Level: intermediate

1330: .seealso: PetscSFComputeDegreeEnd(), PetscSFScatterEnd()
1331: @*/
1332: PetscErrorCode PetscSFGatherEnd(PetscSF sf,MPI_Datatype unit,const void *leafdata,void *multirootdata)
1333: {
1335:   PetscSF        multi;

1339:   PetscSFSetUp(sf);
1340:   PetscSFGetMultiSF(sf,&multi);
1341:   PetscSFReduceEnd(multi,unit,leafdata,multirootdata,MPIU_REPLACE);
1342:   return(0);
1343: }

1345: /*@C
1346:    PetscSFScatterBegin - begin pointwise scatter operation from multi-roots to leaves, to be completed with PetscSFScatterEnd()

1348:    Collective

1350:    Input Arguments:
1351: +  sf - star forest
1352: .  unit - data type
1353: -  multirootdata - root buffer to send to each leaf, one unit of data per leaf

1355:    Output Argument:
1356: .  leafdata - leaf data to be update with personal data from each respective root

1358:    Level: intermediate

1360: .seealso: PetscSFComputeDegreeBegin(), PetscSFScatterBegin()
1361: @*/
1362: PetscErrorCode PetscSFScatterBegin(PetscSF sf,MPI_Datatype unit,const void *multirootdata,void *leafdata)
1363: {
1365:   PetscSF        multi;

1369:   PetscSFSetUp(sf);
1370:   PetscSFGetMultiSF(sf,&multi);
1371:   PetscSFBcastBegin(multi,unit,multirootdata,leafdata);
1372:   return(0);
1373: }

1375: /*@C
1376:    PetscSFScatterEnd - ends pointwise scatter operation that was started with PetscSFScatterBegin()

1378:    Collective

1380:    Input Arguments:
1381: +  sf - star forest
1382: .  unit - data type
1383: -  multirootdata - root buffer to send to each leaf, one unit of data per leaf

1385:    Output Argument:
1386: .  leafdata - leaf data to be update with personal data from each respective root

1388:    Level: intermediate

1390: .seealso: PetscSFComputeDegreeEnd(), PetscSFScatterEnd()
1391: @*/
1392: PetscErrorCode PetscSFScatterEnd(PetscSF sf,MPI_Datatype unit,const void *multirootdata,void *leafdata)
1393: {
1395:   PetscSF        multi;

1399:   PetscSFSetUp(sf);
1400:   PetscSFGetMultiSF(sf,&multi);
1401:   PetscSFBcastEnd(multi,unit,multirootdata,leafdata);
1402:   return(0);
1403: }

1405: /*@
1406:   PetscSFCompose - Compose a new PetscSF equivalent to action to PetscSFs

1408:   Input Parameters:
1409: + sfA - The first PetscSF
1410: - sfB - The second PetscSF

1412:   Output Parameters:
1413: . sfBA - equvalent PetscSF for applying A then B

1415:   Level: developer

1417: .seealso: PetscSF, PetscSFGetGraph(), PetscSFSetGraph()
1418: @*/
1419: PetscErrorCode PetscSFCompose(PetscSF sfA, PetscSF sfB, PetscSF *sfBA)
1420: {
1421:   MPI_Comm           comm;
1422:   const PetscSFNode *remotePointsA, *remotePointsB;
1423:   PetscSFNode       *remotePointsBA;
1424:   const PetscInt    *localPointsA, *localPointsB;
1425:   PetscInt           numRootsA, numLeavesA, numRootsB, numLeavesB;
1426:   PetscErrorCode     ierr;

1430:   PetscSFCheckGraphSet(sfA, 1);
1432:   PetscSFCheckGraphSet(sfB, 2);
1434:   PetscObjectGetComm((PetscObject) sfA, &comm);
1435:   PetscSFGetGraph(sfA, &numRootsA, &numLeavesA, &localPointsA, &remotePointsA);
1436:   PetscSFGetGraph(sfB, &numRootsB, &numLeavesB, &localPointsB, &remotePointsB);
1437:   PetscMalloc1(numLeavesB, &remotePointsBA);
1438:   PetscSFBcastBegin(sfB, MPIU_2INT, remotePointsA, remotePointsBA);
1439:   PetscSFBcastEnd(sfB, MPIU_2INT, remotePointsA, remotePointsBA);
1440:   PetscSFCreate(comm, sfBA);
1441:   PetscSFSetGraph(*sfBA, numRootsA, numLeavesB, localPointsB, PETSC_COPY_VALUES, remotePointsBA, PETSC_OWN_POINTER);
1442:   return(0);
1443: }