Actual source code: iscoloring.c

      1: #define PETSCVEC_DLL

     3:  #include petscsys.h
     4:  #include petscis.h

      8: /*@C
      9:    ISColoringDestroy - Destroys a coloring context.

     11:    Collective on ISColoring

     13:    Input Parameter:
     14: .  iscoloring - the coloring context

     16:    Level: advanced

     18: .seealso: ISColoringView(), MatGetColoring()
     19: @*/
     20: PetscErrorCode PETSCVEC_DLLEXPORT ISColoringDestroy(ISColoring iscoloring)
     21: {
     22:   PetscInt i;

     27:   if (--iscoloring->refct > 0) return(0);

     29:   if (iscoloring->is) {
     30:     for (i=0; i<iscoloring->n; i++) {
     31:       ISDestroy(iscoloring->is[i]);
     32:     }
     33:     PetscFree(iscoloring->is);
     34:   }
     35:   if (iscoloring->colors) {
     36:     PetscFree(iscoloring->colors);
     37:   }
     38:   PetscCommDestroy(&iscoloring->comm);
     39:   PetscFree(iscoloring);
     40:   return(0);
     41: }

     45: /*@C
     46:    ISColoringView - Views a coloring context.

     48:    Collective on ISColoring

     50:    Input Parameters:
     51: +  iscoloring - the coloring context
     52: -  viewer - the viewer

     54:    Level: advanced

     56: .seealso: ISColoringDestroy(), ISColoringGetIS(), MatGetColoring()
     57: @*/
     58: PetscErrorCode PETSCVEC_DLLEXPORT ISColoringView(ISColoring iscoloring,PetscViewer viewer)
     59: {
     60:   PetscInt       i;
     62:   PetscTruth     iascii;
     63:   IS             *is;

     67:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(iscoloring->comm);

     70:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
     71:   if (iascii) {
     72:     MPI_Comm    comm;
     73:     PetscMPIInt rank;
     74:     PetscObjectGetComm((PetscObject)viewer,&comm);
     75:     MPI_Comm_rank(comm,&rank);
     76:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of colors %d\n",rank,iscoloring->n);
     77:     PetscViewerFlush(viewer);
     78:   } else {
     79:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for ISColoring",((PetscObject)viewer)->type_name);
     80:   }

     82:   ISColoringGetIS(iscoloring,PETSC_IGNORE,&is);
     83:   for (i=0; i<iscoloring->n; i++) {
     84:     ISView(iscoloring->is[i],viewer);
     85:   }
     86:   ISColoringRestoreIS(iscoloring,&is);
     87:   return(0);
     88: }

     92: /*@C
     93:    ISColoringGetIS - Extracts index sets from the coloring context

     95:    Collective on ISColoring

     97:    Input Parameter:
     98: .  iscoloring - the coloring context

    100:    Output Parameters:
    101: +  nn - number of index sets in the coloring context
    102: -  is - array of index sets

    104:    Level: advanced

    106: .seealso: ISColoringRestoreIS(), ISColoringView()
    107: @*/
    108: PetscErrorCode PETSCVEC_DLLEXPORT ISColoringGetIS(ISColoring iscoloring,PetscInt *nn,IS *isis[])
    109: {

    115:   if (nn)  *nn  = iscoloring->n;
    116:   if (isis) {
    117:     if (!iscoloring->is) {
    118:       PetscInt        *mcolors,**ii,nc = iscoloring->n,i,base, n = iscoloring->N;
    119:       ISColoringValue *colors = iscoloring->colors;
    120:       IS              *is;

    122: #if defined(PETSC_USE_DEBUG)
    123:       for (i=0; i<n; i++) {
    124:         if (colors[i] >= (ISColoringValue) nc) {
    125:           SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"Coloring is our of range index %d value %d number colors %d",(int)i,(int)colors[i],(int)nc);
    126:         }
    127:       }
    128: #endif
    129:
    130:       /* generate the lists of nodes for each color */
    131:       PetscMalloc(nc*sizeof(PetscInt),&mcolors);
    132:       PetscMemzero(mcolors,nc*sizeof(PetscInt));
    133:       for (i=0; i<n; i++) {
    134:         mcolors[colors[i]]++;
    135:       }

    137:       PetscMalloc(nc*sizeof(PetscInt*),&ii);
    138:       PetscMalloc(n*sizeof(PetscInt),&ii[0]);
    139:       for (i=1; i<nc; i++) {
    140:         ii[i] = ii[i-1] + mcolors[i-1];
    141:       }
    142:
    143:       MPI_Scan(&iscoloring->N,&base,1,MPIU_INT,MPI_SUM,iscoloring->comm);
    144:       base -= iscoloring->N;
    145:       PetscMemzero(mcolors,nc*sizeof(PetscInt));
    146:       for (i=0; i<n; i++) {
    147:         ii[colors[i]][mcolors[colors[i]]++] = i + base;
    148:       }
    149:       PetscMalloc(nc*sizeof(IS),&is);
    150:       for (i=0; i<nc; i++) {
    151:         ISCreateGeneral(iscoloring->comm,mcolors[i],ii[i],is+i);
    152:       }

    154:       iscoloring->is   = is;
    155:       PetscFree(ii[0]);
    156:       PetscFree(ii);
    157:       PetscFree(mcolors);
    158:     }
    159:     *isis = iscoloring->is;
    160:   }

    162:   return(0);
    163: }

    167: /*@C
    168:    ISColoringGetIS - Restores the index sets extracted from the coloring context

    170:    Collective on ISColoring

    172:    Input Parameter:
    173: +  iscoloring - the coloring context
    174: -  is - array of index sets

    176:    Level: advanced

    178: .seealso: ISColoringGetIS(), ISColoringView()
    179: @*/
    180: PetscErrorCode PETSCVEC_DLLEXPORT ISColoringRestoreIS(ISColoring iscoloring,IS *is[])
    181: {
    184:
    185:   /* currently nothing is done here */

    187:   return(0);
    188: }

    193: /*@C
    194:     ISColoringCreate - Generates an ISColoring context from lists (provided
    195:     by each processor) of colors for each node.

    197:     Collective on MPI_Comm

    199:     Input Parameters:
    200: +   comm - communicator for the processors creating the coloring
    201: .   n - number of nodes on this processor
    202: -   colors - array containing the colors for this processor, color
    203:              numbers begin at 0. In C/C++ this array must have been obtained with PetscMalloc()
    204:              and should NOT be freed (The ISColoringDestroy() will free it).

    206:     Output Parameter:
    207: .   iscoloring - the resulting coloring data structure

    209:     Options Database Key:
    210: .   -is_coloring_view - Activates ISColoringView()

    212:    Level: advanced
    213:
    214:     Notes: By default sets coloring type to  IS_COLORING_LOCAL

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

    218: @*/
    219: PetscErrorCode PETSCVEC_DLLEXPORT ISColoringCreate(MPI_Comm comm,PetscInt n,const ISColoringValue colors[],ISColoring *iscoloring)
    220: {
    222:   PetscMPIInt    size,rank,tag;
    223:   PetscInt       base,top,i;
    224:   PetscInt       nc,ncwork;
    225:   PetscTruth     flg;
    226:   MPI_Status     status;

    229:   PetscNew(struct _p_ISColoring,iscoloring);
    230:   PetscCommDuplicate(comm,&(*iscoloring)->comm,&tag);
    231:   comm = (*iscoloring)->comm;

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

    236:   /* should use MPI_Scan() */
    237:   MPI_Comm_rank(comm,&rank);
    238:   if (!rank) {
    239:     base = 0;
    240:     top  = n;
    241:   } else {
    242:     MPI_Recv(&base,1,MPIU_INT,rank-1,tag,comm,&status);
    243:     top = base+n;
    244:   }
    245:   if (rank < size-1) {
    246:     MPI_Send(&top,1,MPIU_INT,rank+1,tag,comm);
    247:   }

    249:   /* compute the total number of colors */
    250:   ncwork = 0;
    251:   for (i=0; i<n; i++) {
    252:     if (ncwork < colors[i]) ncwork = colors[i];
    253:   }
    254:   ncwork++;
    255:   MPI_Allreduce(&ncwork,&nc,1,MPIU_INT,MPI_MAX,comm);
    256:   (*iscoloring)->n      = nc;
    257:   (*iscoloring)->is     = 0;
    258:   (*iscoloring)->colors = (ISColoringValue *)colors;
    259:   (*iscoloring)->N      = n;
    260:   (*iscoloring)->refct  = 1;
    261:   (*iscoloring)->ctype  = IS_COLORING_LOCAL;

    263:   PetscOptionsHasName(PETSC_NULL,"-is_coloring_view",&flg);
    264:   if (flg) {
    265:     ISColoringView(*iscoloring,PETSC_VIEWER_STDOUT_((*iscoloring)->comm));
    266:   }
    267:   PetscLogInfo((0,"ISColoringCreate: Number of colors %d\n",nc));
    268:   return(0);
    269: }

    273: /*@C
    274:     ISPartitioningToNumbering - Takes an ISPartitioning and on each processor
    275:     generates an IS that contains a new global node number for each index based
    276:     on the partitioing.

    278:     Collective on IS

    280:     Input Parameters
    281: .   partitioning - a partitioning as generated by MatPartitioningApply()

    283:     Output Parameter:
    284: .   is - on each processor the index set that defines the global numbers
    285:          (in the new numbering) for all the nodes currently (before the partitioning)
    286:          on that processor

    288:    Level: advanced

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

    292: @*/
    293: PetscErrorCode PETSCVEC_DLLEXPORT ISPartitioningToNumbering(IS part,IS *is)
    294: {
    295:   MPI_Comm       comm;
    296:   PetscInt       i,*indices,np,npt,n,*starts,*sums,*lsizes,*newi;
    298:   PetscMPIInt    size;

    301:   PetscObjectGetComm((PetscObject)part,&comm);
    302:   MPI_Comm_size(comm,&size);

    304:   /* count the number of partitions, i.e., virtual processors */
    305:   ISGetLocalSize(part,&n);
    306:   ISGetIndices(part,&indices);
    307:   np = 0;
    308:   for (i=0; i<n; i++) {
    309:     np = PetscMax(np,indices[i]);
    310:   }
    311:   MPI_Allreduce(&np,&npt,1,MPIU_INT,MPI_MAX,comm);
    312:   np = npt+1; /* so that it looks like a MPI_Comm_size output */

    314:   /*
    315:         lsizes - number of elements of each partition on this particular processor
    316:         sums - total number of "previous" nodes for any particular partition
    317:         starts - global number of first element in each partition on this processor
    318:   */
    319:   PetscMalloc3(np,PetscInt,&lsizes,np,PetscInt,&starts,np,PetscInt,&sums);
    320:   PetscMemzero(lsizes,np*sizeof(PetscInt));
    321:   for (i=0; i<n; i++) {
    322:     lsizes[indices[i]]++;
    323:   }
    324:   MPI_Allreduce(lsizes,sums,np,MPIU_INT,MPI_SUM,comm);
    325:   MPI_Scan(lsizes,starts,np,MPIU_INT,MPI_SUM,comm);
    326:   for (i=0; i<np; i++) {
    327:     starts[i] -= lsizes[i];
    328:   }
    329:   for (i=1; i<np; i++) {
    330:     sums[i]    += sums[i-1];
    331:     starts[i]  += sums[i-1];
    332:   }

    334:   /*
    335:       For each local index give it the new global number
    336:   */
    337:   PetscMalloc(n*sizeof(PetscInt),&newi);
    338:   for (i=0; i<n; i++) {
    339:     newi[i] = starts[indices[i]]++;
    340:   }
    341:   PetscFree3(lsizes,starts,sums);

    343:   ISRestoreIndices(part,&indices);
    344:   ISCreateGeneral(comm,n,newi,is);
    345:   PetscFree(newi);
    346:   ISSetPermutation(*is);
    347:   return(0);
    348: }

    352: /*@C
    353:     ISPartitioningCount - Takes a ISPartitioning and determines the number of
    354:     resulting elements on each processor

    356:     Collective on IS

    358:     Input Parameters:
    359: .   partitioning - a partitioning as generated by MatPartitioningApply()

    361:     Output Parameter:
    362: .   count - array of length size, to contain the number of elements assigned
    363:         to each partition, where size is the number of partitions generated
    364:          (see notes below).

    366:    Level: advanced

    368:     Notes:
    369:         By default the number of partitions generated (and thus the length
    370:         of count) is the size of the communicator associated with IS,
    371:         but it can be set by MatPartitioningSetNParts. The resulting array
    372:         of lengths can for instance serve as input of PCBJacobiSetTotalBlocks.

    375: .seealso: MatPartitioningCreate(), AOCreateBasic(), ISPartitioningToNumbering(),
    376:         MatPartitioningSetNParts()

    378: @*/
    379: PetscErrorCode PETSCVEC_DLLEXPORT ISPartitioningCount(IS part,PetscInt count[])
    380: {
    381:   MPI_Comm       comm;
    382:   PetscInt            i,*indices,np,npt,n,*lsizes;
    384:   PetscMPIInt    size;

    387:   PetscObjectGetComm((PetscObject)part,&comm);
    388:   MPI_Comm_size(comm,&size);

    390:   /* count the number of partitions */
    391:   ISGetLocalSize(part,&n);
    392:   ISGetIndices(part,&indices);
    393:   np = 0;
    394:   for (i=0; i<n; i++) {
    395:     np = PetscMax(np,indices[i]);
    396:   }
    397:   MPI_Allreduce(&np,&npt,1,MPIU_INT,MPI_MAX,comm);
    398:   np = npt+1; /* so that it looks like a MPI_Comm_size output */

    400:   /*
    401:         lsizes - number of elements of each partition on this particular processor
    402:         sums - total number of "previous" nodes for any particular partition
    403:         starts - global number of first element in each partition on this processor
    404:   */
    405:   PetscMalloc(np*sizeof(PetscInt),&lsizes);
    406:   PetscMemzero(lsizes,np*sizeof(PetscInt));
    407:   for (i=0; i<n; i++) {
    408:     lsizes[indices[i]]++;
    409:   }
    410:   ISRestoreIndices(part,&indices);
    411:   MPI_Allreduce(lsizes,count,np,MPIU_INT,MPI_SUM,comm);
    412:   PetscFree(lsizes);
    413:   return(0);
    414: }

    418: /*@C
    419:     ISAllGather - Given an index set (IS) on each processor, generates a large
    420:     index set (same on each processor) by concatenating together each
    421:     processors index set.

    423:     Collective on IS

    425:     Input Parameter:
    426: .   is - the distributed index set

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

    431:     Notes:
    432:     ISAllGather() is clearly not scalable for large index sets.

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

    440:     Level: intermediate

    442:     Concepts: gather^index sets
    443:     Concepts: index sets^gathering to all processors
    444:     Concepts: IS^gathering to all processors

    446: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGatherIndices()
    447: @*/
    448: PetscErrorCode PETSCVEC_DLLEXPORT ISAllGather(IS is,IS *isout)
    449: {
    451:   PetscInt       *indices,n,*lindices,i,N;
    452:   MPI_Comm       comm;
    453:   PetscMPIInt    size,*sizes,*offsets,nn;

    459:   PetscObjectGetComm((PetscObject)is,&comm);
    460:   MPI_Comm_size(comm,&size);
    461:   PetscMalloc2(size,PetscMPIInt,&sizes,size,PetscMPIInt,&offsets);
    462:
    463:   ISGetLocalSize(is,&n);
    464:   nn   = (PetscMPIInt)n;
    465:   MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
    466:   offsets[0] = 0;
    467:   for (i=1;i<size; i++) offsets[i] = offsets[i-1] + sizes[i-1];
    468:   N = offsets[size-1] + sizes[size-1];

    470:   PetscMalloc(N*sizeof(PetscInt),&indices);
    471:   ISGetIndices(is,&lindices);
    472:   MPI_Allgatherv(lindices,nn,MPIU_INT,indices,sizes,offsets,MPIU_INT,comm);
    473:   ISRestoreIndices(is,&lindices);
    474:   PetscFree(sizes);

    476:   ISCreateGeneral(PETSC_COMM_SELF,N,indices,isout);
    477:   PetscFree2(indices,offsets);
    478:   return(0);
    479: }

    483: /*@C
    484:     ISAllGatherIndices - Given a a set of integers on each processor, generates a large
    485:     set (same on each processor) by concatenating together each processors integers

    487:     Collective on MPI_Comm

    489:     Input Parameter:
    490: +   comm - communicator to share the indices
    491: .   n - local size of set
    492: -   lindices - local indices

    494:     Output Parameter:
    495: +   outN - total number of indices
    496: -   outindices - all of the integers

    498:     Notes:
    499:     ISAllGatherIndices() is clearly not scalable for large index sets.

    502:     Level: intermediate

    504:     Concepts: gather^index sets
    505:     Concepts: index sets^gathering to all processors
    506:     Concepts: IS^gathering to all processors

    508: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
    509: @*/
    510: PetscErrorCode PETSCVEC_DLLEXPORT ISAllGatherIndices(MPI_Comm comm,PetscInt n,const PetscInt lindices[],PetscInt *outN,PetscInt *outindices[])
    511: {
    513:   PetscInt       *indices,i,N;
    514:   PetscMPIInt    size,*sizes,*offsets,nn;

    517:   MPI_Comm_size(comm,&size);
    518:   PetscMalloc2(size,PetscMPIInt,&sizes,size,PetscMPIInt,&offsets);
    519:
    520:   nn   = n;
    521:   MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
    522:   offsets[0] = 0;
    523:   for (i=1;i<size; i++) offsets[i] = offsets[i-1] + sizes[i-1];
    524:   N    = offsets[size-1] + sizes[size-1];

    526:   PetscMalloc(N*sizeof(PetscInt),&indices);
    527:   MPI_Allgatherv((void*)lindices,nn,MPIU_INT,indices,sizes,offsets,MPIU_INT,comm);
    528:   PetscFree2(sizes,offsets);

    530:   *outindices = indices;
    531:   if (outN) *outN = N;
    532:   return(0);
    533: }

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

    543:     Collective on MPI_Comm

    545:     Input Parameter:
    546: +   comm - communicator to share the indices
    547: .   n - local size of set
    548: -   lindices - local colors

    550:     Output Parameter:
    551: +   outN - total number of indices
    552: -   outindices - all of the colors

    554:     Notes:
    555:     ISAllGatherColors() is clearly not scalable for large index sets.

    558:     Level: intermediate

    560:     Concepts: gather^index sets
    561:     Concepts: index sets^gathering to all processors
    562:     Concepts: IS^gathering to all processors

    564: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather(), ISAllGatherIndices()
    565: @*/
    566: PetscErrorCode PETSCVEC_DLLEXPORT ISAllGatherColors(MPI_Comm comm,PetscInt n,ISColoringValue *lindices,PetscInt *outN,ISColoringValue *outindices[])
    567: {
    568:   ISColoringValue *indices;
    569:   PetscErrorCode  ierr;
    570:   PetscInt        i,N;
    571:   PetscMPIInt     size,*offsets,*sizes, nn = n;

    574:   MPI_Comm_size(comm,&size);
    575:   PetscMalloc2(size,PetscMPIInt,&sizes,size,PetscMPIInt,&offsets);
    576:
    577:   MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
    578:   offsets[0] = 0;
    579:   for (i=1;i<size; i++) offsets[i] = offsets[i-1] + sizes[i-1];
    580:   N    = offsets[size-1] + sizes[size-1];
    581:   PetscFree2(sizes,offsets);

    583:   PetscMalloc((N+1)*sizeof(ISColoringValue),&indices);
    584:   MPI_Allgatherv(lindices,(PetscMPIInt)n,MPIU_COLORING_VALUE,indices,sizes,offsets,MPIU_COLORING_VALUE,comm);

    586:   *outindices = indices;
    587:   if (outN) *outN = N;
    588:   return(0);
    589: }