Actual source code: iscoloring.c

petsc-3.11.3 2019-06-26
Report Typos and Errors

  2:  #include <petsc/private/isimpl.h>
  3:  #include <petscviewer.h>
  4:  #include <petscsf.h>

  6: const char *const ISColoringTypes[] = {"global","ghosted","ISColoringType","IS_COLORING_",0};

  8: PetscErrorCode ISColoringReference(ISColoring coloring)
  9: {
 11:   (coloring)->refct++;
 12:   return(0);
 13: }

 15: PetscErrorCode ISColoringSetType(ISColoring coloring,ISColoringType type)
 16: {
 18:   (coloring)->ctype = type;
 19:   return(0);
 20: }

 22: /*@
 23:    ISColoringDestroy - Destroys a coloring context.

 25:    Collective on ISColoring

 27:    Input Parameter:
 28: .  iscoloring - the coloring context

 30:    Level: advanced

 32: .seealso: ISColoringView(), MatColoring
 33: @*/
 34: PetscErrorCode  ISColoringDestroy(ISColoring *iscoloring)
 35: {
 36:   PetscInt       i;

 40:   if (!*iscoloring) return(0);
 42:   if (--(*iscoloring)->refct > 0) {*iscoloring = 0; return(0);}

 44:   if ((*iscoloring)->is) {
 45:     for (i=0; i<(*iscoloring)->n; i++) {
 46:       ISDestroy(&(*iscoloring)->is[i]);
 47:     }
 48:     PetscFree((*iscoloring)->is);
 49:   }
 50:   if ((*iscoloring)->allocated) {PetscFree((*iscoloring)->colors);}
 51:   PetscCommDestroy(&(*iscoloring)->comm);
 52:   PetscFree((*iscoloring));
 53:   return(0);
 54: }

 56: /*
 57:   ISColoringViewFromOptions - Processes command line options to determine if/how an ISColoring object is to be viewed.

 59:   Collective on ISColoring

 61:   Input Parameters:
 62: + obj   - the ISColoring object
 63: . prefix - prefix to use for viewing, or NULL to use prefix of 'mat'
 64: - optionname - option to activate viewing

 66:   Level: intermediate

 68:   Developer Note: This cannot use PetscObjectViewFromOptions() because ISColoring is not a PetscObject

 70: */
 71: PetscErrorCode ISColoringViewFromOptions(ISColoring obj,PetscObject bobj,const char optionname[])
 72: {
 73:   PetscErrorCode    ierr;
 74:   PetscViewer       viewer;
 75:   PetscBool         flg;
 76:   PetscViewerFormat format;
 77:   char              *prefix;

 80:   prefix = bobj ? bobj->prefix : NULL;
 81:   PetscOptionsGetViewer(obj->comm,NULL,prefix,optionname,&viewer,&format,&flg);
 82:   if (flg) {
 83:     PetscViewerPushFormat(viewer,format);
 84:     ISColoringView(obj,viewer);
 85:     PetscViewerPopFormat(viewer);
 86:     PetscViewerDestroy(&viewer);
 87:   }
 88:   return(0);
 89: }

 91: /*@C
 92:    ISColoringView - Views a coloring context.

 94:    Collective on ISColoring

 96:    Input Parameters:
 97: +  iscoloring - the coloring context
 98: -  viewer - the viewer

100:    Level: advanced

102: .seealso: ISColoringDestroy(), ISColoringGetIS(), MatColoring
103: @*/
104: PetscErrorCode  ISColoringView(ISColoring iscoloring,PetscViewer viewer)
105: {
106:   PetscInt       i;
108:   PetscBool      iascii;
109:   IS             *is;

113:   if (!viewer) {
114:     PetscViewerASCIIGetStdout(iscoloring->comm,&viewer);
115:   }

118:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
119:   if (iascii) {
120:     MPI_Comm    comm;
121:     PetscMPIInt size,rank;

123:     PetscObjectGetComm((PetscObject)viewer,&comm);
124:     MPI_Comm_size(comm,&size);
125:     MPI_Comm_rank(comm,&rank);
126:     PetscViewerASCIIPrintf(viewer,"ISColoring Object: %d MPI processes\n",size);
127:     PetscViewerASCIIPushSynchronized(viewer);
128:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of colors %d\n",rank,iscoloring->n);
129:     PetscViewerFlush(viewer);
130:     PetscViewerASCIIPopSynchronized(viewer);
131:   }

133:   ISColoringGetIS(iscoloring,PETSC_IGNORE,&is);
134:   for (i=0; i<iscoloring->n; i++) {
135:     ISView(iscoloring->is[i],viewer);
136:   }
137:   ISColoringRestoreIS(iscoloring,&is);
138:   return(0);
139: }

