Actual source code: iscoloring.c


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

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

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

 14: /*@C

 16:     ISColoringSetType - indicates if the coloring is for the local representation (including ghost points) or the global representation

 18:    Collective on coloring

 20:    Input Parameters:
 21: +    coloring - the coloring object
 22: -    type - either IS_COLORING_LOCAL or IS_COLORING_GLOBAL

 24:    Notes:
 25:      With IS_COLORING_LOCAL the coloring is in the numbering of the local vector, for IS_COLORING_GLOBAL it is in the number of the global vector

 27:    Level: intermediate

 29: .seealso: MatFDColoringCreate(), ISColoring, ISColoringCreate(), IS_COLORING_LOCAL, IS_COLORING_GLOBAL, ISColoringGetType()

 31: @*/
 32: PetscErrorCode ISColoringSetType(ISColoring coloring,ISColoringType type)
 33: {
 34:   coloring->ctype = type;
 35:   return 0;
 36: }

 38: /*@C

 40:     ISColoringGetType - gets if the coloring is for the local representation (including ghost points) or the global representation

 42:    Collective on coloring

 44:    Input Parameter:
 45: .   coloring - the coloring object

 47:    Output Parameter:
 48: .    type - either IS_COLORING_LOCAL or IS_COLORING_GLOBAL

 50:    Level: intermediate

 52: .seealso: MatFDColoringCreate(), ISColoring, ISColoringCreate(), IS_COLORING_LOCAL, IS_COLORING_GLOBAL, ISColoringSetType()

 54: @*/
 55: PetscErrorCode ISColoringGetType(ISColoring coloring,ISColoringType *type)
 56: {
 57:   *type = coloring->ctype;
 58:   return 0;
 59: }

 61: /*@
 62:    ISColoringDestroy - Destroys a coloring context.

 64:    Collective on ISColoring

 66:    Input Parameter:
 67: .  iscoloring - the coloring context

 69:    Level: advanced

 71: .seealso: ISColoringView(), MatColoring
 72: @*/
 73: PetscErrorCode  ISColoringDestroy(ISColoring *iscoloring)
 74: {
 75:   PetscInt       i;

 77:   if (!*iscoloring) return 0;
 79:   if (--(*iscoloring)->refct > 0) {*iscoloring = NULL; return 0;}

 81:   if ((*iscoloring)->is) {
 82:     for (i=0; i<(*iscoloring)->n; i++) {
 83:       ISDestroy(&(*iscoloring)->is[i]);
 84:     }
 85:     PetscFree((*iscoloring)->is);
 86:   }
 87:   if ((*iscoloring)->allocated) PetscFree((*iscoloring)->colors);
 88:   PetscCommDestroy(&(*iscoloring)->comm);
 89:   PetscFree((*iscoloring));
 90:   return 0;
 91: }

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

 96:   Collective on ISColoring

 98:   Input Parameters:
 99: + obj   - the ISColoring object
100: . prefix - prefix to use for viewing, or NULL to use prefix of 'mat'
101: - optionname - option to activate viewing

103:   Level: intermediate

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

107: */
108: PetscErrorCode ISColoringViewFromOptions(ISColoring obj,PetscObject bobj,const char optionname[])
109: {
110:   PetscViewer       viewer;
111:   PetscBool         flg;
112:   PetscViewerFormat format;
113:   char              *prefix;

115:   prefix = bobj ? bobj->prefix : NULL;
116:   PetscOptionsGetViewer(obj->comm,NULL,prefix,optionname,&viewer,&format,&flg);
117:   if (flg) {
118:     PetscViewerPushFormat(viewer,format);
119:     ISColoringView(obj,viewer);
120:     PetscViewerPopFormat(viewer);
121:     PetscViewerDestroy(&viewer);
122:   }
123:   return 0;
124: }

