Actual source code: iscoloring.c

petsc-master 2019-11-11
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: /*@C

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

 19:    Collective on coloring

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

 25:    Notes:
 26:      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

 28:    Level: intermediate

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

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

 40: /*@C

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

 44:    Collective on coloring

 46:    Input Parameter:
 47: .   coloring - the coloring object

 49:    Output Parameter:
 50: .    type - either IS_COLORING_LOCAL or IS_COLORING_GLOBAL

 52:    Level: intermediate

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

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

 64: /*@
 65:    ISColoringDestroy - Destroys a coloring context.

 67:    Collective on ISColoring

 69:    Input Parameter:
 70: .  iscoloring - the coloring context

 72:    Level: advanced

 74: .seealso: ISColoringView(), MatColoring
 75: @*/
 76: PetscErrorCode  ISColoringDestroy(ISColoring *iscoloring)
 77: {
 78:   PetscInt       i;

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

 86:   if ((*iscoloring)->is) {
 87:     for (i=0; i<(*iscoloring)->n; i++) {
 88:       ISDestroy(&(*iscoloring)->is[i]);
 89:     }
 90:     PetscFree((*iscoloring)->is);
 91:   }
 92:   if ((*iscoloring)->allocated) {PetscFree((*iscoloring)->colors);}
 93:   PetscCommDestroy(&(*iscoloring)->comm);
 94:   PetscFree((*iscoloring));
 95:   return(0);
 96: }

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

101:   Collective on ISColoring

103:   Input Parameters:
104: + obj   - the ISColoring object
105: . prefix - prefix to use for viewing, or NULL to use prefix of 'mat'
106: - optionname - option to activate viewing

108:   Level: intermediate

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

112: */
113: PetscErrorCode ISColoringViewFromOptions(ISColoring obj,PetscObject bobj,const char optionname[])
114: {
115:   PetscErrorCode    ierr;
116:   PetscViewer       viewer;
117:   PetscBool         flg;
118:   PetscViewerFormat format;
119:   char              *prefix;

122:   prefix = bobj ? bobj->prefix : NULL;
123:   PetscOptionsGetViewer(obj->comm,NULL,prefix,optionname,&viewer,&format,&flg);
124:   if (flg) {
125:     PetscViewerPushFormat(viewer,format);
126:     ISColoringView(obj,viewer);
127:     PetscViewerPopFormat(viewer);
128:     PetscViewerDestroy(&viewer);
129:   }
130:   return(0);
131: }

133: /*@C
134:    ISColoringView - Views a coloring context.

136:    Collective on ISColoring

138:    Input Parameters:
139: +  iscoloring - the coloring context
140: -  viewer - the viewer

142:    Level: advanced

144: .seealso: ISColoringDestroy(), ISColoringGetIS(), MatColoring
145: @*/
146: PetscErrorCode  ISColoringView(ISColoring iscoloring,PetscViewer viewer)
147: {
148:   PetscInt       i;
150:   PetscBool      iascii;
151:   IS             *is;

155:   if (!viewer) {
156:     PetscViewerASCIIGetStdout(iscoloring->comm,&viewer);
157:   }

160:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
161:   if (iascii) {
162:     MPI_Comm    comm;
163:     PetscMPIInt size,rank;

165:     PetscObjectGetComm((PetscObject)viewer,&comm);
166:     MPI_Comm_size(comm,&size);
167:     MPI_Comm_rank(comm,&rank);
168:     PetscViewerASCIIPrintf(viewer,"ISColoring Object: %d MPI processes\n",size);
169:     PetscViewerASCIIPrintf(viewer,"ISColoringType: %s\n",ISColoringTypes[iscoloring->ctype]);
170:     PetscViewerASCIIPushSynchronized(viewer);
171:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of colors %d\n",rank,iscoloring->n);
172:     PetscViewerFlush(viewer);
173:     PetscViewerASCIIPopSynchronized(viewer);
174:   }

176:   ISColoringGetIS(iscoloring,PETSC_USE_POINTER,PETSC_IGNORE,&is);
177:   for (i=0; i<iscoloring->n; i++) {
178:     ISView(iscoloring->is[i],viewer);
179:   }
180:   ISColoringRestoreIS(iscoloring,PETSC_USE_POINTER,&is);
181:   return(0);
182: }

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