141: /*@C
142:    ISColoringGetIS - Extracts index sets from the coloring context

144:    Collective on ISColoring

146:    Input Parameter:
147: .  iscoloring - the coloring context

149:    Output Parameters:
150: +  nn - number of index sets in the coloring context
151: -  is - array of index sets

153:    Level: advanced

155: .seealso: ISColoringRestoreIS(), ISColoringView()
156: @*/
157: PetscErrorCode  ISColoringGetIS(ISColoring iscoloring,PetscInt *nn,IS *isis[])
158: {


164:   if (nn) *nn = iscoloring->n;
165:   if (isis) {
166:     if (!iscoloring->is) {
167:       PetscInt        *mcolors,**ii,nc = iscoloring->n,i,base, n = iscoloring->N;
168:       ISColoringValue *colors = iscoloring->colors;
169:       IS              *is;

171: #if defined(PETSC_USE_DEBUG)
172:       for (i=0; i<n; i++) {
173:         if (((PetscInt)colors[i]) >= nc) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coloring is our of range index %d value %d number colors %d",(int)i,(int)colors[i],(int)nc);
174:       }
175: #endif

177:       /* generate the lists of nodes for each color */
178:       PetscCalloc1(nc,&mcolors);
179:       for (i=0; i<n; i++) mcolors[colors[i]]++;

181:       PetscMalloc1(nc,&ii);
182:       PetscMalloc1(n,&ii[0]);
183:       for (i=1; i<nc; i++) ii[i] = ii[i-1] + mcolors[i-1];
184:       PetscMemzero(mcolors,nc*sizeof(PetscInt));

186:       if (iscoloring->ctype == IS_COLORING_GLOBAL) {
187:         MPI_Scan(&iscoloring->N,&base,1,MPIU_INT,MPI_SUM,iscoloring->comm);
188:         base -= iscoloring->N;
189:         for (i=0; i<n; i++) ii[colors[i]][mcolors[colors[i]]++] = i + base; /* global idx */
190:       } else if (iscoloring->ctype == IS_COLORING_LOCAL) {
191:         for (i=0; i<n; i++) ii[colors[i]][mcolors[colors[i]]++] = i;   /* local idx */
192:       } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not provided for this ISColoringType type");

194:       PetscMalloc1(nc,&is);
195:       for (i=0; i<nc; i++) {
196:         ISCreateGeneral(iscoloring->comm,mcolors[i],ii[i],PETSC_COPY_VALUES,is+i);
197:       }

199:       iscoloring->is = is;
200:       PetscFree(ii[0]);
201:       PetscFree(ii);
202:       PetscFree(mcolors);
203:     }
204:     *isis = iscoloring->is;
205:   }
206:   return(0);
207: }

209: /*@C
210:    ISColoringRestoreIS - Restores the index sets extracted from the coloring context

212:    Collective on ISColoring

214:    Input Parameter:
215: +  iscoloring - the coloring context
216: -  is - array of index sets

218:    Level: advanced

220: .seealso: ISColoringGetIS(), ISColoringView()
221: @*/
222: PetscErrorCode  ISColoringRestoreIS(ISColoring iscoloring,IS *is[])
223: {

227:   /* currently nothing is done here */
228:   return(0);
229: }