126: /*@C
127:    ISColoringView - Views a coloring context.

129:    Collective on ISColoring

131:    Input Parameters:
132: +  iscoloring - the coloring context
133: -  viewer - the viewer

135:    Level: advanced

137: .seealso: ISColoringDestroy(), ISColoringGetIS(), MatColoring
138: @*/
139: PetscErrorCode  ISColoringView(ISColoring iscoloring,PetscViewer viewer)
140: {
141:   PetscInt       i;
142:   PetscBool      iascii;
143:   IS             *is;

146:   if (!viewer) {
147:     PetscViewerASCIIGetStdout(iscoloring->comm,&viewer);
148:   }

151:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
152:   if (iascii) {
153:     MPI_Comm    comm;
154:     PetscMPIInt size,rank;

156:     PetscObjectGetComm((PetscObject)viewer,&comm);
157:     MPI_Comm_size(comm,&size);
158:     MPI_Comm_rank(comm,&rank);
159:     PetscViewerASCIIPrintf(viewer,"ISColoring Object: %d MPI processes\n",size);
160:     PetscViewerASCIIPrintf(viewer,"ISColoringType: %s\n",ISColoringTypes[iscoloring->ctype]);
161:     PetscViewerASCIIPushSynchronized(viewer);
162:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of colors %" PetscInt_FMT "\n",rank,iscoloring->n);
163:     PetscViewerFlush(viewer);
164:     PetscViewerASCIIPopSynchronized(viewer);
165:   }

167:   ISColoringGetIS(iscoloring,PETSC_USE_POINTER,PETSC_IGNORE,&is);
168:   for (i=0; i<iscoloring->n; i++) {
169:     ISView(iscoloring->is[i],viewer);
170:   }
171:   ISColoringRestoreIS(iscoloring,PETSC_USE_POINTER,&is);
172:   return 0;
173: }

175: /*@C
176:    ISColoringGetColors - Returns an array with the color for each node

178:    Not Collective

180:    Input Parameter:
181: .  iscoloring - the coloring context

183:    Output Parameters:
184: +  n - number of nodes
185: .  nc - number of colors
186: -  colors - color for each node

188:    Level: advanced

190: .seealso: ISColoringRestoreIS(), ISColoringView(), ISColoringGetIS()
191: @*/
192: PetscErrorCode  ISColoringGetColors(ISColoring iscoloring,PetscInt *n,PetscInt *nc,const ISColoringValue **colors)
193: {

196:   if (n) *n = iscoloring->N;
197:   if (nc) *nc = iscoloring->n;
198:   if (colors) *colors = iscoloring->colors;
199:   return 0;
200: }

202: /*@C
203:    ISColoringGetIS - Extracts index sets from the coloring context. Each is contains the nodes of one color

205:    Collective on ISColoring

207:    Input Parameters:
208: +  iscoloring - the coloring context
209: -  mode - if this value is PETSC_OWN_POINTER then the caller owns the pointer and must free the array of IS and each IS in the array

211:    Output Parameters:
212: +  nn - number of index sets in the coloring context
213: -  is - array of index sets

215:    Level: advanced

217: .seealso: ISColoringRestoreIS(), ISColoringView(), ISColoringGetColoring()
218: @*/
219: PetscErrorCode  ISColoringGetIS(ISColoring iscoloring,PetscCopyMode mode, PetscInt *nn,IS *isis[])
220: {

223:   if (nn) *nn = iscoloring->n;
224:   if (isis) {
225:     if (!iscoloring->is) {
226:       PetscInt        *mcolors,**ii,nc = iscoloring->n,i,base, n = iscoloring->N;
227:       ISColoringValue *colors = iscoloring->colors;
228:       IS              *is;

230:       if (PetscDefined(USE_DEBUG)) {
231:         for (i=0; i<n; i++) {
233:         }
234:       }

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

240:       PetscMalloc1(nc,&ii);
241:       PetscMalloc1(n,&ii[0]);
242:       for (i=1; i<nc; i++) ii[i] = ii[i-1] + mcolors[i-1];
243:       PetscArrayzero(mcolors,nc);

245:       if (iscoloring->ctype == IS_COLORING_GLOBAL) {
246:         MPI_Scan(&iscoloring->N,&base,1,MPIU_INT,MPI_SUM,iscoloring->comm);
247:         base -= iscoloring->N;
248:         for (i=0; i<n; i++) ii[colors[i]][mcolors[colors[i]]++] = i + base; /* global idx */
249:       } else if (iscoloring->ctype == IS_COLORING_LOCAL) {
250:         for (i=0; i<n; i++) ii[colors[i]][mcolors[colors[i]]++] = i;   /* local idx */
251:       } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not provided for this ISColoringType type");

253:       PetscMalloc1(nc,&is);
254:       for (i=0; i<nc; i++) {
255:         ISCreateGeneral(iscoloring->comm,mcolors[i],ii[i],PETSC_COPY_VALUES,is+i);
256:       }

258:       if (mode != PETSC_OWN_POINTER) iscoloring->is = is;
259:       *isis = is;
260:       PetscFree(ii[0]);
261:       PetscFree(ii);
262:       PetscFree(mcolors);
263:     } else {
264:       *isis = iscoloring->is;
265:     }
266:   }
267:   return 0;
268: }

