Actual source code: sf.c

petsc-master 2019-05-25
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 on MPI_Comm

 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:
297:     In Fortran you must use PETSC_COPY_VALUES for localmode and remotemode

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

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

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

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

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

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

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

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

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

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

389:    Collective

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

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

397:    Level: advanced

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

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

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

416:   PetscSFCheckGraphSet(sf,1);

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

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

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

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

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

459:    Collective

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

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

468:    Level: beginner

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

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

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

499:    Not Collective

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

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

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

513:    Level: intermediate

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

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

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

532:    Not Collective

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

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

541:    Level: developer

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

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

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

559:    Collective

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

565:    Level: beginner

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

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

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

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

628:    Not Collective

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

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

640:    Level: developer

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

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

658: /*@C
659:    PetscSFGetLeafRanks - Get leaf ranks referencing roots on this process

661:    Not Collective

663:    Input Arguments:
664: .  sf - star forest

666:    Output Arguments:
667: +  niranks - number of leaf ranks referencing roots on this process
668: .  iranks - array of ranks
669: .  ioffset - offset in irootloc for each rank (length niranks+1)
670: -  irootloc - concatenated array holding local indices of roots referenced by each leaf rank

672:    Level: developer

674: .seealso: PetscSFGetRanks()
675: @*/
676: PetscErrorCode PetscSFGetLeafRanks(PetscSF sf,PetscInt *niranks,const PetscMPIInt **iranks,const PetscInt **ioffset,const PetscInt **irootloc)
677: {

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

693: static PetscBool InList(PetscMPIInt needle,PetscMPIInt n,const PetscMPIInt *list) {
694:   PetscInt i;
695:   for (i=0; i<n; i++) {
696:     if (needle == list[i]) return PETSC_TRUE;
697:   }
698:   return PETSC_FALSE;
699: }

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

704:    Collective

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

710:    Level: developer

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

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

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

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

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

805: /*@C
806:    PetscSFGetGroups - gets incoming and outgoing process groups

808:    Collective

810:    Input Argument:
811: .  sf - star forest

813:    Output Arguments:
814: +  incoming - group of origin processes for incoming edges (leaves that reference my roots)
815: -  outgoing - group of destination processes for outgoing edges (roots that I reference)

817:    Level: developer

819: .seealso: PetscSFGetWindow(), PetscSFRestoreWindow()
820: @*/
821: PetscErrorCode PetscSFGetGroups(PetscSF sf,MPI_Group *incoming,MPI_Group *outgoing)
822: {
824:   MPI_Group      group = MPI_GROUP_NULL;

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

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

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

860:   if (sf->outgroup == MPI_GROUP_NULL) {
861:     MPI_Comm_group(PetscObjectComm((PetscObject)sf),&group);
862:     MPI_Group_incl(group,sf->nranks,sf->ranks,&sf->outgroup);
863:     MPI_Group_free(&group);
864:   }
865:   *outgoing = sf->outgroup;
866:   return(0);
867: }

869: /*@
870:    PetscSFGetMultiSF - gets the inner SF implemeting gathers and scatters

872:    Collective

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

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

880:    Level: developer

882:    Notes:

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

888: .seealso: PetscSFSetGraph(), PetscSFGatherBegin(), PetscSFScatterBegin(), PetscSFComputeMultiRootOriginalNumbering()
889: @*/
890: PetscErrorCode PetscSFGetMultiSF(PetscSF sf,PetscSF *multi)
891: {

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

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

966:    Collective

968:    Input Arguments:
969: +  sf - original star forest
970: .  nselected  - number of selected roots on this process
971: -  selected   - indices of the selected roots on this process

973:    Output Arguments:
974: .  newsf - new star forest

976:    Level: advanced

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

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

996:   PetscSFCheckGraphSet(sf,1);
999:   PetscLogEventBegin(PETSCSF_EmbedSF,sf,0,0,0);

1001:   /* Find out which leaves (not leaf data items) are still connected to roots in the embedded sf */
1002:   PetscSFGetGraph(sf,&nroots,&nleaves,&ilocal,&iremote);
1003:   PetscSFDuplicate(sf,PETSCSF_DUPLICATE_RANKS,&tmpsf);
1004:   PetscSFSetGraph(tmpsf,nroots,nleaves,NULL/*contiguous*/,PETSC_USE_POINTER,iremote,PETSC_USE_POINTER);
1005:   PetscCalloc2(nroots,&rootdata,nleaves,&leafdata);
1006:   for (i=0; i<nselected; ++i) {
1007:     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);
1008:     rootdata[selected[i]] = 1;
1009:   }

1011:   PetscSFBcastBegin(tmpsf,MPIU_INT,rootdata,leafdata);
1012:   PetscSFBcastEnd(tmpsf,MPIU_INT,rootdata,leafdata);
1013:   PetscSFDestroy(&tmpsf);

1015:   /* Build newsf with leaves that are still connected */
1016:   for (i = 0; i < nleaves; ++i) new_nleaves += leafdata[i];
1017:   PetscMalloc1(new_nleaves,&new_ilocal);
1018:   PetscMalloc1(new_nleaves,&new_iremote);
1019:   for (i = 0, n = 0; i < nleaves; ++i) {
1020:     if (leafdata[i]) {
1021:       new_ilocal[n]        = sf->mine ? sf->mine[i] : i;
1022:       new_iremote[n].rank  = sf->remote[i].rank;
1023:       new_iremote[n].index = sf->remote[i].index;
1024:       ++n;
1025:     }
1026:   }
1027:   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);
1028:   PetscSFDuplicate(sf,PETSCSF_DUPLICATE_CONFONLY,newsf);
1029:   PetscSFSetGraph(*newsf,nroots,new_nleaves,new_ilocal,PETSC_OWN_POINTER,new_iremote,PETSC_OWN_POINTER);
1030:   PetscFree2(rootdata,leafdata);
1031:   PetscSFSetUp(*newsf);
1032:   PetscLogEventEnd(PETSCSF_EmbedSF,sf,0,0,0);
1033:   return(0);
1034: }

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