232: /*@
233:     ISColoringCreate - Generates an ISColoring context from lists (provided
234:     by each processor) of colors for each node.

236:     Collective on MPI_Comm

238:     Input Parameters:
239: +   comm - communicator for the processors creating the coloring
240: .   ncolors - max color value
241: .   n - number of nodes on this processor
242: .   colors - array containing the colors for this processor, color numbers begin at 0.
243: -   mode - see PetscCopyMode for meaning of this flag.

245:     Output Parameter:
246: .   iscoloring - the resulting coloring data structure

248:     Options Database Key:
249: .   -is_coloring_view - Activates ISColoringView()

251:    Level: advanced

253:     Notes:
254:     By default sets coloring type to  IS_COLORING_GLOBAL

256: .seealso: MatColoringCreate(), ISColoringView(), ISColoringDestroy(), ISColoringSetType()

258: @*/
259: PetscErrorCode  ISColoringCreate(MPI_Comm comm,PetscInt ncolors,PetscInt n,const ISColoringValue colors[],PetscCopyMode mode,ISColoring *iscoloring)
260: {
262:   PetscMPIInt    size,rank,tag;
263:   PetscInt       base,top,i;
264:   PetscInt       nc,ncwork;
265:   MPI_Status     status;

268:   if (ncolors != PETSC_DECIDE && ncolors > IS_COLORING_MAX) {
269:     if (ncolors > 65535) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Max color value exeeds 65535 limit. This number is unrealistic. Perhaps a bug in code?\nCurrent max: %d user rewuested: %d",IS_COLORING_MAX,ncolors);
270:     else                 SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Max color value exeeds limit. Perhaps reconfigure PETSc with --with-is-color-value-type=short?\n Current max: %d user rewuested: %d",IS_COLORING_MAX,ncolors);
271:   }
272:   PetscNew(iscoloring);
273:   PetscCommDuplicate(comm,&(*iscoloring)->comm,&tag);
274:   comm = (*iscoloring)->comm;

276:   /* compute the number of the first node on my processor */
277:   MPI_Comm_size(comm,&size);

279:   /* should use MPI_Scan() */
280:   MPI_Comm_rank(comm,&rank);
281:   if (!rank) {
282:     base = 0;
283:     top  = n;
284:   } else {
285:     MPI_Recv(&base,1,MPIU_INT,rank-1,tag,comm,&status);
286:     top  = base+n;
287:   }
288:   if (rank < size-1) {
289:     MPI_Send(&top,1,MPIU_INT,rank+1,tag,comm);
290:   }

292:   /* compute the total number of colors */
293:   ncwork = 0;
294:   for (i=0; i<n; i++) {
295:     if (ncwork < colors[i]) ncwork = colors[i];
296:   }
297:   ncwork++;
298:   MPIU_Allreduce(&ncwork,&nc,1,MPIU_INT,MPI_MAX,comm);
299:   if (nc > ncolors) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of colors passed in %D is less then the actual number of colors in array %D",ncolors,nc);
300:   (*iscoloring)->n      = nc;
301:   (*iscoloring)->is     = 0;
302:   (*iscoloring)->N      = n;
303:   (*iscoloring)->refct  = 1;
304:   (*iscoloring)->ctype  = IS_COLORING_GLOBAL;
305:   if (mode == PETSC_COPY_VALUES) {
306:     PetscMalloc1(n,&(*iscoloring)->colors);
307:     PetscLogObjectMemory((PetscObject)(*iscoloring),n*sizeof(ISColoringValue));
308:     PetscMemcpy((*iscoloring)->colors,colors,n*sizeof(ISColoringValue));
309:     (*iscoloring)->allocated = PETSC_TRUE;
310:   } else if (mode == PETSC_OWN_POINTER) {
311:     (*iscoloring)->colors    = (ISColoringValue*)colors;
312:     (*iscoloring)->allocated = PETSC_TRUE;
313:   } else {
314:     (*iscoloring)->colors    = (ISColoringValue*)colors;
315:     (*iscoloring)->allocated = PETSC_FALSE;
316:   }
317:   ISColoringViewFromOptions(*iscoloring,NULL,"-is_coloring_view");
318:   PetscInfo1(0,"Number of colors %D\n",nc);
319:   return(0);
320: }