187:    Not Collective

189:    Input Parameter:
190: .  iscoloring - the coloring context

192:    Output Parameters:
193: +  n - number of nodes
194: .  nc - number of colors
195: -  colors - color for each node

197:    Level: advanced

199: .seealso: ISColoringRestoreIS(), ISColoringView(), ISColoringGetIS()
200: @*/
201: PetscErrorCode  ISColoringGetColors(ISColoring iscoloring,PetscInt *n,PetscInt *nc,const ISColoringValue **colors)
202: {

206:   if (n) *n = iscoloring->N;
207:   if (nc) *nc = iscoloring->n;
208:   if (colors) *colors = iscoloring->colors;
209:   return(0);
210: }

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

215:    Collective on ISColoring

217:    Input Parameter:
218: +  iscoloring - the coloring context
219: -  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

221:    Output Parameters:
222: +  nn - number of index sets in the coloring context
223: -  is - array of index sets

225:    Level: advanced

227: .seealso: ISColoringRestoreIS(), ISColoringView(), ISColoringGetColoring()
228: @*/
229: PetscErrorCode  ISColoringGetIS(ISColoring iscoloring,PetscCopyMode mode, PetscInt *nn,IS *isis[])
230: {


236:   if (nn) *nn = iscoloring->n;
237:   if (isis) {
238:     if (!iscoloring->is) {
239:       PetscInt        *mcolors,**ii,nc = iscoloring->n,i,base, n = iscoloring->N;
240:       ISColoringValue *colors = iscoloring->colors;
241:       IS              *is;

243: #if defined(PETSC_USE_DEBUG)
244:       for (i=0; i<n; i++) {
245:         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);
246:       }
247: #endif

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

253:       PetscMalloc1(nc,&ii);
254:       PetscMalloc1(n,&ii[0]);
255:       for (i=1; i<nc; i++) ii[i] = ii[i-1] + mcolors[i-1];
256:       PetscArrayzero(mcolors,nc);

258:       if (iscoloring->ctype == IS_COLORING_GLOBAL) {
259:         MPI_Scan(&iscoloring->N,&base,1,MPIU_INT,MPI_SUM,iscoloring->comm);
260:         base -= iscoloring->N;
261:         for (i=0; i<n; i++) ii[colors[i]][mcolors[colors[i]]++] = i + base; /* global idx */
262:       } else if (iscoloring->ctype == IS_COLORING_LOCAL) {
263:         for (i=0; i<n; i++) ii[colors[i]][mcolors[colors[i]]++] = i;   /* local idx */
264:       } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not provided for this ISColoringType type");

266:       PetscMalloc1(nc,&is);
267:       for (i=0; i<nc; i++) {
268:         ISCreateGeneral(iscoloring->comm,mcolors[i],ii[i],PETSC_COPY_VALUES,is+i);
269:       }

271:       if (mode != PETSC_OWN_POINTER) iscoloring->is = is;
272:       *isis = is;
273:       PetscFree(ii[0]);
274:       PetscFree(ii);
275:       PetscFree(mcolors);
276:     } else {
277:       *isis = iscoloring->is;
278:     }
279:   }
280:   return(0);
281: }

283: /*@C
284:    ISColoringRestoreIS - Restores the index sets extracted from the coloring context

286:    Collective on ISColoring

288:    Input Parameter:
289: +  iscoloring - the coloring context
290: .  mode - who retains ownership of the is
291: -  is - array of index sets

293:    Level: advanced

295: .seealso: ISColoringGetIS(), ISColoringView()
296: @*/
297: PetscErrorCode  ISColoringRestoreIS(ISColoring iscoloring,PetscCopyMode mode,IS *is[])
298: {

302:   /* currently nothing is done here */
303:   return(0);
304: }