270: /*@C
271:    ISColoringRestoreIS - Restores the index sets extracted from the coloring context

273:    Collective on ISColoring

275:    Input Parameters:
276: +  iscoloring - the coloring context
277: .  mode - who retains ownership of the is
278: -  is - array of index sets

280:    Level: advanced

282: .seealso: ISColoringGetIS(), ISColoringView()
283: @*/
284: PetscErrorCode  ISColoringRestoreIS(ISColoring iscoloring,PetscCopyMode mode,IS *is[])
285: {

288:   /* currently nothing is done here */
289:   return 0;
290: }

292: /*@
293:     ISColoringCreate - Generates an ISColoring context from lists (provided
294:     by each processor) of colors for each node.

296:     Collective

298:     Input Parameters:
299: +   comm - communicator for the processors creating the coloring
300: .   ncolors - max color value
301: .   n - number of nodes on this processor
302: .   colors - array containing the colors for this processor, color numbers begin at 0.
303: -   mode - see PetscCopyMode for meaning of this flag.

305:     Output Parameter:
306: .   iscoloring - the resulting coloring data structure

308:     Options Database Key:
309: .   -is_coloring_view - Activates ISColoringView()

311:    Level: advanced

313:     Notes:
314:     By default sets coloring type to  IS_COLORING_GLOBAL

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

318: @*/
319: PetscErrorCode  ISColoringCreate(MPI_Comm comm,PetscInt ncolors,PetscInt n,const ISColoringValue colors[],PetscCopyMode mode,ISColoring *iscoloring)
320: {
321:   PetscMPIInt    size,rank,tag;
322:   PetscInt       base,top,i;
323:   PetscInt       nc,ncwork;
324:   MPI_Status     status;

326:   if (ncolors != PETSC_DECIDE && ncolors > IS_COLORING_MAX) {
328:     else                 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Max color value exceeds limit. Perhaps reconfigure PETSc with --with-is-color-value-type=short?\n Current max: %d user requested: %" PetscInt_FMT,PETSC_IS_COLORING_MAX,ncolors);
329:   }
330:   PetscNew(iscoloring);
331:   PetscCommDuplicate(comm,&(*iscoloring)->comm,&tag);
332:   comm = (*iscoloring)->comm;

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

337:   /* should use MPI_Scan() */
338:   MPI_Comm_rank(comm,&rank);
339:   if (rank == 0) {
340:     base = 0;
341:     top  = n;
342:   } else {
343:     MPI_Recv(&base,1,MPIU_INT,rank-1,tag,comm,&status);
344:     top  = base+n;
345:   }
346:   if (rank < size-1) {
347:     MPI_Send(&top,1,MPIU_INT,rank+1,tag,comm);
348:   }

350:   /* compute the total number of colors */
351:   ncwork = 0;
352:   for (i=0; i<n; i++) {
353:     if (ncwork < colors[i]) ncwork = colors[i];
354:   }
355:   ncwork++;
356:   MPIU_Allreduce(&ncwork,&nc,1,MPIU_INT,MPI_MAX,comm);
358:   (*iscoloring)->n      = nc;
359:   (*iscoloring)->is     = NULL;
360:   (*iscoloring)->N      = n;
361:   (*iscoloring)->refct  = 1;
362:   (*iscoloring)->ctype  = IS_COLORING_GLOBAL;
363:   if (mode == PETSC_COPY_VALUES) {
364:     PetscMalloc1(n,&(*iscoloring)->colors);
365:     PetscLogObjectMemory((PetscObject)(*iscoloring),n*sizeof(ISColoringValue));
366:     PetscArraycpy((*iscoloring)->colors,colors,n);
367:     (*iscoloring)->allocated = PETSC_TRUE;
368:   } else if (mode == PETSC_OWN_POINTER) {
369:     (*iscoloring)->colors    = (ISColoringValue*)colors;
370:     (*iscoloring)->allocated = PETSC_TRUE;
371:   } else {
372:     (*iscoloring)->colors    = (ISColoringValue*)colors;
373:     (*iscoloring)->allocated = PETSC_FALSE;
374:   }
375:   ISColoringViewFromOptions(*iscoloring,NULL,"-is_coloring_view");
376:   PetscInfo(0,"Number of colors %" PetscInt_FMT "\n",nc);
377:   return 0;
378: }

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

