Actual source code: sf.c

petsc-master 2019-07-23
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:    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: .seealso: PetscSFType, PetscSFCreate()
115: @*/
116: PetscErrorCode PetscSFSetType(PetscSF sf,PetscSFType type)
117: {
118:   PetscErrorCode ierr,(*r)(PetscSF);
119:   PetscBool      match;


125:   PetscObjectTypeCompare((PetscObject)sf,type,&match);
126:   if (match) return(0);

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

138: /*@C
139:   PetscSFGetType - Get the PetscSF communication implementation

141:   Not Collective

143:   Input Parameter:
144: . sf  - the PetscSF context

146:   Output Parameter:
147: . type - the PetscSF type name

149:   Level: intermediate

151: .seealso: PetscSFSetType(), PetscSFCreate()
152: @*/
153: PetscErrorCode PetscSFGetType(PetscSF sf, PetscSFType *type)
154: {
158:   *type = ((PetscObject)sf)->type_name;
159:   return(0);
160: }

162: /*@
163:    PetscSFDestroy - destroy star forest

165:    Collective

167:    Input Arguments:
168: .  sf - address of star forest

170:    Level: intermediate

172: .seealso: PetscSFCreate(), PetscSFReset()
173: @*/
174: PetscErrorCode PetscSFDestroy(PetscSF *sf)
175: {

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

188: /*@
189:    PetscSFSetUp - set up communication structures

191:    Collective

193:    Input Arguments:
194: .  sf - star forest communication object

196:    Level: beginner

198: .seealso: PetscSFSetFromOptions(), PetscSFSetType()
199: @*/
200: PetscErrorCode PetscSFSetUp(PetscSF sf)
201: {

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

216: /*@
217:    PetscSFSetFromOptions - set PetscSF options using the options database

219:    Logically Collective

221:    Input Arguments:
222: .  sf - star forest

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

228:    Level: intermediate

230: .seealso: PetscSFWindowSetSyncType()
231: @*/
232: PetscErrorCode PetscSFSetFromOptions(PetscSF sf)
233: {
234:   PetscSFType    deft;
235:   char           type[256];
237:   PetscBool      flg;

241:   PetscObjectOptionsBegin((PetscObject)sf);
242:   deft = ((PetscObject)sf)->type_name ? ((PetscObject)sf)->type_name : PETSCSFBASIC;
243:   PetscOptionsFList("-sf_type","PetscSF implementation type","PetscSFSetType",PetscSFList,deft,type,sizeof(type),&flg);
244:   PetscSFSetType(sf,flg ? type : deft);
245:   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);
246:   if (sf->ops->SetFromOptions) {(*sf->ops->SetFromOptions)(PetscOptionsObject,sf);}
247:   PetscOptionsEnd();
248:   return(0);
249: }

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

254:    Logically Collective

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

260:    Level: advanced

262: .seealso: PetscSFGatherBegin(), PetscSFScatterBegin()
263: @*/
264: PetscErrorCode PetscSFSetRankOrder(PetscSF sf,PetscBool flg)
265: {

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

275: /*@
276:    PetscSFSetGraph - Set a parallel star forest

278:    Collective

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

289:    Level: intermediate

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

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

298: .seealso: PetscSFCreate(), PetscSFView(), PetscSFGetGraph()
299: @*/
300: PetscErrorCode PetscSFSetGraph(PetscSF sf,PetscInt nroots,PetscInt nleaves,const PetscInt *ilocal,PetscCopyMode localmode,const PetscSFNode *iremote,PetscCopyMode remotemode)
301: {

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

311:   PetscSFReset(sf);
312:   PetscLogEventBegin(PETSCSF_SetGraph,sf,0,0,0);

314:   sf->nroots  = nroots;
315:   sf->nleaves = nleaves;

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

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

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

376:   PetscLogEventEnd(PETSCSF_SetGraph,sf,0,0,0);
377:   sf->graphset = PETSC_TRUE;
378:   return(0);
379: }

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

384:    Collective

386:    Input Arguments:
387: .  sf - star forest to invert

389:    Output Arguments:
390: .  isf - inverse of sf

392:    Level: advanced

394:    Notes:
395:    All roots must have degree 1.

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

399: .seealso: PetscSFSetGraph()
400: @*/
401: PetscErrorCode PetscSFCreateInverseSF(PetscSF sf,PetscSF *isf)
402: {
404:   PetscMPIInt    rank;
405:   PetscInt       i,nroots,nleaves,maxlocal,count,*newilocal;
406:   const PetscInt *ilocal;
407:   PetscSFNode    *roots,*leaves;

411:   PetscSFCheckGraphSet(sf,1);

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

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

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

445:   PetscSFDuplicate(sf,PETSCSF_DUPLICATE_CONFONLY,isf);
446:   PetscSFSetGraph(*isf,maxlocal,count,newilocal,PETSC_OWN_POINTER,roots,PETSC_COPY_VALUES);
447:   PetscFree2(roots,leaves);
448:   return(0);
449: }

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

454:    Collective

456:    Input Arguments:
457: +  sf - communication object to duplicate
458: -  opt - PETSCSF_DUPLICATE_CONFONLY, PETSCSF_DUPLICATE_RANKS, or PETSCSF_DUPLICATE_GRAPH (see PetscSFDuplicateOption)

460:    Output Arguments:
461: .  newsf - new communication object

463:    Level: beginner

465: .seealso: PetscSFCreate(), PetscSFSetType(), PetscSFSetGraph()
466: @*/
467: PetscErrorCode PetscSFDuplicate(PetscSF sf,PetscSFDuplicateOption opt,PetscSF *newsf)
468: {
469:   PetscSFType    type;

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

491: /*@C
492:    PetscSFGetGraph - Get the graph specifying a parallel star forest

494:    Not Collective

496:    Input Arguments:
497: .  sf - star forest

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

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

508:    Level: intermediate

510: .seealso: PetscSFCreate(), PetscSFView(), PetscSFSetGraph()
511: @*/
512: PetscErrorCode PetscSFGetGraph(PetscSF sf,PetscInt *nroots,PetscInt *nleaves,const PetscInt **ilocal,const PetscSFNode **iremote)
513: {

517:   if (nroots) *nroots = sf->nroots;
518:   if (nleaves) *nleaves = sf->nleaves;
519:   if (ilocal) *ilocal = sf->mine;
520:   if (iremote) *iremote = sf->remote;
521:   return(0);
522: }

524: /*@
525:    PetscSFGetLeafRange - Get the active leaf ranges

527:    Not Collective

529:    Input Arguments:
530: .  sf - star forest

532:    Output Arguments:
533: +  minleaf - minimum active leaf on this process
534: -  maxleaf - maximum active leaf on this process

536:    Level: developer

538: .seealso: PetscSFCreate(), PetscSFView(), PetscSFSetGraph(), PetscSFGetGraph()
539: @*/
540: PetscErrorCode PetscSFGetLeafRange(PetscSF sf,PetscInt *minleaf,PetscInt *maxleaf)
541: {

545:   PetscSFCheckGraphSet(sf,1);
546:   if (minleaf) *minleaf = sf->minleaf;
547:   if (maxleaf) *maxleaf = sf->maxleaf;
548:   return(0);
549: }

551: /*@C
552:    PetscSFView - view a star forest

554:    Collective

556:    Input Arguments:
557: +  sf - star forest
558: -  viewer - viewer to display graph, for example PETSC_VIEWER_STDOUT_WORLD

560:    Level: beginner

562: .seealso: PetscSFCreate(), PetscSFSetGraph()
563: @*/
564: PetscErrorCode PetscSFView(PetscSF sf,PetscViewer viewer)
565: {
566:   PetscErrorCode    ierr;
567:   PetscBool         iascii;
568:   PetscViewerFormat format;

572:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sf),&viewer);}
575:   if (sf->graphset) {PetscSFSetUp(sf);}
576:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
577:   if (iascii) {
578:     PetscMPIInt rank;
579:     PetscInt    ii,i,j;

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

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

623:    Not Collective

625:    Input Arguments:
626: .  sf - star forest

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

635:    Level: developer

637: .seealso: PetscSFSetGraph()
638: @*/
639: PetscErrorCode PetscSFGetRanks(PetscSF sf,PetscInt *nranks,const PetscMPIInt **ranks,const PetscInt **roffset,const PetscInt **rmine,const PetscInt **rremote)
640: {

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

653: /*@C
654:    PetscSFGetLeafRanks - Get leaf ranks referencing roots on this process

656:    Not Collective

658:    Input Arguments:
659: .  sf - star forest

661:    Output Arguments:
662: +  niranks - number of leaf ranks referencing roots on this process
663: .  iranks - array of ranks
664: .  ioffset - offset in irootloc for each rank (length niranks+1)
665: -  irootloc - concatenated array holding local indices of roots referenced by each leaf rank

667:    Level: developer

669: .seealso: PetscSFGetRanks()
670: @*/
671: PetscErrorCode PetscSFGetLeafRanks(PetscSF sf,PetscInt *niranks,const PetscMPIInt **iranks,const PetscInt **ioffset,const PetscInt **irootloc)
672: {

677:   if (!sf->setupcalled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call PetscSFSetUp() before obtaining ranks");
678:   if (sf->ops->GetLeafRanks) {
679:     (sf->ops->GetLeafRanks)(sf,niranks,iranks,ioffset,irootloc);
680:   } else {
681:     PetscSFType type;
682:     PetscSFGetType(sf,&type);
683:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"PetscSFGetLeafRanks() is not supported on this StarForest type: %s", type);
684:   }
685:   return(0);
686: }

688: static PetscBool InList(PetscMPIInt needle,PetscMPIInt n,const PetscMPIInt *list) {
689:   PetscInt i;
690:   for (i=0; i<n; i++) {
691:     if (needle == list[i]) return PETSC_TRUE;
692:   }
693:   return PETSC_FALSE;
694: }

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

699:    Collective

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

705:    Level: developer

707: .seealso: PetscSFGetRanks()
708: @*/
709: PetscErrorCode PetscSFSetUpRanks(PetscSF sf,MPI_Group dgroup)
710: {
711:   PetscErrorCode     ierr;
712:   PetscTable         table;
713:   PetscTablePosition pos;
714:   PetscMPIInt        size,groupsize,*groupranks;
715:   PetscInt           *rcount,*ranks;
716:   PetscInt           i, irank = -1,orank = -1;

720:   PetscSFCheckGraphSet(sf,1);
721:   MPI_Comm_size(PetscObjectComm((PetscObject)sf),&size);
722:   PetscTableCreate(10,size,&table);
723:   for (i=0; i<sf->nleaves; i++) {
724:     /* Log 1-based rank */
725:     PetscTableAdd(table,sf->remote[i].rank+1,1,ADD_VALUES);
726:   }
727:   PetscTableGetCount(table,&sf->nranks);
728:   PetscMalloc4(sf->nranks,&sf->ranks,sf->nranks+1,&sf->roffset,sf->nleaves,&sf->rmine,sf->nleaves,&sf->rremote);
729:   PetscMalloc2(sf->nranks,&rcount,sf->nranks,&ranks);
730:   PetscTableGetHeadPosition(table,&pos);
731:   for (i=0; i<sf->nranks; i++) {
732:     PetscTableGetNext(table,&pos,&ranks[i],&rcount[i]);
733:     ranks[i]--;             /* Convert back to 0-based */
734:   }
735:   PetscTableDestroy(&table);

737:   /* We expect that dgroup is reliably "small" while nranks could be large */
738:   {
739:     MPI_Group group = MPI_GROUP_NULL;
740:     PetscMPIInt *dgroupranks;
741:     MPI_Comm_group(PetscObjectComm((PetscObject)sf),&group);
742:     MPI_Group_size(dgroup,&groupsize);
743:     PetscMalloc1(groupsize,&dgroupranks);
744:     PetscMalloc1(groupsize,&groupranks);
745:     for (i=0; i<groupsize; i++) dgroupranks[i] = i;
746:     if (groupsize) {MPI_Group_translate_ranks(dgroup,groupsize,dgroupranks,group,groupranks);}
747:     MPI_Group_free(&group);
748:     PetscFree(dgroupranks);
749:   }

751:   /* Partition ranks[] into distinguished (first sf->ndranks) followed by non-distinguished */
752:   for (sf->ndranks=0,i=sf->nranks; sf->ndranks<i; ) {
753:     for (i--; sf->ndranks<i; i--) { /* Scan i backward looking for distinguished rank */
754:       if (InList(ranks[i],groupsize,groupranks)) break;
755:     }
756:     for ( ; sf->ndranks<=i; sf->ndranks++) { /* Scan sf->ndranks forward looking for non-distinguished rank */
757:       if (!InList(ranks[sf->ndranks],groupsize,groupranks)) break;
758:     }
759:     if (sf->ndranks < i) {                         /* Swap ranks[sf->ndranks] with ranks[i] */
760:       PetscInt    tmprank,tmpcount;

762:       tmprank             = ranks[i];
763:       tmpcount            = rcount[i];
764:       ranks[i]            = ranks[sf->ndranks];
765:       rcount[i]           = rcount[sf->ndranks];
766:       ranks[sf->ndranks]  = tmprank;
767:       rcount[sf->ndranks] = tmpcount;
768:       sf->ndranks++;
769:     }
770:   }
771:   PetscFree(groupranks);
772:   PetscSortIntWithArray(sf->ndranks,ranks,rcount);
773:   PetscSortIntWithArray(sf->nranks-sf->ndranks,ranks+sf->ndranks,rcount+sf->ndranks);
774:   sf->roffset[0] = 0;
775:   for (i=0; i<sf->nranks; i++) {
776:     PetscMPIIntCast(ranks[i],sf->ranks+i);
777:     sf->roffset[i+1] = sf->roffset[i] + rcount[i];
778:     rcount[i]        = 0;
779:   }
780:   for (i=0, irank = -1, orank = -1; i<sf->nleaves; i++) {
781:     /* short circuit */
782:     if (orank != sf->remote[i].rank) {
783:       /* Search for index of iremote[i].rank in sf->ranks */
784:       PetscFindMPIInt(sf->remote[i].rank,sf->ndranks,sf->ranks,&irank);
785:       if (irank < 0) {
786:         PetscFindMPIInt(sf->remote[i].rank,sf->nranks-sf->ndranks,sf->ranks+sf->ndranks,&irank);
787:         if (irank >= 0) irank += sf->ndranks;
788:       }
789:       orank = sf->remote[i].rank;
790:     }
791:     if (irank < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Could not find rank %D in array",sf->remote[i].rank);
792:     sf->rmine[sf->roffset[irank] + rcount[irank]]   = sf->mine ? sf->mine[i] : i;
793:     sf->rremote[sf->roffset[irank] + rcount[irank]] = sf->remote[i].index;
794:     rcount[irank]++;
795:   }
796:   PetscFree2(rcount,ranks);
797:   return(0);
798: }

800: /*@C
801:    PetscSFGetGroups - gets incoming and outgoing process groups

803:    Collective

805:    Input Argument:
806: .  sf - star forest

808:    Output Arguments:
809: +  incoming - group of origin processes for incoming edges (leaves that reference my roots)
810: -  outgoing - group of destination processes for outgoing edges (roots that I reference)

812:    Level: developer

814: .seealso: PetscSFGetWindow(), PetscSFRestoreWindow()
815: @*/
816: PetscErrorCode PetscSFGetGroups(PetscSF sf,MPI_Group *incoming,MPI_Group *outgoing)
817: {
819:   MPI_Group      group = MPI_GROUP_NULL;

822:   if (!sf->setupcalled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call PetscSFSetUp() before obtaining groups");
823:   if (sf->ingroup == MPI_GROUP_NULL) {
824:     PetscInt       i;
825:     const PetscInt *indegree;
826:     PetscMPIInt    rank,*outranks,*inranks;
827:     PetscSFNode    *remote;
828:     PetscSF        bgcount;

830:     /* Compute the number of incoming ranks */
831:     PetscMalloc1(sf->nranks,&remote);
832:     for (i=0; i<sf->nranks; i++) {
833:       remote[i].rank  = sf->ranks[i];
834:       remote[i].index = 0;
835:     }
836:     PetscSFDuplicate(sf,PETSCSF_DUPLICATE_CONFONLY,&bgcount);
837:     PetscSFSetGraph(bgcount,1,sf->nranks,NULL,PETSC_COPY_VALUES,remote,PETSC_OWN_POINTER);
838:     PetscSFComputeDegreeBegin(bgcount,&indegree);
839:     PetscSFComputeDegreeEnd(bgcount,&indegree);

841:     /* Enumerate the incoming ranks */
842:     PetscMalloc2(indegree[0],&inranks,sf->nranks,&outranks);
843:     MPI_Comm_rank(PetscObjectComm((PetscObject)sf),&rank);
844:     for (i=0; i<sf->nranks; i++) outranks[i] = rank;
845:     PetscSFGatherBegin(bgcount,MPI_INT,outranks,inranks);
846:     PetscSFGatherEnd(bgcount,MPI_INT,outranks,inranks);
847:     MPI_Comm_group(PetscObjectComm((PetscObject)sf),&group);
848:     MPI_Group_incl(group,indegree[0],inranks,&sf->ingroup);
849:     MPI_Group_free(&group);
850:     PetscFree2(inranks,outranks);
851:     PetscSFDestroy(&bgcount);
852:   }
853:   *incoming = sf->ingroup;

855:   if (sf->outgroup == MPI_GROUP_NULL) {
856:     MPI_Comm_group(PetscObjectComm((PetscObject)sf),&group);
857:     MPI_Group_incl(group,sf->nranks,sf->ranks,&sf->outgroup);
858:     MPI_Group_free(&group);
859:   }
860:   *outgoing = sf->outgroup;
861:   return(0);
862: }

864: /*@
865:    PetscSFGetMultiSF - gets the inner SF implemeting gathers and scatters

867:    Collective

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

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

875:    Level: developer

877:    Notes:

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

883: .seealso: PetscSFSetGraph(), PetscSFGatherBegin(), PetscSFScatterBegin(), PetscSFComputeMultiRootOriginalNumbering()
884: @*/
885: PetscErrorCode PetscSFGetMultiSF(PetscSF sf,PetscSF *multi)
886: {

892:   if (sf->nroots < 0) {         /* Graph has not been set yet; why do we need this? */
893:     PetscSFDuplicate(sf,PETSCSF_DUPLICATE_RANKS,&sf->multi);
894:     *multi = sf->multi;
895:     return(0);
896:   }
897:   if (!sf->multi) {
898:     const PetscInt *indegree;
899:     PetscInt       i,*inoffset,*outones,*outoffset,maxlocal;
900:     PetscSFNode    *remote;
901:     maxlocal = sf->maxleaf+1; /* TODO: We should use PetscSFGetLeafRange() */
902:     PetscSFComputeDegreeBegin(sf,&indegree);
903:     PetscSFComputeDegreeEnd(sf,&indegree);
904:     PetscMalloc3(sf->nroots+1,&inoffset,maxlocal,&outones,maxlocal,&outoffset);
905:     inoffset[0] = 0;
906:     for (i=0; i<sf->nroots; i++) inoffset[i+1] = inoffset[i] + indegree[i];
907:     for (i=0; i<maxlocal; i++) outones[i] = 1;
908:     PetscSFFetchAndOpBegin(sf,MPIU_INT,inoffset,outones,outoffset,MPI_SUM);
909:     PetscSFFetchAndOpEnd(sf,MPIU_INT,inoffset,outones,outoffset,MPI_SUM);
910:     for (i=0; i<sf->nroots; i++) inoffset[i] -= indegree[i]; /* Undo the increment */
911: #if 0
912: #if defined(PETSC_USE_DEBUG)                                 /* Check that the expected number of increments occurred */
913:     for (i=0; i<sf->nroots; i++) {
914:       if (inoffset[i] + indegree[i] != inoffset[i+1]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Incorrect result after PetscSFFetchAndOp");
915:     }
916: #endif
917: #endif
918:     PetscMalloc1(sf->nleaves,&remote);
919:     for (i=0; i<sf->nleaves; i++) {
920:       remote[i].rank  = sf->remote[i].rank;
921:       remote[i].index = outoffset[sf->mine ? sf->mine[i] : i];
922:     }
923:     PetscSFDuplicate(sf,PETSCSF_DUPLICATE_RANKS,&sf->multi);
924:     PetscSFSetGraph(sf->multi,inoffset[sf->nroots],sf->nleaves,sf->mine,PETSC_COPY_VALUES,remote,PETSC_OWN_POINTER);
925:     if (sf->rankorder) {        /* Sort the ranks */
926:       PetscMPIInt rank;
927:       PetscInt    *inranks,*newoffset,*outranks,*newoutoffset,*tmpoffset,maxdegree;
928:       PetscSFNode *newremote;
929:       MPI_Comm_rank(PetscObjectComm((PetscObject)sf),&rank);
930:       for (i=0,maxdegree=0; i<sf->nroots; i++) maxdegree = PetscMax(maxdegree,indegree[i]);
931:       PetscMalloc5(sf->multi->nroots,&inranks,sf->multi->nroots,&newoffset,maxlocal,&outranks,maxlocal,&newoutoffset,maxdegree,&tmpoffset);
932:       for (i=0; i<maxlocal; i++) outranks[i] = rank;
933:       PetscSFReduceBegin(sf->multi,MPIU_INT,outranks,inranks,MPIU_REPLACE);
934:       PetscSFReduceEnd(sf->multi,MPIU_INT,outranks,inranks,MPIU_REPLACE);
935:       /* Sort the incoming ranks at each vertex, build the inverse map */
936:       for (i=0; i<sf->nroots; i++) {
937:         PetscInt j;
938:         for (j=0; j<indegree[i]; j++) tmpoffset[j] = j;
939:         PetscSortIntWithArray(indegree[i],inranks+inoffset[i],tmpoffset);
940:         for (j=0; j<indegree[i]; j++) newoffset[inoffset[i] + tmpoffset[j]] = inoffset[i] + j;
941:       }
942:       PetscSFBcastBegin(sf->multi,MPIU_INT,newoffset,newoutoffset);
943:       PetscSFBcastEnd(sf->multi,MPIU_INT,newoffset,newoutoffset);
944:       PetscMalloc1(sf->nleaves,&newremote);
945:       for (i=0; i<sf->nleaves; i++) {
946:         newremote[i].rank  = sf->remote[i].rank;
947:         newremote[i].index = newoutoffset[sf->mine ? sf->mine[i] : i];
948:       }
949:       PetscSFSetGraph(sf->multi,inoffset[sf->nroots],sf->nleaves,sf->mine,PETSC_COPY_VALUES,newremote,PETSC_OWN_POINTER);
950:       PetscFree5(inranks,newoffset,outranks,newoutoffset,tmpoffset);
951:     }
952:     PetscFree3(inoffset,outones,outoffset);
953:   }
954:   *multi = sf->multi;
955:   return(0);
956: }

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

961:    Collective

963:    Input Arguments:
964: +  sf - original star forest
965: .  nselected  - number of selected roots on this process
966: -  selected   - indices of the selected roots on this process

968:    Output Arguments:
969: .  newsf - new star forest

971:    Level: advanced

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

977: .seealso: PetscSFSetGraph(), PetscSFGetGraph()
978: @*/
979: PetscErrorCode PetscSFCreateEmbeddedSF(PetscSF sf,PetscInt nselected,const PetscInt *selected,PetscSF *newsf)
980: {
981:   PetscInt          *rootdata,*leafdata,*new_ilocal;
982:   PetscSFNode       *new_iremote;
983:   const PetscInt    *ilocal;
984:   const PetscSFNode *iremote;
985:   PetscInt          nleaves,nroots,n,i,new_nleaves = 0;
986:   PetscSF           tmpsf;
987:   PetscErrorCode    ierr;

991:   PetscSFCheckGraphSet(sf,1);
994:   PetscLogEventBegin(PETSCSF_EmbedSF,sf,0,0,0);

996:   /* Find out which leaves (not leaf data items) are still connected to roots in the embedded sf */
997:   PetscSFGetGraph(sf,&nroots,&nleaves,&ilocal,&iremote);
998:   PetscSFDuplicate(sf,PETSCSF_DUPLICATE_RANKS,&tmpsf);
999:   PetscSFSetGraph(tmpsf,nroots,nleaves,NULL/*contiguous*/,PETSC_USE_POINTER,iremote,PETSC_USE_POINTER);
1000:   PetscCalloc2(nroots,&rootdata,nleaves,&leafdata);
1001:   for (i=0; i<nselected; ++i) {
1002:     if (selected[i] < 0 || selected[i] >= nroots) SETERRQ2(PetscObjectComm((PetscObject)sf),PETSC_ERR_ARG_OUTOFRANGE,"Root index %D is not in [0,%D)",selected[i],nroots);
1003:     rootdata[selected[i]] = 1;
1004:   }

1006:   PetscSFBcastBegin(tmpsf,MPIU_INT,rootdata,leafdata);
1007:   PetscSFBcastEnd(tmpsf,MPIU_INT,rootdata,leafdata);
1008:   PetscSFDestroy(&tmpsf);

1010:   /* Build newsf with leaves that are still connected */
1011:   for (i = 0; i < nleaves; ++i) new_nleaves += leafdata[i];
1012:   PetscMalloc1(new_nleaves,&new_ilocal);
1013:   PetscMalloc1(new_nleaves,&new_iremote);
1014:   for (i = 0, n = 0; i < nleaves; ++i) {
1015:     if (leafdata[i]) {
1016:       new_ilocal[n]        = sf->mine ? sf->mine[i] : i;
1017:       new_iremote[n].rank  = sf->remote[i].rank;
1018:       new_iremote[n].index = sf->remote[i].index;
1019:       ++n;
1020:     }
1021:   }
1022:   if (n != new_nleaves) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"There is a size mismatch in the SF embedding, %D != %D",n,new_nleaves);
1023:   PetscSFDuplicate(sf,PETSCSF_DUPLICATE_CONFONLY,newsf);
1024:   PetscSFSetGraph(*newsf,nroots,new_nleaves,new_ilocal,PETSC_OWN_POINTER,new_iremote,PETSC_OWN_POINTER);
1025:   PetscFree2(rootdata,leafdata);
1026:   PetscSFSetUp(*newsf);
1027:   PetscLogEventEnd(PETSCSF_EmbedSF,sf,0,0,0);
1028:   return(0);
1029: }

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

1034:   Collective

1036:   Input Arguments:
1037: + sf - original star forest
1038: . nleaves - number of leaves to select on this process
1039: - selected - selected leaves on this process

1041:   Output Arguments:
1042: .  newsf - new star forest

1044:   Level: advanced

1046: .seealso: PetscSFCreateEmbeddedSF(), PetscSFSetGraph(), PetscSFGetGraph()
1047: @*/
1048: PetscErrorCode PetscSFCreateEmbeddedLeafSF(PetscSF sf, PetscInt nleaves, const PetscInt *selected, PetscSF *newsf)
1049: {
1050:   PetscSFNode   *iremote;
1051:   PetscInt      *ilocal;
1052:   PetscInt       i;

1057:   PetscSFCheckGraphSet(sf, 1);
1060:   PetscMalloc1(nleaves, &ilocal);
1061:   PetscMalloc1(nleaves, &iremote);
1062:   for (i = 0; i < nleaves; ++i) {
1063:     const PetscInt l = selected[i];

1065:     ilocal[i]        = sf->mine ? sf->mine[l] : l;
1066:     iremote[i].rank  = sf->remote[l].rank;
1067:     iremote[i].index = sf->remote[l].index;
1068:   }
1069:   PetscSFDuplicate(sf, PETSCSF_DUPLICATE_RANKS, newsf);
1070:   PetscSFSetGraph(*newsf, sf->nroots, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
1071:   return(0);
1072: }

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

1077:    Collective on PetscSF

1079:    Input Arguments:
1080: +  sf - star forest on which to communicate
1081: .  unit - data type associated with each node
1082: -  rootdata - buffer to broadcast

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

1087:    Level: intermediate

1089: .seealso: PetscSFCreate(), PetscSFSetGraph(), PetscSFView(), PetscSFBcastEnd(), PetscSFReduceBegin()
1090: @*/
1091: PetscErrorCode PetscSFBcastBegin(PetscSF sf,MPI_Datatype unit,const void *rootdata,void *leafdata)
1092: {

1097:   PetscSFSetUp(sf);
1098:   PetscLogEventBegin(PETSCSF_BcastBegin,sf,0,0,0);
1099:   (*sf->ops->BcastBegin)(sf,unit,rootdata,leafdata);
1100:   PetscLogEventEnd(PETSCSF_BcastBegin,sf,0,0,0);
1101:   return(0);
1102: }

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

1107:    Collective

1109:    Input Arguments:
1110: +  sf - star forest
1111: .  unit - data type
1112: -  rootdata - buffer to broadcast

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

1117:    Level: intermediate

1119: .seealso: PetscSFSetGraph(), PetscSFReduceEnd()
1120: @*/
1121: PetscErrorCode PetscSFBcastEnd(PetscSF sf,MPI_Datatype unit,const void *rootdata,void *leafdata)
1122: {

1127:   PetscSFSetUp(sf);
1128:   PetscLogEventBegin(PETSCSF_BcastEnd,sf,0,0,0);
1129:   (*sf->ops->BcastEnd)(sf,unit,rootdata,leafdata);
1130:   PetscLogEventEnd(PETSCSF_BcastEnd,sf,0,0,0);
1131:   return(0);
1132: }

1134: /*@C
1135:    PetscSFBcastAndOpBegin - begin pointwise broadcast with root value being reduced to leaf value, to be concluded with call to PetscSFBcastAndOpEnd()

1137:    Collective on PetscSF

1139:    Input Arguments:
1140: +  sf - star forest on which to communicate
1141: .  unit - data type associated with each node
1142: .  rootdata - buffer to broadcast
1143: -  op - operation to use for reduction

1145:    Output Arguments:
1146: .  leafdata - buffer to be reduced with values from each leaf's respective root

1148:    Level: intermediate

1150: .seealso: PetscSFBcastAndOpEnd(), PetscSFBcastBegin(), PetscSFBcastEnd()
1151: @*/
1152: PetscErrorCode PetscSFBcastAndOpBegin(PetscSF sf,MPI_Datatype unit,const void *rootdata,void *leafdata,MPI_Op op)
1153: {

1158:   PetscSFSetUp(sf);
1159:   PetscLogEventBegin(PETSCSF_BcastAndOpBegin,sf,0,0,0);
1160:   (*sf->ops->BcastAndOpBegin)(sf,unit,rootdata,leafdata,op);
1161:   PetscLogEventEnd(PETSCSF_BcastAndOpBegin,sf,0,0,0);
1162:   return(0);
1163: }

1165: /*@C
1166:    PetscSFBcastAndOpEnd - end a broadcast & reduce operation started with PetscSFBcastAndOpBegin()

1168:    Collective

1170:    Input Arguments:
1171: +  sf - star forest
1172: .  unit - data type
1173: .  rootdata - buffer to broadcast
1174: -  op - operation to use for reduction

1176:    Output Arguments:
1177: .  leafdata - buffer to be reduced with values from each leaf's respective root

1179:    Level: intermediate

1181: .seealso: PetscSFSetGraph(), PetscSFReduceEnd()
1182: @*/
1183: PetscErrorCode PetscSFBcastAndOpEnd(PetscSF sf,MPI_Datatype unit,const void *rootdata,void *leafdata,MPI_Op op)
1184: {

1189:   PetscSFSetUp(sf);
1190:   PetscLogEventBegin(PETSCSF_BcastAndOpEnd,sf,0,0,0);
1191:   (*sf->ops->BcastAndOpEnd)(sf,unit,rootdata,leafdata,op);
1192:   PetscLogEventEnd(PETSCSF_BcastAndOpEnd,sf,0,0,0);
1193:   return(0);
1194: }

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

1199:    Collective

1201:    Input Arguments:
1202: +  sf - star forest
1203: .  unit - data type
1204: .  leafdata - values to reduce
1205: -  op - reduction operation

1207:    Output Arguments:
1208: .  rootdata - result of reduction of values from all leaves of each root

1210:    Level: intermediate

1212: .seealso: PetscSFBcastBegin()
1213: @*/
1214: PetscErrorCode PetscSFReduceBegin(PetscSF sf,MPI_Datatype unit,const void *leafdata,void *rootdata,MPI_Op op)
1215: {

1220:   PetscSFSetUp(sf);
1221:   PetscLogEventBegin(PETSCSF_ReduceBegin,sf,0,0,0);
1222:   (sf->ops->ReduceBegin)(sf,unit,leafdata,rootdata,op);
1223:   PetscLogEventEnd(PETSCSF_ReduceBegin,sf,0,0,0);
1224:   return(0);
1225: }

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

1230:    Collective

1232:    Input Arguments:
1233: +  sf - star forest
1234: .  unit - data type
1235: .  leafdata - values to reduce
1236: -  op - reduction operation

1238:    Output Arguments:
1239: .  rootdata - result of reduction of values from all leaves of each root

1241:    Level: intermediate

1243: .seealso: PetscSFSetGraph(), PetscSFBcastEnd()
1244: @*/
1245: PetscErrorCode PetscSFReduceEnd(PetscSF sf,MPI_Datatype unit,const void *leafdata,void *rootdata,MPI_Op op)
1246: {

1251:   PetscSFSetUp(sf);
1252:   PetscLogEventBegin(PETSCSF_ReduceEnd,sf,0,0,0);
1253:   (*sf->ops->ReduceEnd)(sf,unit,leafdata,rootdata,op);
1254:   PetscLogEventEnd(PETSCSF_ReduceEnd,sf,0,0,0);
1255:   return(0);
1256: }

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

1261:    Collective

1263:    Input Arguments:
1264: .  sf - star forest

1266:    Output Arguments:
1267: .  degree - degree of each root vertex

1269:    Level: advanced

1271:    Notes:
1272:    The returned array is owned by PetscSF and automatically freed by PetscSFDestroy(). Hence no need to call PetscFree() on it.

1274: .seealso: PetscSFGatherBegin()
1275: @*/
1276: PetscErrorCode PetscSFComputeDegreeBegin(PetscSF sf,const PetscInt **degree)
1277: {

1282:   PetscSFCheckGraphSet(sf,1);
1284:   if (!sf->degreeknown) {
1285:     PetscInt i, nroots = sf->nroots, maxlocal = sf->maxleaf+1;  /* TODO: We should use PetscSFGetLeafRange() */
1286:     if (sf->degree) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Calls to PetscSFComputeDegreeBegin() cannot be nested.");
1287:     PetscMalloc1(nroots,&sf->degree);
1288:     PetscMalloc1(PetscMax(maxlocal,1),&sf->degreetmp); /* allocate at least one entry, see check in PetscSFComputeDegreeEnd() */
1289:     for (i=0; i<nroots; i++) sf->degree[i] = 0;
1290:     for (i=0; i<maxlocal; i++) sf->degreetmp[i] = 1;
1291:     PetscSFReduceBegin(sf,MPIU_INT,sf->degreetmp,sf->degree,MPI_SUM);
1292:   }
1293:   *degree = NULL;
1294:   return(0);
1295: }

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

1300:    Collective

1302:    Input Arguments:
1303: .  sf - star forest

1305:    Output Arguments:
1306: .  degree - degree of each root vertex

1308:    Level: developer

1310:    Notes:
1311:    The returned array is owned by PetscSF and automatically freed by PetscSFDestroy(). Hence no need to call PetscFree() on it.

1313: .seealso:
1314: @*/
1315: PetscErrorCode PetscSFComputeDegreeEnd(PetscSF sf,const PetscInt **degree)
1316: {

1321:   PetscSFCheckGraphSet(sf,1);
1323:   if (!sf->degreeknown) {
1324:     if (!sf->degreetmp) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PetscSFComputeDegreeBegin() before PetscSFComputeDegreeEnd()");
1325:     PetscSFReduceEnd(sf,MPIU_INT,sf->degreetmp,sf->degree,MPI_SUM);
1326:     PetscFree(sf->degreetmp);
1327:     sf->degreeknown = PETSC_TRUE;
1328:   }
1329:   *degree = sf->degree;
1330:   return(0);
1331: }


1334: /*@C
1335:    PetscSFComputeMultiRootOriginalNumbering - Returns original numbering of multi-roots (roots of multi-SF returned by PetscSFGetMultiSF()).
1336:    Each multi-root is assigned index of the corresponding original root.

1338:    Collective

1340:    Input Arguments:
1341: +  sf - star forest
1342: -  degree - degree of each root vertex, computed with PetscSFComputeDegreeBegin()/PetscSFComputeDegreeEnd()

1344:    Output Arguments:
1345: +  nMultiRoots - (optional) number of multi-roots (roots of multi-SF)
1346: -  multiRootsOrigNumbering - original indices of multi-roots; length of this array is nMultiRoots

1348:    Level: developer
1349:    
1350:    Notes:
1351:    The returned array multiRootsOrigNumbering is newly allocated and should be destroyed with PetscFree() when no longer needed.

1353: .seealso: PetscSFComputeDegreeBegin(), PetscSFComputeDegreeEnd(), PetscSFGetMultiSF()
1354: @*/
1355: PetscErrorCode PetscSFComputeMultiRootOriginalNumbering(PetscSF sf, const PetscInt degree[], PetscInt *nMultiRoots, PetscInt *multiRootsOrigNumbering[])
1356: {
1357:   PetscSF             msf;
1358:   PetscInt            i, j, k, nroots, nmroots;
1359:   PetscErrorCode      ierr;

1363:   PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
1367:   PetscSFGetMultiSF(sf,&msf);
1368:   PetscSFGetGraph(msf, &nmroots, NULL, NULL, NULL);
1369:   PetscMalloc1(nmroots, multiRootsOrigNumbering);
1370:   for (i=0,j=0,k=0; i<nroots; i++) {
1371:     if (!degree[i]) continue;
1372:     for (j=0; j<degree[i]; j++,k++) {
1373:       (*multiRootsOrigNumbering)[k] = i;
1374:     }
1375:   }
1376: #if defined(PETSC_USE_DEBUG)
1377:   if (PetscUnlikely(k != nmroots)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"sanity check fail");
1378: #endif
1379:   if (nMultiRoots) *nMultiRoots = nmroots;
1380:   return(0);
1381: }

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

1386:    Collective

1388:    Input Arguments:
1389: +  sf - star forest
1390: .  unit - data type
1391: .  leafdata - leaf values to use in reduction
1392: -  op - operation to use for reduction

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

1398:    Level: advanced

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

1406: .seealso: PetscSFComputeDegreeBegin(), PetscSFReduceBegin(), PetscSFSetGraph()
1407: @*/
1408: PetscErrorCode PetscSFFetchAndOpBegin(PetscSF sf,MPI_Datatype unit,void *rootdata,const void *leafdata,void *leafupdate,MPI_Op op)
1409: {

1414:   PetscSFSetUp(sf);
1415:   PetscLogEventBegin(PETSCSF_FetchAndOpBegin,sf,0,0,0);
1416:   (*sf->ops->FetchAndOpBegin)(sf,unit,rootdata,leafdata,leafupdate,op);
1417:   PetscLogEventEnd(PETSCSF_FetchAndOpBegin,sf,0,0,0);
1418:   return(0);
1419: }

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

1424:    Collective

1426:    Input Arguments:
1427: +  sf - star forest
1428: .  unit - data type
1429: .  leafdata - leaf values to use in reduction
1430: -  op - operation to use for reduction

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

1436:    Level: advanced

1438: .seealso: PetscSFComputeDegreeEnd(), PetscSFReduceEnd(), PetscSFSetGraph()
1439: @*/
1440: PetscErrorCode PetscSFFetchAndOpEnd(PetscSF sf,MPI_Datatype unit,void *rootdata,const void *leafdata,void *leafupdate,MPI_Op op)
1441: {

1446:   PetscSFSetUp(sf);
1447:   PetscLogEventBegin(PETSCSF_FetchAndOpEnd,sf,0,0,0);
1448:   (*sf->ops->FetchAndOpEnd)(sf,unit,rootdata,leafdata,leafupdate,op);
1449:   PetscLogEventEnd(PETSCSF_FetchAndOpEnd,sf,0,0,0);
1450:   return(0);
1451: }

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

1456:    Collective

1458:    Input Arguments:
1459: +  sf - star forest
1460: .  unit - data type
1461: -  leafdata - leaf data to gather to roots

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

1466:    Level: intermediate

1468: .seealso: PetscSFComputeDegreeBegin(), PetscSFScatterBegin()
1469: @*/
1470: PetscErrorCode PetscSFGatherBegin(PetscSF sf,MPI_Datatype unit,const void *leafdata,void *multirootdata)
1471: {
1473:   PetscSF        multi;

1477:   PetscSFSetUp(sf);
1478:   PetscSFGetMultiSF(sf,&multi);
1479:   PetscSFReduceBegin(multi,unit,leafdata,multirootdata,MPIU_REPLACE);
1480:   return(0);
1481: }

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

1486:    Collective

1488:    Input Arguments:
1489: +  sf - star forest
1490: .  unit - data type
1491: -  leafdata - leaf data to gather to roots

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

1496:    Level: intermediate

1498: .seealso: PetscSFComputeDegreeEnd(), PetscSFScatterEnd()
1499: @*/
1500: PetscErrorCode PetscSFGatherEnd(PetscSF sf,MPI_Datatype unit,const void *leafdata,void *multirootdata)
1501: {
1503:   PetscSF        multi;

1507:   PetscSFSetUp(sf);
1508:   PetscSFGetMultiSF(sf,&multi);
1509:   PetscSFReduceEnd(multi,unit,leafdata,multirootdata,MPIU_REPLACE);
1510:   return(0);
1511: }

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

1516:    Collective

1518:    Input Arguments:
1519: +  sf - star forest
1520: .  unit - data type
1521: -  multirootdata - root buffer to send to each leaf, one unit of data per leaf

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

1526:    Level: intermediate

1528: .seealso: PetscSFComputeDegreeBegin(), PetscSFScatterBegin()
1529: @*/
1530: PetscErrorCode PetscSFScatterBegin(PetscSF sf,MPI_Datatype unit,const void *multirootdata,void *leafdata)
1531: {
1533:   PetscSF        multi;

1537:   PetscSFSetUp(sf);
1538:   PetscSFGetMultiSF(sf,&multi);
1539:   PetscSFBcastBegin(multi,unit,multirootdata,leafdata);
1540:   return(0);
1541: }

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

1546:    Collective

1548:    Input Arguments:
1549: +  sf - star forest
1550: .  unit - data type
1551: -  multirootdata - root buffer to send to each leaf, one unit of data per leaf

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

1556:    Level: intermediate

1558: .seealso: PetscSFComputeDegreeEnd(), PetscSFScatterEnd()
1559: @*/
1560: PetscErrorCode PetscSFScatterEnd(PetscSF sf,MPI_Datatype unit,const void *multirootdata,void *leafdata)
1561: {
1563:   PetscSF        multi;

1567:   PetscSFSetUp(sf);
1568:   PetscSFGetMultiSF(sf,&multi);
1569:   PetscSFBcastEnd(multi,unit,multirootdata,leafdata);
1570:   return(0);
1571: }

1573: /*@
1574:   PetscSFCompose - Compose a new PetscSF equivalent to action to PetscSFs

1576:   Input Parameters:
1577: + sfA - The first PetscSF
1578: - sfB - The second PetscSF

1580:   Output Parameters:
1581: . sfBA - equvalent PetscSF for applying A then B

1583:   Level: developer

1585: .seealso: PetscSF, PetscSFGetGraph(), PetscSFSetGraph()
1586: @*/
1587: PetscErrorCode PetscSFCompose(PetscSF sfA, PetscSF sfB, PetscSF *sfBA)
1588: {
1589:   MPI_Comm           comm;
1590:   const PetscSFNode *remotePointsA, *remotePointsB;
1591:   PetscSFNode       *remotePointsBA;
1592:   const PetscInt    *localPointsA, *localPointsB;
1593:   PetscInt           numRootsA, numLeavesA, numRootsB, numLeavesB;
1594:   PetscErrorCode     ierr;

1598:   PetscSFCheckGraphSet(sfA, 1);
1600:   PetscSFCheckGraphSet(sfB, 2);
1602:   PetscObjectGetComm((PetscObject) sfA, &comm);
1603:   PetscSFGetGraph(sfA, &numRootsA, &numLeavesA, &localPointsA, &remotePointsA);
1604:   PetscSFGetGraph(sfB, &numRootsB, &numLeavesB, &localPointsB, &remotePointsB);
1605:   PetscMalloc1(numLeavesB, &remotePointsBA);
1606:   PetscSFBcastBegin(sfB, MPIU_2INT, remotePointsA, remotePointsBA);
1607:   PetscSFBcastEnd(sfB, MPIU_2INT, remotePointsA, remotePointsBA);
1608:   PetscSFCreate(comm, sfBA);
1609:   PetscSFSetGraph(*sfBA, numRootsA, numLeavesB, localPointsB, PETSC_COPY_VALUES, remotePointsBA, PETSC_OWN_POINTER);
1610:   return(0);
1611: }