322: /*@
323:     ISBuildTwoSided - Takes an IS that describes where we will go. Generates an IS that contains new numbers from remote or local
324:     on the IS.

326:     Collective on IS

328:     Input Parameters
329: .   ito - an IS describes where we will go. Negative target rank will be ignored
330: .   toindx - an IS describes what indices should send. NULL means sending natural numbering

332:     Output Parameter:
333: .   rows - contains new numbers from remote or local

335:    Level: advanced

337: .seealso: MatPartitioningCreate(), ISPartitioningToNumbering(), ISPartitioningCount()

339: @*/
340: PetscErrorCode  ISBuildTwoSided(IS ito,IS toindx, IS *rows)
341: {
342:    const PetscInt *ito_indices,*toindx_indices;
343:    PetscInt       *send_indices,rstart,*recv_indices,nrecvs,nsends;
344:    PetscInt       *tosizes,*fromsizes,i,j,*tosizes_tmp,*tooffsets_tmp,ito_ln;
345:    PetscMPIInt    *toranks,*fromranks,size,target_rank,*fromperm_newtoold,nto,nfrom;
346:    PetscLayout     isrmap;
347:    MPI_Comm        comm;
348:    PetscSF         sf;
349:    PetscSFNode    *iremote;
350:    PetscErrorCode  ierr;

353:    PetscObjectGetComm((PetscObject)ito,&comm);
354:    MPI_Comm_size(comm,&size);
355:    ISGetLocalSize(ito,&ito_ln);
356:    /* why we do not have ISGetLayout? */
357:    isrmap = ito->map;
358:    PetscLayoutGetRange(isrmap,&rstart,NULL);
359:    ISGetIndices(ito,&ito_indices);
360:    PetscCalloc2(size,&tosizes_tmp,size+1,&tooffsets_tmp);
361:    for (i=0; i<ito_ln; i++) {
362:      if (ito_indices[i]<0) continue;
363: #if defined(PETSC_USE_DEBUG)
364:      if (ito_indices[i]>=size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"target rank %d is larger than communicator size %d ",ito_indices[i],size);
365: #endif
366:      tosizes_tmp[ito_indices[i]]++;
367:    }
368:    nto = 0;
369:    for (i=0; i<size; i++) {
370:      tooffsets_tmp[i+1] = tooffsets_tmp[i]+tosizes_tmp[i];
371:      if (tosizes_tmp[i]>0) nto++;
372:    }
373:    PetscCalloc2(nto,&toranks,2*nto,&tosizes);
374:    nto  = 0;
375:    for (i=0; i<size; i++) {
376:      if (tosizes_tmp[i]>0) {
377:        toranks[nto]     = i;
378:        tosizes[2*nto]   = tosizes_tmp[i];/* size */
379:        tosizes[2*nto+1] = tooffsets_tmp[i];/* offset */
380:        nto++;
381:      }
382:    }
383:    nsends = tooffsets_tmp[size];
384:    PetscCalloc1(nsends,&send_indices);
385:    if (toindx) {
386:      ISGetIndices(toindx,&toindx_indices);
387:    }
388:    for (i=0; i<ito_ln; i++) {
389:      if (ito_indices[i]<0) continue;
390:      target_rank = ito_indices[i];
391:      send_indices[tooffsets_tmp[target_rank]] = toindx? toindx_indices[i]:(i+rstart);
392:      tooffsets_tmp[target_rank]++;
393:    }
394:    if (toindx) {
395:      ISRestoreIndices(toindx,&toindx_indices);
396:    }
397:    ISRestoreIndices(ito,&ito_indices);
398:    PetscFree2(tosizes_tmp,tooffsets_tmp);
399:    PetscCommBuildTwoSided(comm,2,MPIU_INT,nto,toranks,tosizes,&nfrom,&fromranks,&fromsizes);
400:    PetscFree2(toranks,tosizes);
401:    PetscCalloc1(nfrom,&fromperm_newtoold);
402:    for (i=0; i<nfrom; i++) fromperm_newtoold[i] = i;
403:    PetscSortMPIIntWithArray(nfrom,fromranks,fromperm_newtoold);
404:    nrecvs = 0;
405:    for (i=0; i<nfrom; i++) nrecvs += fromsizes[i*2];
406:    PetscCalloc1(nrecvs,&recv_indices);
407:    PetscCalloc1(nrecvs,&iremote);
408:    nrecvs = 0;
409:    for (i=0; i<nfrom; i++) {
410:      for (j=0; j<fromsizes[2*fromperm_newtoold[i]]; j++) {
411:        iremote[nrecvs].rank    = fromranks[i];
412:        iremote[nrecvs++].index = fromsizes[2*fromperm_newtoold[i]+1]+j;
413:      }
414:    }
415:    PetscSFCreate(comm,&sf);
416:    PetscSFSetGraph(sf,nsends,nrecvs,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
417:    PetscSFSetType(sf,PETSCSFBASIC);
418:    /* how to put a prefix ? */
419:    PetscSFSetFromOptions(sf);
420:    PetscSFBcastBegin(sf,MPIU_INT,send_indices,recv_indices);
421:    PetscSFBcastEnd(sf,MPIU_INT,send_indices,recv_indices);
422:    PetscSFDestroy(&sf);
423:    PetscFree(fromranks);
424:    PetscFree(fromsizes);
425:    PetscFree(fromperm_newtoold);
426:    PetscFree(send_indices);
427:    if (rows) {
428:      PetscSortInt(nrecvs,recv_indices);
429:      ISCreateGeneral(comm,nrecvs,recv_indices,PETSC_OWN_POINTER,rows);
430:    } else {
431:      PetscFree(recv_indices);
432:    }
433:    return(0);
434: }


437: /*@
438:     ISPartitioningToNumbering - Takes an ISPartitioning and on each processor
439:     generates an IS that contains a new global node number for each index based
440:     on the partitioing.

442:     Collective on IS

444:     Input Parameters
445: .   partitioning - a partitioning as generated by MatPartitioningApply()
446:                    or MatPartitioningApplyND()

448:     Output Parameter:
449: .   is - on each processor the index set that defines the global numbers
450:          (in the new numbering) for all the nodes currently (before the partitioning)
451:          on that processor

453:    Level: advanced

455: .seealso: MatPartitioningCreate(), AOCreateBasic(), ISPartitioningCount()

457: @*/
458: PetscErrorCode  ISPartitioningToNumbering(IS part,IS *is)
459: {
460:   MPI_Comm       comm;
461:   IS             ndorder;
462:   PetscInt       i,np,npt,n,*starts = NULL,*sums = NULL,*lsizes = NULL,*newi = NULL;
463:   const PetscInt *indices = NULL;

469:   /* see if the partitioning comes from nested dissection */
470:   PetscObjectQuery((PetscObject)part,"_petsc_matpartitioning_ndorder",(PetscObject*)&ndorder);
471:   if (ndorder) {
472:     PetscObjectReference((PetscObject)ndorder);
473:     *is  = ndorder;
474:     return(0);
475:   }

477:   PetscObjectGetComm((PetscObject)part,&comm);
478:   /* count the number of partitions, i.e., virtual processors */
479:   ISGetLocalSize(part,&n);
480:   ISGetIndices(part,&indices);
481:   np   = 0;
482:   for (i=0; i<n; i++) np = PetscMax(np,indices[i]);
483:   MPIU_Allreduce(&np,&npt,1,MPIU_INT,MPI_MAX,comm);
484:   np   = npt+1; /* so that it looks like a MPI_Comm_size output */

486:   /*
487:         lsizes - number of elements of each partition on this particular processor
488:         sums - total number of "previous" nodes for any particular partition
489:         starts - global number of first element in each partition on this processor
490:   */
491:   PetscMalloc3(np,&lsizes,np,&starts,np,&sums);
492:   PetscMemzero(lsizes,np*sizeof(PetscInt));
493:   for (i=0; i<n; i++) lsizes[indices[i]]++;
494:   MPIU_Allreduce(lsizes,sums,np,MPIU_INT,MPI_SUM,comm);
495:   MPI_Scan(lsizes,starts,np,MPIU_INT,MPI_SUM,comm);
496:   for (i=0; i<np; i++) starts[i] -= lsizes[i];
497:   for (i=1; i<np; i++) {
498:     sums[i]   += sums[i-1];
499:     starts[i] += sums[i-1];
500:   }

502:   /*
503:       For each local index give it the new global number
504:   */
505:   PetscMalloc1(n,&newi);
506:   for (i=0; i<n; i++) newi[i] = starts[indices[i]]++;
507:   PetscFree3(lsizes,starts,sums);

509:   ISRestoreIndices(part,&indices);
510:   ISCreateGeneral(comm,n,newi,PETSC_OWN_POINTER,is);
511:   ISSetPermutation(*is);
512:   return(0);
513: }

515: /*@
516:     ISPartitioningCount - Takes a ISPartitioning and determines the number of
517:     resulting elements on each (partition) process

519:     Collective on IS

521:     Input Parameters:
522: +   partitioning - a partitioning as generated by MatPartitioningApply() or
523:                    MatPartitioningApplyND()
524: -   len - length of the array count, this is the total number of partitions

526:     Output Parameter:
527: .   count - array of length size, to contain the number of elements assigned
528:         to each partition, where size is the number of partitions generated
529:          (see notes below).

531:    Level: advanced

533:     Notes:
534:         By default the number of partitions generated (and thus the length
535:         of count) is the size of the communicator associated with IS,
536:         but it can be set by MatPartitioningSetNParts. The resulting array
537:         of lengths can for instance serve as input of PCBJacobiSetTotalBlocks.
538:         If the partitioning has been obtained by MatPartitioningApplyND(),
539:         the returned count does not include the separators.

541: .seealso: MatPartitioningCreate(), AOCreateBasic(), ISPartitioningToNumbering(),
542:         MatPartitioningSetNParts(), MatPartitioningApply(), MatPartitioningApplyND()

544: @*/
545: PetscErrorCode  ISPartitioningCount(IS part,PetscInt len,PetscInt count[])
546: {
547:   MPI_Comm       comm;
548:   PetscInt       i,n,*lsizes;
549:   const PetscInt *indices;
551:   PetscMPIInt    npp;

554:   PetscObjectGetComm((PetscObject)part,&comm);
555:   if (len == PETSC_DEFAULT) {
556:     PetscMPIInt size;
557:     MPI_Comm_size(comm,&size);
558:     len  = (PetscInt) size;
559:   }

561:   /* count the number of partitions */
562:   ISGetLocalSize(part,&n);
563:   ISGetIndices(part,&indices);
564: #if defined(PETSC_USE_DEBUG)
565:   {
566:     PetscInt np = 0,npt;
567:     for (i=0; i<n; i++) np = PetscMax(np,indices[i]);
568:     MPIU_Allreduce(&np,&npt,1,MPIU_INT,MPI_MAX,comm);
569:     np   = npt+1; /* so that it looks like a MPI_Comm_size output */
570:     if (np > len) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Length of count array %D is less than number of partitions %D",len,np);
571:   }
572: #endif

574:   /*
575:         lsizes - number of elements of each partition on this particular processor
576:         sums - total number of "previous" nodes for any particular partition
577:         starts - global number of first element in each partition on this processor
578:   */
579:   PetscCalloc1(len,&lsizes);
580:   for (i=0; i<n; i++) {
581:     if (indices[i] > -1) lsizes[indices[i]]++;
582:   }
583:   ISRestoreIndices(part,&indices);
584:   PetscMPIIntCast(len,&npp);
585:   MPIU_Allreduce(lsizes,count,npp,MPIU_INT,MPI_SUM,comm);
586:   PetscFree(lsizes);
587:   return(0);
588: }

590: /*@
591:     ISAllGather - Given an index set (IS) on each processor, generates a large
592:     index set (same on each processor) by concatenating together each
593:     processors index set.

595:     Collective on IS

597:     Input Parameter:
598: .   is - the distributed index set

600:     Output Parameter:
601: .   isout - the concatenated index set (same on all processors)

603:     Notes:
604:     ISAllGather() is clearly not scalable for large index sets.

606:     The IS created on each processor must be created with a common
607:     communicator (e.g., PETSC_COMM_WORLD). If the index sets were created
608:     with PETSC_COMM_SELF, this routine will not work as expected, since
609:     each process will generate its own new IS that consists only of
610:     itself.

612:     The communicator for this new IS is PETSC_COMM_SELF

614:     Level: intermediate

616:     Concepts: gather^index sets
617:     Concepts: index sets^gathering to all processors
618:     Concepts: IS^gathering to all processors

620: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock()
621: @*/
622: PetscErrorCode  ISAllGather(IS is,IS *isout)
623: {
625:   PetscInt       *indices,n,i,N,step,first;
626:   const PetscInt *lindices;
627:   MPI_Comm       comm;
628:   PetscMPIInt    size,*sizes = NULL,*offsets = NULL,nn;
629:   PetscBool      stride;


635:   PetscObjectGetComm((PetscObject)is,&comm);
636:   MPI_Comm_size(comm,&size);
637:   ISGetLocalSize(is,&n);
638:   PetscObjectTypeCompare((PetscObject)is,ISSTRIDE,&stride);
639:   if (size == 1 && stride) { /* should handle parallel ISStride also */
640:     ISStrideGetInfo(is,&first,&step);
641:     ISCreateStride(PETSC_COMM_SELF,n,first,step,isout);
642:   } else {
643:     PetscMalloc2(size,&sizes,size,&offsets);

645:     PetscMPIIntCast(n,&nn);
646:     MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
647:     offsets[0] = 0;
648:     for (i=1; i<size; i++) {
649:       PetscInt s = offsets[i-1] + sizes[i-1];
650:       PetscMPIIntCast(s,&offsets[i]);
651:     }
652:     N = offsets[size-1] + sizes[size-1];

654:     PetscMalloc1(N,&indices);
655:     ISGetIndices(is,&lindices);
656:     MPI_Allgatherv((void*)lindices,nn,MPIU_INT,indices,sizes,offsets,MPIU_INT,comm);
657:     ISRestoreIndices(is,&lindices);
658:     PetscFree2(sizes,offsets);

660:     ISCreateGeneral(PETSC_COMM_SELF,N,indices,PETSC_OWN_POINTER,isout);
661:   }
662:   return(0);
663: }

665: /*@C
666:     ISAllGatherColors - Given a a set of colors on each processor, generates a large
667:     set (same on each processor) by concatenating together each processors colors

669:     Collective on MPI_Comm

671:     Input Parameter:
672: +   comm - communicator to share the indices
673: .   n - local size of set
674: -   lindices - local colors

676:     Output Parameter:
677: +   outN - total number of indices
678: -   outindices - all of the colors

680:     Notes:
681:     ISAllGatherColors() is clearly not scalable for large index sets.


684:     Level: intermediate

686:     Concepts: gather^index sets
687:     Concepts: index sets^gathering to all processors
688:     Concepts: IS^gathering to all processors

690: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
691: @*/
692: PetscErrorCode  ISAllGatherColors(MPI_Comm comm,PetscInt n,ISColoringValue *lindices,PetscInt *outN,ISColoringValue *outindices[])
693: {
694:   ISColoringValue *indices;
695:   PetscErrorCode  ierr;
696:   PetscInt        i,N;
697:   PetscMPIInt     size,*offsets = NULL,*sizes = NULL, nn = n;

700:   MPI_Comm_size(comm,&size);
701:   PetscMalloc2(size,&sizes,size,&offsets);

703:   MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
704:   offsets[0] = 0;
705:   for (i=1; i<size; i++) offsets[i] = offsets[i-1] + sizes[i-1];
706:   N    = offsets[size-1] + sizes[size-1];
707:   PetscFree2(sizes,offsets);

709:   PetscMalloc1(N+1,&indices);
710:   MPI_Allgatherv(lindices,(PetscMPIInt)n,MPIU_COLORING_VALUE,indices,sizes,offsets,MPIU_COLORING_VALUE,comm);

712:   *outindices = indices;
713:   if (outN) *outN = N;
714:   return(0);
715: }

717: /*@
718:     ISComplement - Given an index set (IS) generates the complement index set. That is all
719:        all indices that are NOT in the given set.

721:     Collective on IS

723:     Input Parameter:
724: +   is - the index set
725: .   nmin - the first index desired in the local part of the complement
726: -   nmax - the largest index desired in the local part of the complement (note that all indices in is must be greater or equal to nmin and less than nmax)

728:     Output Parameter:
729: .   isout - the complement

731:     Notes:
732:     The communicator for this new IS is the same as for the input IS

734:       For a parallel IS, this will generate the local part of the complement on each process

736:       To generate the entire complement (on each process) of a parallel IS, first call ISAllGather() and then
737:     call this routine.

739:     Level: intermediate

741:     Concepts: gather^index sets
742:     Concepts: index sets^gathering to all processors
743:     Concepts: IS^gathering to all processors

745: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
746: @*/
747: PetscErrorCode  ISComplement(IS is,PetscInt nmin,PetscInt nmax,IS *isout)
748: {
750:   const PetscInt *indices;
751:   PetscInt       n,i,j,unique,cnt,*nindices;
752:   PetscBool      sorted;

757:   if (nmin < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nmin %D cannot be negative",nmin);
758:   if (nmin > nmax) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nmin %D cannot be greater than nmax %D",nmin,nmax);
759:   ISSorted(is,&sorted);
760:   if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Index set must be sorted");

762:   ISGetLocalSize(is,&n);
763:   ISGetIndices(is,&indices);
764: #if defined(PETSC_USE_DEBUG)
765:   for (i=0; i<n; i++) {
766:     if (indices[i] <  nmin) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Index %D's value %D is smaller than minimum given %D",i,indices[i],nmin);
767:     if (indices[i] >= nmax) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Index %D's value %D is larger than maximum given %D",i,indices[i],nmax);
768:   }
769: #endif
770:   /* Count number of unique entries */
771:   unique = (n>0);
772:   for (i=0; i<n-1; i++) {
773:     if (indices[i+1] != indices[i]) unique++;
774:   }
775:   PetscMalloc1(nmax-nmin-unique,&nindices);
776:   cnt  = 0;
777:   for (i=nmin,j=0; i<nmax; i++) {
778:     if (j<n && i==indices[j]) do { j++; } while (j<n && i==indices[j]);
779:     else nindices[cnt++] = i;
780:   }
781:   if (cnt != nmax-nmin-unique) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Number of entries found in complement %D does not match expected %D",cnt,nmax-nmin-unique);
782:   ISCreateGeneral(PetscObjectComm((PetscObject)is),cnt,nindices,PETSC_OWN_POINTER,isout);
783:   ISRestoreIndices(is,&indices);
784:   return(0);
785: }