307: /*@
308:     ISColoringCreate - Generates an ISColoring context from lists (provided
309:     by each processor) of colors for each node.

311:     Collective

313:     Input Parameters:
314: +   comm - communicator for the processors creating the coloring
315: .   ncolors - max color value
316: .   n - number of nodes on this processor
317: .   colors - array containing the colors for this processor, color numbers begin at 0.
318: -   mode - see PetscCopyMode for meaning of this flag.

320:     Output Parameter:
321: .   iscoloring - the resulting coloring data structure

323:     Options Database Key:
324: .   -is_coloring_view - Activates ISColoringView()

326:    Level: advanced

328:     Notes:
329:     By default sets coloring type to  IS_COLORING_GLOBAL

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

333: @*/
334: PetscErrorCode  ISColoringCreate(MPI_Comm comm,PetscInt ncolors,PetscInt n,const ISColoringValue colors[],PetscCopyMode mode,ISColoring *iscoloring)
335: {
337:   PetscMPIInt    size,rank,tag;
338:   PetscInt       base,top,i;
339:   PetscInt       nc,ncwork;
340:   MPI_Status     status;

343:   if (ncolors != PETSC_DECIDE && ncolors > IS_COLORING_MAX) {
344:     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);
345:     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);
346:   }
347:   PetscNew(iscoloring);
348:   PetscCommDuplicate(comm,&(*iscoloring)->comm,&tag);
349:   comm = (*iscoloring)->comm;

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

354:   /* should use MPI_Scan() */
355:   MPI_Comm_rank(comm,&rank);
356:   if (!rank) {
357:     base = 0;
358:     top  = n;
359:   } else {
360:     MPI_Recv(&base,1,MPIU_INT,rank-1,tag,comm,&status);
361:     top  = base+n;
362:   }
363:   if (rank < size-1) {
364:     MPI_Send(&top,1,MPIU_INT,rank+1,tag,comm);
365:   }

367:   /* compute the total number of colors */
368:   ncwork = 0;
369:   for (i=0; i<n; i++) {
370:     if (ncwork < colors[i]) ncwork = colors[i];
371:   }
372:   ncwork++;
373:   MPIU_Allreduce(&ncwork,&nc,1,MPIU_INT,MPI_MAX,comm);
374:   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);
375:   (*iscoloring)->n      = nc;
376:   (*iscoloring)->is     = 0;
377:   (*iscoloring)->N      = n;
378:   (*iscoloring)->refct  = 1;
379:   (*iscoloring)->ctype  = IS_COLORING_GLOBAL;
380:   if (mode == PETSC_COPY_VALUES) {
381:     PetscMalloc1(n,&(*iscoloring)->colors);
382:     PetscLogObjectMemory((PetscObject)(*iscoloring),n*sizeof(ISColoringValue));
383:     PetscArraycpy((*iscoloring)->colors,colors,n);
384:     (*iscoloring)->allocated = PETSC_TRUE;
385:   } else if (mode == PETSC_OWN_POINTER) {
386:     (*iscoloring)->colors    = (ISColoringValue*)colors;
387:     (*iscoloring)->allocated = PETSC_TRUE;
388:   } else {
389:     (*iscoloring)->colors    = (ISColoringValue*)colors;
390:     (*iscoloring)->allocated = PETSC_FALSE;
391:   }
392:   ISColoringViewFromOptions(*iscoloring,NULL,"-is_coloring_view");
393:   PetscInfo1(0,"Number of colors %D\n",nc);
394:   return(0);
395: }

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

401:     Collective on IS

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

407:     Output Parameter:
408: .   rows - contains new numbers from remote or local

410:    Level: advanced

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