1039:   Collective

1041:   Input Arguments:
1042: + sf - original star forest
1043: . nleaves - number of leaves to select on this process
1044: - selected - selected leaves on this process

1046:   Output Arguments:
1047: .  newsf - new star forest

1049:   Level: advanced

1051: .seealso: PetscSFCreateEmbeddedSF(), PetscSFSetGraph(), PetscSFGetGraph()
1052: @*/
1053: PetscErrorCode PetscSFCreateEmbeddedLeafSF(PetscSF sf, PetscInt nleaves, const PetscInt *selected, PetscSF *newsf)
1054: {
1055:   PetscSFNode   *iremote;
1056:   PetscInt      *ilocal;
1057:   PetscInt       i;

1062:   PetscSFCheckGraphSet(sf, 1);
1065:   PetscMalloc1(nleaves, &ilocal);
1066:   PetscMalloc1(nleaves, &iremote);
1067:   for (i = 0; i < nleaves; ++i) {
1068:     const PetscInt l = selected[i];

1070:     ilocal[i]        = sf->mine ? sf->mine[l] : l;
1071:     iremote[i].rank  = sf->remote[l].rank;
1072:     iremote[i].index = sf->remote[l].index;
1073:   }
1074:   PetscSFDuplicate(sf, PETSCSF_DUPLICATE_RANKS, newsf);
1075:   PetscSFSetGraph(*newsf, sf->nroots, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
1076:   return(0);
1077: }

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

1082:    Collective on PetscSF

1084:    Input Arguments:
1085: +  sf - star forest on which to communicate
1086: .  unit - data type associated with each node
1087: -  rootdata - buffer to broadcast

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

1092:    Level: intermediate

1094: .seealso: PetscSFCreate(), PetscSFSetGraph(), PetscSFView(), PetscSFBcastEnd(), PetscSFReduceBegin()
1095: @*/
1096: PetscErrorCode PetscSFBcastBegin(PetscSF sf,MPI_Datatype unit,const void *rootdata,void *leafdata)
1097: {

1102:   PetscSFSetUp(sf);
1103:   PetscLogEventBegin(PETSCSF_BcastBegin,sf,0,0,0);
1104:   (*sf->ops->BcastBegin)(sf,unit,rootdata,leafdata);
1105:   PetscLogEventEnd(PETSCSF_BcastBegin,sf,0,0,0);
1106:   return(0);
1107: }

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

1112:    Collective

1114:    Input Arguments:
1115: +  sf - star forest
1116: .  unit - data type
1117: -  rootdata - buffer to broadcast

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

1122:    Level: intermediate

1124: .seealso: PetscSFSetGraph(), PetscSFReduceEnd()
1125: @*/
1126: PetscErrorCode PetscSFBcastEnd(PetscSF sf,MPI_Datatype unit,const void *rootdata,void *leafdata)
1127: {

1132:   PetscSFSetUp(sf);
1133:   PetscLogEventBegin(PETSCSF_BcastEnd,sf,0,0,0);
1134:   (*sf->ops->BcastEnd)(sf,unit,rootdata,leafdata);
1135:   PetscLogEventEnd(PETSCSF_BcastEnd,sf,0,0,0);
1136:   return(0);
1137: }

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

1142:    Collective on PetscSF

1144:    Input Arguments:
1145: +  sf - star forest on which to communicate
1146: .  unit - data type associated with each node
1147: .  rootdata - buffer to broadcast
1148: -  op - operation to use for reduction

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

1153:    Level: intermediate

1155: .seealso: PetscSFBcastAndOpEnd(), PetscSFBcastBegin(), PetscSFBcastEnd()
1156: @*/
1157: PetscErrorCode PetscSFBcastAndOpBegin(PetscSF sf,MPI_Datatype unit,const void *rootdata,void *leafdata,MPI_Op op)
1158: {

1163:   PetscSFSetUp(sf);
1164:   PetscLogEventBegin(PETSCSF_BcastAndOpBegin,sf,0,0,0);
1165:   (*sf->ops->BcastAndOpBegin)(sf,unit,rootdata,leafdata,op);
1166:   PetscLogEventEnd(PETSCSF_BcastAndOpBegin,sf,0,0,0);
1167:   return(0);
1168: }

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

1173:    Collective

1175:    Input Arguments:
1176: +  sf - star forest
1177: .  unit - data type
1178: .  rootdata - buffer to broadcast
1179: -  op - operation to use for reduction

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

1184:    Level: intermediate

1186: .seealso: PetscSFSetGraph(), PetscSFReduceEnd()
1187: @*/
1188: PetscErrorCode PetscSFBcastAndOpEnd(PetscSF sf,MPI_Datatype unit,const void *rootdata,void *leafdata,MPI_Op op)
1189: {

1194:   PetscSFSetUp(sf);
1195:   PetscLogEventBegin(PETSCSF_BcastAndOpEnd,sf,0,0,0);
1196:   (*sf->ops->BcastAndOpEnd)(sf,unit,rootdata,leafdata,op);
1197:   PetscLogEventEnd(PETSCSF_BcastAndOpEnd,sf,0,0,0);
1198:   return(0);
1199: }

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

1204:    Collective

1206:    Input Arguments:
1207: +  sf - star forest
1208: .  unit - data type
1209: .  leafdata - values to reduce
1210: -  op - reduction operation

1212:    Output Arguments:
1213: .  rootdata - result of reduction of values from all leaves of each root

1215:    Level: intermediate

1217: .seealso: PetscSFBcastBegin()
1218: @*/
1219: PetscErrorCode PetscSFReduceBegin(PetscSF sf,MPI_Datatype unit,const void *leafdata,void *rootdata,MPI_Op op)
1220: {

1225:   PetscSFSetUp(sf);
1226:   PetscLogEventBegin(PETSCSF_ReduceBegin,sf,0,0,0);
1227:   (sf->ops->ReduceBegin)(sf,unit,leafdata,rootdata,op);
1228:   PetscLogEventEnd(PETSCSF_ReduceBegin,sf,0,0,0);
1229:   return(0);
1230: }

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

1235:    Collective

1237:    Input Arguments:
1238: +  sf - star forest
1239: .  unit - data type
1240: .  leafdata - values to reduce
1241: -  op - reduction operation

1243:    Output Arguments:
1244: .  rootdata - result of reduction of values from all leaves of each root

1246:    Level: intermediate

1248: .seealso: PetscSFSetGraph(), PetscSFBcastEnd()
1249: @*/
1250: PetscErrorCode PetscSFReduceEnd(PetscSF sf,MPI_Datatype unit,const void *leafdata,void *rootdata,MPI_Op op)
1251: {

1256:   PetscSFSetUp(sf);
1257:   PetscLogEventBegin(PETSCSF_ReduceEnd,sf,0,0,0);
1258:   (*sf->ops->ReduceEnd)(sf,unit,leafdata,rootdata,op);
1259:   PetscLogEventEnd(PETSCSF_ReduceEnd,sf,0,0,0);
1260:   return(0);
1261: }

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

1266:    Collective

1268:    Input Arguments:
1269: .  sf - star forest

1271:    Output Arguments:
1272: .  degree - degree of each root vertex

1274:    Level: advanced

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

1279: .seealso: PetscSFGatherBegin()
1280: @*/
1281: PetscErrorCode PetscSFComputeDegreeBegin(PetscSF sf,const PetscInt **degree)
1282: {

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

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

1305:    Collective

1307:    Input Arguments:
1308: .  sf - star forest

1310:    Output Arguments:
1311: .  degree - degree of each root vertex

1313:    Level: developer

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

1318: .seealso:
1319: @*/
1320: PetscErrorCode PetscSFComputeDegreeEnd(PetscSF sf,const PetscInt **degree)
1321: {

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


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

1343:    Collective

1345:    Input Arguments:
1346: +  sf - star forest
1347: -  degree - degree of each root vertex, computed with PetscSFComputeDegreeBegin()/PetscSFComputeDegreeEnd()

1349:    Output Arguments:
1350: +  nMultiRoots - (optional) number of multi-roots (roots of multi-SF)
1351: -  multiRootsOrigNumbering - original indices of multi-roots; length of this array is nMultiRoots

1353:    Level: developer
1354:    
1355:    Notes:
1356:    The returned array multiRootsOrigNumbering is newly allocated and should be destroyed with PetscFree() when no longer needed.

1358: .seealso: PetscSFComputeDegreeBegin(), PetscSFComputeDegreeEnd(), PetscSFGetMultiSF()
1359: @*/
1360: PetscErrorCode PetscSFComputeMultiRootOriginalNumbering(PetscSF sf, const PetscInt degree[], PetscInt *nMultiRoots, PetscInt *multiRootsOrigNumbering[])
1361: {
1362:   PetscSF             msf;
1363:   PetscInt            i, j, k, nroots, nmroots;
1364:   PetscErrorCode      ierr;

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

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

1391:    Collective

1393:    Input Arguments:
1394: +  sf - star forest
1395: .  unit - data type
1396: .  leafdata - leaf values to use in reduction
1397: -  op - operation to use for reduction

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

1403:    Level: advanced

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

1411: .seealso: PetscSFComputeDegreeBegin(), PetscSFReduceBegin(), PetscSFSetGraph()
1412: @*/
1413: PetscErrorCode PetscSFFetchAndOpBegin(PetscSF sf,MPI_Datatype unit,void *rootdata,const void *leafdata,void *leafupdate,MPI_Op op)
1414: {

1419:   PetscSFSetUp(sf);
1420:   PetscLogEventBegin(PETSCSF_FetchAndOpBegin,sf,0,0,0);
1421:   (*sf->ops->FetchAndOpBegin)(sf,unit,rootdata,leafdata,leafupdate,op);
1422:   PetscLogEventEnd(PETSCSF_FetchAndOpBegin,sf,0,0,0);
1423:   return(0);
1424: }

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

1429:    Collective

1431:    Input Arguments:
1432: +  sf - star forest
1433: .  unit - data type
1434: .  leafdata - leaf values to use in reduction
1435: -  op - operation to use for reduction

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

1441:    Level: advanced

1443: .seealso: PetscSFComputeDegreeEnd(), PetscSFReduceEnd(), PetscSFSetGraph()
1444: @*/
1445: PetscErrorCode PetscSFFetchAndOpEnd(PetscSF sf,MPI_Datatype unit,void *rootdata,const void *leafdata,void *leafupdate,MPI_Op op)
1446: {

1451:   PetscSFSetUp(sf);
1452:   PetscLogEventBegin(PETSCSF_FetchAndOpEnd,sf,0,0,0);
1453:   (*sf->ops->FetchAndOpEnd)(sf,unit,rootdata,leafdata,leafupdate,op);
1454:   PetscLogEventEnd(PETSCSF_FetchAndOpEnd,sf,0,0,0);
1455:   return(0);
1456: }

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

1461:    Collective

1463:    Input Arguments:
1464: +  sf - star forest
1465: .  unit - data type
1466: -  leafdata - leaf data to gather to roots

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

1471:    Level: intermediate

1473: .seealso: PetscSFComputeDegreeBegin(), PetscSFScatterBegin()
1474: @*/
1475: PetscErrorCode PetscSFGatherBegin(PetscSF sf,MPI_Datatype unit,const void *leafdata,void *multirootdata)
1476: {
1478:   PetscSF        multi;

1482:   PetscSFSetUp(sf);
1483:   PetscSFGetMultiSF(sf,&multi);
1484:   PetscSFReduceBegin(multi,unit,leafdata,multirootdata,MPIU_REPLACE);
1485:   return(0);
1486: }

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

1491:    Collective

1493:    Input Arguments:
1494: +  sf - star forest
1495: .  unit - data type
1496: -  leafdata - leaf data to gather to roots

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

1501:    Level: intermediate

1503: .seealso: PetscSFComputeDegreeEnd(), PetscSFScatterEnd()
1504: @*/
1505: PetscErrorCode PetscSFGatherEnd(PetscSF sf,MPI_Datatype unit,const void *leafdata,void *multirootdata)
1506: {
1508:   PetscSF        multi;

1512:   PetscSFSetUp(sf);
1513:   PetscSFGetMultiSF(sf,&multi);
1514:   PetscSFReduceEnd(multi,unit,leafdata,multirootdata,MPIU_REPLACE);
1515:   return(0);
1516: }

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

1521:    Collective

1523:    Input Arguments:
1524: +  sf - star forest
1525: .  unit - data type
1526: -  multirootdata - root buffer to send to each leaf, one unit of data per leaf

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

1531:    Level: intermediate

1533: .seealso: PetscSFComputeDegreeBegin(), PetscSFScatterBegin()
1534: @*/
1535: PetscErrorCode PetscSFScatterBegin(PetscSF sf,MPI_Datatype unit,const void *multirootdata,void *leafdata)
1536: {
1538:   PetscSF        multi;

1542:   PetscSFSetUp(sf);
1543:   PetscSFGetMultiSF(sf,&multi);
1544:   PetscSFBcastBegin(multi,unit,multirootdata,leafdata);
1545:   return(0);
1546: }

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

1551:    Collective

1553:    Input Arguments:
1554: +  sf - star forest
1555: .  unit - data type
1556: -  multirootdata - root buffer to send to each leaf, one unit of data per leaf

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

1561:    Level: intermediate

1563: .seealso: PetscSFComputeDegreeEnd(), PetscSFScatterEnd()
1564: @*/
1565: PetscErrorCode PetscSFScatterEnd(PetscSF sf,MPI_Datatype unit,const void *multirootdata,void *leafdata)
1566: {
1568:   PetscSF        multi;

1572:   PetscSFSetUp(sf);
1573:   PetscSFGetMultiSF(sf,&multi);
1574:   PetscSFBcastEnd(multi,unit,multirootdata,leafdata);
1575:   return(0);
1576: }

1578: /*@
1579:   PetscSFCompose - Compose a new PetscSF equivalent to action to PetscSFs

1581:   Input Parameters:
1582: + sfA - The first PetscSF
1583: - sfB - The second PetscSF

1585:   Output Parameters:
1586: . sfBA - equvalent PetscSF for applying A then B

1588:   Level: developer

1590: .seealso: PetscSF, PetscSFGetGraph(), PetscSFSetGraph()
1591: @*/
1592: PetscErrorCode PetscSFCompose(PetscSF sfA, PetscSF sfB, PetscSF *sfBA)
1593: {
1594:   MPI_Comm           comm;
1595:   const PetscSFNode *remotePointsA, *remotePointsB;
1596:   PetscSFNode       *remotePointsBA;
1597:   const PetscInt    *localPointsA, *localPointsB;
1598:   PetscInt           numRootsA, numLeavesA, numRootsB, numLeavesB;
1599:   PetscErrorCode     ierr;

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