384:     Collective on IS

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

390:     Output Parameter:
391: .   rows - contains new numbers from remote or local

393:    Level: advanced

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

397: @*/
398: PetscErrorCode  ISBuildTwoSided(IS ito,IS toindx, IS *rows)
399: {
400:    const PetscInt *ito_indices,*toindx_indices;
401:    PetscInt       *send_indices,rstart,*recv_indices,nrecvs,nsends;
402:    PetscInt       *tosizes,*fromsizes,i,j,*tosizes_tmp,*tooffsets_tmp,ito_ln;
403:    PetscMPIInt    *toranks,*fromranks,size,target_rank,*fromperm_newtoold,nto,nfrom;
404:    PetscLayout     isrmap;
405:    MPI_Comm        comm;
406:    PetscSF         sf;
407:    PetscSFNode    *iremote;

409:    PetscObjectGetComm((PetscObject)ito,&comm);
410:    MPI_Comm_size(comm,&size);
411:    ISGetLocalSize(ito,&ito_ln);
412:    ISGetLayout(ito,&isrmap);
413:    PetscLayoutGetRange(isrmap,&rstart,NULL);
414:    ISGetIndices(ito,&ito_indices);
415:    PetscCalloc2(size,&tosizes_tmp,size+1,&tooffsets_tmp);
416:    for (i=0; i<ito_ln; i++) {
417:      if (ito_indices[i]<0) continue;
419:      tosizes_tmp[ito_indices[i]]++;
420:    }
421:    nto = 0;
422:    for (i=0; i<size; i++) {
423:      tooffsets_tmp[i+1] = tooffsets_tmp[i]+tosizes_tmp[i];
424:      if (tosizes_tmp[i]>0) nto++;
425:    }
426:    PetscCalloc2(nto,&toranks,2*nto,&tosizes);
427:    nto  = 0;
428:    for (i=0; i<size; i++) {
429:      if (tosizes_tmp[i]>0) {
430:        toranks[nto]     = i;
431:        tosizes[2*nto]   = tosizes_tmp[i];/* size */
432:        tosizes[2*nto+1] = tooffsets_tmp[i];/* offset */
433:        nto++;
434:      }
435:    }
436:    nsends = tooffsets_tmp[size];
437:    PetscCalloc1(nsends,&send_indices);
438:    if (toindx) {
439:      ISGetIndices(toindx,&toindx_indices);
440:    }
441:    for (i=0; i<ito_ln; i++) {
442:      if (ito_indices[i]<0) continue;
443:      target_rank = ito_indices[i];
444:      send_indices[tooffsets_tmp[target_rank]] = toindx? toindx_indices[i]:(i+rstart);
445:      tooffsets_tmp[target_rank]++;
446:    }
447:    if (toindx) {
448:      ISRestoreIndices(toindx,&toindx_indices);
449:    }
450:    ISRestoreIndices(ito,&ito_indices);
451:    PetscFree2(tosizes_tmp,tooffsets_tmp);
452:    PetscCommBuildTwoSided(comm,2,MPIU_INT,nto,toranks,tosizes,&nfrom,&fromranks,&fromsizes);
453:    PetscFree2(toranks,tosizes);
454:    PetscMalloc1(nfrom,&fromperm_newtoold);
455:    for (i=0; i<nfrom; i++) fromperm_newtoold[i] = i;
456:    PetscSortMPIIntWithArray(nfrom,fromranks,fromperm_newtoold);
457:    nrecvs = 0;
458:    for (i=0; i<nfrom; i++) nrecvs += fromsizes[i*2];
459:    PetscCalloc1(nrecvs,&recv_indices);
460:    PetscMalloc1(nrecvs,&iremote);
461:    nrecvs = 0;
462:    for (i=0; i<nfrom; i++) {
463:      for (j=0; j<fromsizes[2*fromperm_newtoold[i]]; j++) {
464:        iremote[nrecvs].rank    = fromranks[i];
465:        iremote[nrecvs++].index = fromsizes[2*fromperm_newtoold[i]+1]+j;
466:      }
467:    }
468:    PetscSFCreate(comm,&sf);
469:    PetscSFSetGraph(sf,nsends,nrecvs,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
470:    PetscSFSetType(sf,PETSCSFBASIC);
471:    /* how to put a prefix ? */
472:    PetscSFSetFromOptions(sf);
473:    PetscSFBcastBegin(sf,MPIU_INT,send_indices,recv_indices,MPI_REPLACE);
474:    PetscSFBcastEnd(sf,MPIU_INT,send_indices,recv_indices,MPI_REPLACE);
475:    PetscSFDestroy(&sf);
476:    PetscFree(fromranks);
477:    PetscFree(fromsizes);
478:    PetscFree(fromperm_newtoold);
479:    PetscFree(send_indices);
480:    if (rows) {
481:      PetscSortInt(nrecvs,recv_indices);
482:      ISCreateGeneral(comm,nrecvs,recv_indices,PETSC_OWN_POINTER,rows);
483:    } else {
484:      PetscFree(recv_indices);
485:    }
486:    return 0;
487: }

489: /*@
490:     ISPartitioningToNumbering - Takes an ISPartitioning and on each processor
491:     generates an IS that contains a new global node number for each index based
492:     on the partitioing.

494:     Collective on IS

496:     Input Parameters:
497: .   partitioning - a partitioning as generated by MatPartitioningApply()
498:                    or MatPartitioningApplyND()

500:     Output Parameter:
501: .   is - on each processor the index set that defines the global numbers
502:          (in the new numbering) for all the nodes currently (before the partitioning)
503:          on that processor

505:    Level: advanced

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

509: @*/
510: PetscErrorCode  ISPartitioningToNumbering(IS part,IS *is)
511: {
512:   MPI_Comm       comm;
513:   IS             ndorder;
514:   PetscInt       i,np,npt,n,*starts = NULL,*sums = NULL,*lsizes = NULL,*newi = NULL;
515:   const PetscInt *indices = NULL;

519:   /* see if the partitioning comes from nested dissection */
520:   PetscObjectQuery((PetscObject)part,"_petsc_matpartitioning_ndorder",(PetscObject*)&ndorder);
521:   if (ndorder) {
522:     PetscObjectReference((PetscObject)ndorder);
523:     *is  = ndorder;
524:     return 0;
525:   }

527:   PetscObjectGetComm((PetscObject)part,&comm);
528:   /* count the number of partitions, i.e., virtual processors */
529:   ISGetLocalSize(part,&n);
530:   ISGetIndices(part,&indices);
531:   np   = 0;
532:   for (i=0; i<n; i++) np = PetscMax(np,indices[i]);
533:   MPIU_Allreduce(&np,&npt,1,MPIU_INT,MPI_MAX,comm);
534:   np   = npt+1; /* so that it looks like a MPI_Comm_size output */

536:   /*
537:         lsizes - number of elements of each partition on this particular processor
538:         sums - total number of "previous" nodes for any particular partition
539:         starts - global number of first element in each partition on this processor
540:   */
541:   PetscMalloc3(np,&lsizes,np,&starts,np,&sums);
542:   PetscArrayzero(lsizes,np);
543:   for (i=0; i<n; i++) lsizes[indices[i]]++;
544:   MPIU_Allreduce(lsizes,sums,np,MPIU_INT,MPI_SUM,comm);
545:   MPI_Scan(lsizes,starts,np,MPIU_INT,MPI_SUM,comm);
546:   for (i=0; i<np; i++) starts[i] -= lsizes[i];
547:   for (i=1; i<np; i++) {
548:     sums[i]   += sums[i-1];
549:     starts[i] += sums[i-1];
550:   }

552:   /*
553:       For each local index give it the new global number
554:   */
555:   PetscMalloc1(n,&newi);
556:   for (i=0; i<n; i++) newi[i] = starts[indices[i]]++;
557:   PetscFree3(lsizes,starts,sums);

559:   ISRestoreIndices(part,&indices);
560:   ISCreateGeneral(comm,n,newi,PETSC_OWN_POINTER,is);
561:   ISSetPermutation(*is);
562:   return 0;
563: }

565: /*@
566:     ISPartitioningCount - Takes a ISPartitioning and determines the number of
567:     resulting elements on each (partition) process

569:     Collective on IS

571:     Input Parameters:
572: +   partitioning - a partitioning as generated by MatPartitioningApply() or
573:                    MatPartitioningApplyND()
574: -   len - length of the array count, this is the total number of partitions

576:     Output Parameter:
577: .   count - array of length size, to contain the number of elements assigned
578:         to each partition, where size is the number of partitions generated
579:          (see notes below).

581:    Level: advanced

583:     Notes:
584:         By default the number of partitions generated (and thus the length
585:         of count) is the size of the communicator associated with IS,
586:         but it can be set by MatPartitioningSetNParts. The resulting array
587:         of lengths can for instance serve as input of PCBJacobiSetTotalBlocks.
588:         If the partitioning has been obtained by MatPartitioningApplyND(),
589:         the returned count does not include the separators.

591: .seealso: MatPartitioningCreate(), AOCreateBasic(), ISPartitioningToNumbering(),
592:         MatPartitioningSetNParts(), MatPartitioningApply(), MatPartitioningApplyND()

594: @*/
595: PetscErrorCode  ISPartitioningCount(IS part,PetscInt len,PetscInt count[])
596: {
597:   MPI_Comm       comm;
598:   PetscInt       i,n,*lsizes;
599:   const PetscInt *indices;
600:   PetscMPIInt    npp;

602:   PetscObjectGetComm((PetscObject)part,&comm);
603:   if (len == PETSC_DEFAULT) {
604:     PetscMPIInt size;
605:     MPI_Comm_size(comm,&size);
606:     len  = (PetscInt) size;
607:   }

609:   /* count the number of partitions */
610:   ISGetLocalSize(part,&n);
611:   ISGetIndices(part,&indices);
612:   if (PetscDefined(USE_DEBUG)) {
613:     PetscInt np = 0,npt;
614:     for (i=0; i<n; i++) np = PetscMax(np,indices[i]);
615:     MPIU_Allreduce(&np,&npt,1,MPIU_INT,MPI_MAX,comm);
616:     np   = npt+1; /* so that it looks like a MPI_Comm_size output */
618:   }

620:   /*
621:         lsizes - number of elements of each partition on this particular processor
622:         sums - total number of "previous" nodes for any particular partition
623:         starts - global number of first element in each partition on this processor
624:   */
625:   PetscCalloc1(len,&lsizes);
626:   for (i=0; i<n; i++) {
627:     if (indices[i] > -1) lsizes[indices[i]]++;
628:   }
629:   ISRestoreIndices(part,&indices);
630:   PetscMPIIntCast(len,&npp);
631:   MPIU_Allreduce(lsizes,count,npp,MPIU_INT,MPI_SUM,comm);
632:   PetscFree(lsizes);
633:   return 0;
634: }

636: /*@
637:     ISAllGather - Given an index set (IS) on each processor, generates a large
638:     index set (same on each processor) by concatenating together each
639:     processors index set.

641:     Collective on IS

643:     Input Parameter:
644: .   is - the distributed index set

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

649:     Notes:
650:     ISAllGather() is clearly not scalable for large index sets.

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

658:     The communicator for this new IS is PETSC_COMM_SELF

660:     Level: intermediate

662: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock()
663: @*/
664: PetscErrorCode  ISAllGather(IS is,IS *isout)
665: {
666:   PetscInt       *indices,n,i,N,step,first;
667:   const PetscInt *lindices;
668:   MPI_Comm       comm;
669:   PetscMPIInt    size,*sizes = NULL,*offsets = NULL,nn;
670:   PetscBool      stride;


675:   PetscObjectGetComm((PetscObject)is,&comm);
676:   MPI_Comm_size(comm,&size);
677:   ISGetLocalSize(is,&n);
678:   PetscObjectTypeCompare((PetscObject)is,ISSTRIDE,&stride);
679:   if (size == 1 && stride) { /* should handle parallel ISStride also */
680:     ISStrideGetInfo(is,&first,&step);
681:     ISCreateStride(PETSC_COMM_SELF,n,first,step,isout);
682:   } else {
683:     PetscMalloc2(size,&sizes,size,&offsets);

685:     PetscMPIIntCast(n,&nn);
686:     MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
687:     offsets[0] = 0;
688:     for (i=1; i<size; i++) {
689:       PetscInt s = offsets[i-1] + sizes[i-1];
690:       PetscMPIIntCast(s,&offsets[i]);
691:     }
692:     N = offsets[size-1] + sizes[size-1];

694:     PetscMalloc1(N,&indices);
695:     ISGetIndices(is,&lindices);
696:     MPI_Allgatherv((void*)lindices,nn,MPIU_INT,indices,sizes,offsets,MPIU_INT,comm);
697:     ISRestoreIndices(is,&lindices);
698:     PetscFree2(sizes,offsets);

700:     ISCreateGeneral(PETSC_COMM_SELF,N,indices,PETSC_OWN_POINTER,isout);
701:   }
702:   return 0;
703: }

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

709:     Collective

711:     Input Parameters:
712: +   comm - communicator to share the indices
713: .   n - local size of set
714: -   lindices - local colors

716:     Output Parameters:
717: +   outN - total number of indices
718: -   outindices - all of the colors

720:     Notes:
721:     ISAllGatherColors() is clearly not scalable for large index sets.

723:     Level: intermediate

725: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
726: @*/
727: PetscErrorCode  ISAllGatherColors(MPI_Comm comm,PetscInt n,ISColoringValue *lindices,PetscInt *outN,ISColoringValue *outindices[])
728: {
729:   ISColoringValue *indices;
730:   PetscInt        i,N;
731:   PetscMPIInt     size,*offsets = NULL,*sizes = NULL, nn = n;

733:   MPI_Comm_size(comm,&size);
734:   PetscMalloc2(size,&sizes,size,&offsets);

736:   MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
737:   offsets[0] = 0;
738:   for (i=1; i<size; i++) offsets[i] = offsets[i-1] + sizes[i-1];
739:   N    = offsets[size-1] + sizes[size-1];
740:   PetscFree2(sizes,offsets);

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

745:   *outindices = indices;
746:   if (outN) *outN = N;
747:   return 0;
748: }

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

754:     Collective on IS

756:     Input Parameters:
757: +   is - the index set
758: .   nmin - the first index desired in the local part of the complement
759: -   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)