414: @*/
415: PetscErrorCode  ISBuildTwoSided(IS ito,IS toindx, IS *rows)
416: {
417:    const PetscInt *ito_indices,*toindx_indices;
418:    PetscInt       *send_indices,rstart,*recv_indices,nrecvs,nsends;
419:    PetscInt       *tosizes,*fromsizes,i,j,*tosizes_tmp,*tooffsets_tmp,ito_ln;
420:    PetscMPIInt    *toranks,*fromranks,size,target_rank,*fromperm_newtoold,nto,nfrom;
421:    PetscLayout     isrmap;
422:    MPI_Comm        comm;
423:    PetscSF         sf;
424:    PetscSFNode    *iremote;
425:    PetscErrorCode  ierr;

428:    PetscObjectGetComm((PetscObject)ito,&comm);
429:    MPI_Comm_size(comm,&size);
430:    ISGetLocalSize(ito,&ito_ln);
431:    /* why we do not have ISGetLayout? */
432:    isrmap = ito->map;
433:    PetscLayoutGetRange(isrmap,&rstart,NULL);
434:    ISGetIndices(ito,&ito_indices);
435:    PetscCalloc2(size,&tosizes_tmp,size+1,&tooffsets_tmp);
436:    for (i=0; i<ito_ln; i++) {
437:      if (ito_indices[i]<0) continue;
438: #if defined(PETSC_USE_DEBUG)
439:      if (ito_indices[i]>=size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"target rank %d is larger than communicator size %d ",ito_indices[i],size);
440: #endif
441:      tosizes_tmp[ito_indices[i]]++;
442:    }
443:    nto = 0;
444:    for (i=0; i<size; i++) {
445:      tooffsets_tmp[i+1] = tooffsets_tmp[i]+tosizes_tmp[i];
446:      if (tosizes_tmp[i]>0) nto++;
447:    }
448:    PetscCalloc2(nto,&toranks,2*nto,&tosizes);
449:    nto  = 0;
450:    for (i=0; i<size; i++) {
451:      if (tosizes_tmp[i]>0) {
452:        toranks[nto]     = i;
453:        tosizes[2*nto]   = tosizes_tmp[i];/* size */
454:        tosizes[2*nto+1] = tooffsets_tmp[i];/* offset */
455:        nto++;
456:      }
457:    }
458:    nsends = tooffsets_tmp[size];
459:    PetscCalloc1(nsends,&send_indices);
460:    if (toindx) {
461:      ISGetIndices(toindx,&toindx_indices);
462:    }
463:    for (i=0; i<ito_ln; i++) {
464:      if (ito_indices[i]<0) continue;
465:      target_rank = ito_indices[i];
466:      send_indices[tooffsets_tmp[target_rank]] = toindx? toindx_indices[i]:(i+rstart);
467:      tooffsets_tmp[target_rank]++;
468:    }
469:    if (toindx) {
470:      ISRestoreIndices(toindx,&toindx_indices);
471:    }
472:    ISRestoreIndices(ito,&ito_indices);
473:    PetscFree2(tosizes_tmp,tooffsets_tmp);
474:    PetscCommBuildTwoSided(comm,2,MPIU_INT,nto,toranks,tosizes,&nfrom,&fromranks,&fromsizes);
475:    PetscFree2(toranks,tosizes);
476:    PetscMalloc1(nfrom,&fromperm_newtoold);
477:    for (i=0; i<nfrom; i++) fromperm_newtoold[i] = i;
478:    PetscSortMPIIntWithArray(nfrom,fromranks,fromperm_newtoold);
479:    nrecvs = 0;
480:    for (i=0; i<nfrom; i++) nrecvs += fromsizes[i*2];
481:    PetscCalloc1(nrecvs,&recv_indices);
482:    PetscMalloc1(nrecvs,&iremote);
483:    nrecvs = 0;
484:    for (i=0; i<nfrom; i++) {
485:      for (j=0; j<fromsizes[2*fromperm_newtoold[i]]; j++) {
486:        iremote[nrecvs].rank    = fromranks[i];
487:        iremote[nrecvs++].index = fromsizes[2*fromperm_newtoold[i]+1]+j;
488:      }
489:    }
490:    PetscSFCreate(comm,&sf);
491:    PetscSFSetGraph(sf,nsends,nrecvs,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
492:    PetscSFSetType(sf,PETSCSFBASIC);
493:    /* how to put a prefix ? */
494:    PetscSFSetFromOptions(sf);
495:    PetscSFBcastBegin(sf,MPIU_INT,send_indices,recv_indices);
496:    PetscSFBcastEnd(sf,MPIU_INT,send_indices,recv_indices);
497:    PetscSFDestroy(&sf);
498:    PetscFree(fromranks);
499:    PetscFree(fromsizes);
500:    PetscFree(fromperm_newtoold);
501:    PetscFree(send_indices);
502:    if (rows) {
503:      PetscSortInt(nrecvs,recv_indices);
504:      ISCreateGeneral(comm,nrecvs,recv_indices,PETSC_OWN_POINTER,rows);
505:    } else {
506:      PetscFree(recv_indices);
507:    }
508:    return(0);
509: }


512: /*@
513:     ISPartitioningToNumbering - Takes an ISPartitioning and on each processor
514:     generates an IS that contains a new global node number for each index based
515:     on the partitioing.

517:     Collective on IS

519:     Input Parameters
520: .   partitioning - a partitioning as generated by MatPartitioningApply()
521:                    or MatPartitioningApplyND()

523:     Output Parameter:
524: .   is - on each processor the index set that defines the global numbers
525:          (in the new numbering) for all the nodes currently (before the partitioning)
526:          on that processor

528:    Level: advanced

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

532: @*/
533: PetscErrorCode  ISPartitioningToNumbering(IS part,IS *is)
534: {
535:   MPI_Comm       comm;
536:   IS             ndorder;
537:   PetscInt       i,np,npt,n,*starts = NULL,*sums = NULL,*lsizes = NULL,*newi = NULL;
538:   const PetscInt *indices = NULL;

544:   /* see if the partitioning comes from nested dissection */
545:   PetscObjectQuery((PetscObject)part,"_petsc_matpartitioning_ndorder",(PetscObject*)&ndorder);
546:   if (ndorder) {
547:     PetscObjectReference((PetscObject)ndorder);
548:     *is  = ndorder;
549:     return(0);
550:   }

552:   PetscObjectGetComm((PetscObject)part,&comm);
553:   /* count the number of partitions, i.e., virtual processors */
554:   ISGetLocalSize(part,&n);
555:   ISGetIndices(part,&indices);
556:   np   = 0;
557:   for (i=0; i<n; i++) np = PetscMax(np,indices[i]);
558:   MPIU_Allreduce(&np,&npt,1,MPIU_INT,MPI_MAX,comm);
559:   np   = npt+1; /* so that it looks like a MPI_Comm_size output */

561:   /*
562:         lsizes - number of elements of each partition on this particular processor
563:         sums - total number of "previous" nodes for any particular partition
564:         starts - global number of first element in each partition on this processor
565:   */
566:   PetscMalloc3(np,&lsizes,np,&starts,np,&sums);
567:   PetscArrayzero(lsizes,np);
568:   for (i=0; i<n; i++) lsizes[indices[i]]++;
569:   MPIU_Allreduce(lsizes,sums,np,MPIU_INT,MPI_SUM,comm);
570:   MPI_Scan(lsizes,starts,np,MPIU_INT,MPI_SUM,comm);
571:   for (i=0; i<np; i++) starts[i] -= lsizes[i];
572:   for (i=1; i<np; i++) {
573:     sums[i]   += sums[i-1];
574:     starts[i] += sums[i-1];
575:   }

577:   /*
578:       For each local index give it the new global number
579:   */
580:   PetscMalloc1(n,&newi);
581:   for (i=0; i<n; i++) newi[i] = starts[indices[i]]++;
582:   PetscFree3(lsizes,starts,sums);

584:   ISRestoreIndices(part,&indices);
585:   ISCreateGeneral(comm,n,newi,PETSC_OWN_POINTER,is);
586:   ISSetPermutation(*is);
587:   return(0);
588: }

590: /*@
591:     ISPartitioningCount - Takes a ISPartitioning and determines the number of
592:     resulting elements on each (partition) process

594:     Collective on IS

596:     Input Parameters:
597: +   partitioning - a partitioning as generated by MatPartitioningApply() or
598:                    MatPartitioningApplyND()
599: -   len - length of the array count, this is the total number of partitions

601:     Output Parameter:
602: .   count - array of length size, to contain the number of elements assigned
603:         to each partition, where size is the number of partitions generated
604:          (see notes below).

606:    Level: advanced

608:     Notes:
609:         By default the number of partitions generated (and thus the length
610:         of count) is the size of the communicator associated with IS,
611:         but it can be set by MatPartitioningSetNParts. The resulting array
612:         of lengths can for instance serve as input of PCBJacobiSetTotalBlocks.
613:         If the partitioning has been obtained by MatPartitioningApplyND(),
614:         the returned count does not include the separators.

616: .seealso: MatPartitioningCreate(), AOCreateBasic(), ISPartitioningToNumbering(),
617:         MatPartitioningSetNParts(), MatPartitioningApply(), MatPartitioningApplyND()

619: @*/
620: PetscErrorCode  ISPartitioningCount(IS part,PetscInt len,PetscInt count[])
621: {
622:   MPI_Comm       comm;
623:   PetscInt       i,n,*lsizes;
624:   const PetscInt *indices;
626:   PetscMPIInt    npp;

629:   PetscObjectGetComm((PetscObject)part,&comm);
630:   if (len == PETSC_DEFAULT) {
631:     PetscMPIInt size;
632:     MPI_Comm_size(comm,&size);
633:     len  = (PetscInt) size;
634:   }

636:   /* count the number of partitions */
637:   ISGetLocalSize(part,&n);
638:   ISGetIndices(part,&indices);
639: #if defined(PETSC_USE_DEBUG)
640:   {
641:     PetscInt np = 0,npt;
642:     for (i=0; i<n; i++) np = PetscMax(np,indices[i]);
643:     MPIU_Allreduce(&np,&npt,1,MPIU_INT,MPI_MAX,comm);
644:     np   = npt+1; /* so that it looks like a MPI_Comm_size output */
645:     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);
646:   }
647: #endif

649:   /*
650:         lsizes - number of elements of each partition on this particular processor
651:         sums - total number of "previous" nodes for any particular partition
652:         starts - global number of first element in each partition on this processor
653:   */
654:   PetscCalloc1(len,&lsizes);
655:   for (i=0; i<n; i++) {
656:     if (indices[i] > -1) lsizes[indices[i]]++;
657:   }
658:   ISRestoreIndices(part,&indices);
659:   PetscMPIIntCast(len,&npp);
660:   MPIU_Allreduce(lsizes,count,npp,MPIU_INT,MPI_SUM,comm);
661:   PetscFree(lsizes);
662:   return(0);
663: }