761:     Output Parameter:
762: .   isout - the complement

764:     Notes:
765:     The communicator for this new IS is the same as for the input IS

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

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

772:     Level: intermediate

774: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
775: @*/
776: PetscErrorCode  ISComplement(IS is,PetscInt nmin,PetscInt nmax,IS *isout)
777: {
778:   const PetscInt *indices;
779:   PetscInt       n,i,j,unique,cnt,*nindices;
780:   PetscBool      sorted;

786:   ISSorted(is,&sorted);

789:   ISGetLocalSize(is,&n);
790:   ISGetIndices(is,&indices);
791:   if (PetscDefined(USE_DEBUG)) {
792:     for (i=0; i<n; i++) {
795:     }
796:   }
797:   /* Count number of unique entries */
798:   unique = (n>0);
799:   for (i=0; i<n-1; i++) {
800:     if (indices[i+1] != indices[i]) unique++;
801:   }
802:   PetscMalloc1(nmax-nmin-unique,&nindices);
803:   cnt  = 0;
804:   for (i=nmin,j=0; i<nmax; i++) {
805:     if (j<n && i==indices[j]) do { j++; } while (j<n && i==indices[j]);
806:     else nindices[cnt++] = i;
807:   }
809:   ISCreateGeneral(PetscObjectComm((PetscObject)is),cnt,nindices,PETSC_OWN_POINTER,isout);
810:   ISRestoreIndices(is,&indices);
811:   return 0;
812: }