665: /*@
666:     ISAllGather - Given an index set (IS) on each processor, generates a large
667:     index set (same on each processor) by concatenating together each
668:     processors index set.

670:     Collective on IS

672:     Input Parameter:
673: .   is - the distributed index set

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

678:     Notes:
679:     ISAllGather() is clearly not scalable for large index sets.

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

687:     The communicator for this new IS is PETSC_COMM_SELF

689:     Level: intermediate

691: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock()
692: @*/
693: PetscErrorCode  ISAllGather(IS is,IS *isout)
694: {
696:   PetscInt       *indices,n,i,N,step,first;
697:   const PetscInt *lindices;
698:   MPI_Comm       comm;
699:   PetscMPIInt    size,*sizes = NULL,*offsets = NULL,nn;
700:   PetscBool      stride;


706:   PetscObjectGetComm((PetscObject)is,&comm);
707:   MPI_Comm_size(comm,&size);
708:   ISGetLocalSize(is,&n);
709:   PetscObjectTypeCompare((PetscObject)is,ISSTRIDE,&stride);
710:   if (size == 1 && stride) { /* should handle parallel ISStride also */
711:     ISStrideGetInfo(is,&first,&step);
712:     ISCreateStride(PETSC_COMM_SELF,n,first,step,isout);
713:   } else {
714:     PetscMalloc2(size,&sizes,size,&offsets);

716:     PetscMPIIntCast(n,&nn);
717:     MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
718:     offsets[0] = 0;
719:     for (i=1; i<size; i++) {
720:       PetscInt s = offsets[i-1] + sizes[i-1];
721:       PetscMPIIntCast(s,&offsets[i]);
722:     }
723:     N = offsets[size-1] + sizes[size-1];

725:     PetscMalloc1(N,&indices);
726:     ISGetIndices(is,&lindices);
727:     MPI_Allgatherv((void*)lindices,nn,MPIU_INT,indices,sizes,offsets,MPIU_INT,comm);
728:     ISRestoreIndices(is,&lindices);
729:     PetscFree2(sizes,offsets);

731:     ISCreateGeneral(PETSC_COMM_SELF,N,indices,PETSC_OWN_POINTER,isout);
732:   }
733:   return(0);
734: }

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

740:     Collective

742:     Input Parameter:
743: +   comm - communicator to share the indices
744: .   n - local size of set
745: -   lindices - local colors

747:     Output Parameter:
748: +   outN - total number of indices
749: -   outindices - all of the colors

751:     Notes:
752:     ISAllGatherColors() is clearly not scalable for large index sets.


755:     Level: intermediate

757: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
758: @*/
759: PetscErrorCode  ISAllGatherColors(MPI_Comm comm,PetscInt n,ISColoringValue *lindices,PetscInt *outN,ISColoringValue *outindices[])
760: {
761:   ISColoringValue *indices;
762:   PetscErrorCode  ierr;
763:   PetscInt        i,N;
764:   PetscMPIInt     size,*offsets = NULL,*sizes = NULL, nn = n;

767:   MPI_Comm_size(comm,&size);
768:   PetscMalloc2(size,&sizes,size,&offsets);

770:   MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
771:   offsets[0] = 0;
772:   for (i=1; i<size; i++) offsets[i] = offsets[i-1] + sizes[i-1];
773:   N    = offsets[size-1] + sizes[size-1];
774:   PetscFree2(sizes,offsets);

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

779:   *outindices = indices;
780:   if (outN) *outN = N;
781:   return(0);
782: }

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

788:     Collective on IS

790:     Input Parameter:
791: +   is - the index set
792: .   nmin - the first index desired in the local part of the complement
793: -   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)

795:     Output Parameter:
796: .   isout - the complement

798:     Notes:
799:     The communicator for this new IS is the same as for the input IS

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

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

806:     Level: intermediate

808: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
809: @*/
810: PetscErrorCode  ISComplement(IS is,PetscInt nmin,PetscInt nmax,IS *isout)
811: {
813:   const PetscInt *indices;
814:   PetscInt       n,i,j,unique,cnt,*nindices;
815:   PetscBool      sorted;

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

825:   ISGetLocalSize(is,&n);
826:   ISGetIndices(is,&indices);
827: #if defined(PETSC_USE_DEBUG)
828:   for (i=0; i<n; i++) {
829:     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);
830:     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);
831:   }
832: #endif
833:   /* Count number of unique entries */
834:   unique = (n>0);
835:   for (i=0; i<n-1; i++) {
836:     if (indices[i+1] != indices[i]) unique++;
837:   }
838:   PetscMalloc1(nmax-nmin-unique,&nindices);
839:   cnt  = 0;
840:   for (i=nmin,j=0; i<nmax; i++) {
841:     if (j<n && i==indices[j]) do { j++; } while (j<n && i==indices[j]);
842:     else nindices[cnt++] = i;
843:   }
844:   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);
845:   ISCreateGeneral(PetscObjectComm((PetscObject)is),cnt,nindices,PETSC_OWN_POINTER,isout);
846:   ISRestoreIndices(is,&indices);
847:   return(0);
848: }