Actual source code: block.c

      1: #define PETSCVEC_DLL
      2: /*
      3:      Provides the functions for index sets (IS) defined by a list of integers.
      4:    These are for blocks of data, each block is indicated with a single integer.
      5: */
     6:  #include src/vec/is/isimpl.h
     7:  #include petscsys.h

      9: EXTERN PetscErrorCode PETSCVEC_DLLEXPORT VecInitializePackage(char *);

     11: typedef struct {
     12:   PetscInt        N,n;            /* number of blocks */
     13:   PetscTruth      sorted;       /* are the blocks sorted? */
     14:   PetscInt        *idx;
     15:   PetscInt        bs;           /* blocksize */
     16: } IS_Block;

     20: PetscErrorCode ISDestroy_Block(IS is)
     21: {
     22:   IS_Block       *is_block = (IS_Block*)is->data;

     26:   PetscFree(is_block->idx);
     27:   PetscFree(is_block);
     28:   PetscHeaderDestroy(is);
     29:   return(0);
     30: }

     34: PetscErrorCode ISGetIndices_Block(IS in,PetscInt **idx)
     35: {
     36:   IS_Block       *sub = (IS_Block*)in->data;
     38:   PetscInt       i,j,k,bs = sub->bs,n = sub->n,*ii,*jj;

     41:   if (sub->bs == 1) {
     42:     *idx = sub->idx;
     43:   } else {
     44:     PetscMalloc(sub->bs*sub->n*sizeof(PetscInt),&jj);
     45:     *idx = jj;
     46:     k    = 0;
     47:     ii   = sub->idx;
     48:     for (i=0; i<n; i++) {
     49:       for (j=0; j<bs; j++) {
     50:         jj[k++] = ii[i] + j;
     51:       }
     52:     }
     53:   }
     54:   return(0);
     55: }

     59: PetscErrorCode ISRestoreIndices_Block(IS in,PetscInt **idx)
     60: {
     61:   IS_Block       *sub = (IS_Block*)in->data;

     65:   if (sub->bs != 1) {
     66:     PetscFree(*idx);
     67:   } else {
     68:     if (*idx !=  sub->idx) {
     69:       SETERRQ(PETSC_ERR_ARG_WRONG,"Must restore with value from ISGetIndices()");
     70:     }
     71:   }
     72:   return(0);
     73: }

     77: PetscErrorCode ISGetSize_Block(IS is,PetscInt *size)
     78: {
     79:   IS_Block *sub = (IS_Block *)is->data;

     82:   *size = sub->bs*sub->N;
     83:   return(0);
     84: }

     88: PetscErrorCode ISGetLocalSize_Block(IS is,PetscInt *size)
     89: {
     90:   IS_Block *sub = (IS_Block *)is->data;

     93:   *size = sub->bs*sub->n;
     94:   return(0);
     95: }

     99: PetscErrorCode ISInvertPermutation_Block(IS is,PetscInt nlocal,IS *isout)
    100: {
    101:   IS_Block       *sub = (IS_Block *)is->data;
    102:   PetscInt       i,*ii,n = sub->n,*idx = sub->idx;
    103:   PetscMPIInt    size;

    107:   MPI_Comm_size(is->comm,&size);
    108:   if (size == 1) {
    109:     PetscMalloc(n*sizeof(PetscInt),&ii);
    110:     for (i=0; i<n; i++) {
    111:       ii[idx[i]] = i;
    112:     }
    113:     ISCreateBlock(PETSC_COMM_SELF,sub->bs,n,ii,isout);
    114:     ISSetPermutation(*isout);
    115:     PetscFree(ii);
    116:   } else {
    117:     SETERRQ(PETSC_ERR_SUP,"No inversion written yet for block IS");
    118:   }
    119:   return(0);
    120: }

    124: PetscErrorCode ISView_Block(IS is, PetscViewer viewer)
    125: {
    126:   IS_Block       *sub = (IS_Block *)is->data;
    128:   PetscInt       i,n = sub->n,*idx = sub->idx;
    129:   PetscTruth     iascii;

    132:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    133:   if (iascii) {
    134:     if (is->isperm) {
    135:       PetscViewerASCIISynchronizedPrintf(viewer,"Block Index set is permutation\n");
    136:     }
    137:     PetscViewerASCIISynchronizedPrintf(viewer,"Block size %D\n",sub->bs);
    138:     PetscViewerASCIISynchronizedPrintf(viewer,"Number of block indices in set %D\n",n);
    139:     PetscViewerASCIISynchronizedPrintf(viewer,"The first indices of each block are\n");
    140:     for (i=0; i<n; i++) {
    141:       PetscViewerASCIISynchronizedPrintf(viewer,"%D %D\n",i,idx[i]);
    142:     }
    143:     PetscViewerFlush(viewer);
    144:   } else {
    145:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for this object",((PetscObject)viewer)->type_name);
    146:   }
    147:   return(0);
    148: }

    152: PetscErrorCode ISSort_Block(IS is)
    153: {
    154:   IS_Block       *sub = (IS_Block *)is->data;

    158:   if (sub->sorted) return(0);
    159:   PetscSortInt(sub->n,sub->idx);
    160:   sub->sorted = PETSC_TRUE;
    161:   return(0);
    162: }

    166: PetscErrorCode ISSorted_Block(IS is,PetscTruth *flg)
    167: {
    168:   IS_Block *sub = (IS_Block *)is->data;

    171:   *flg = sub->sorted;
    172:   return(0);
    173: }

    177: PetscErrorCode ISDuplicate_Block(IS is,IS *newIS)
    178: {
    180:   IS_Block       *sub = (IS_Block *)is->data;

    183:   ISCreateBlock(is->comm,sub->bs,sub->n,sub->idx,newIS);
    184:   return(0);
    185: }

    189: PetscErrorCode ISIdentity_Block(IS is,PetscTruth *ident)
    190: {
    191:   IS_Block *is_block = (IS_Block*)is->data;
    192:   PetscInt i,n = is_block->n,*idx = is_block->idx,bs = is_block->bs;

    195:   is->isidentity = PETSC_TRUE;
    196:   *ident         = PETSC_TRUE;
    197:   for (i=0; i<n; i++) {
    198:     if (idx[i] != bs*i) {
    199:       is->isidentity = PETSC_FALSE;
    200:       *ident         = PETSC_FALSE;
    201:       return(0);
    202:     }
    203:   }
    204:   return(0);
    205: }

    207: static struct _ISOps myops = { ISGetSize_Block,
    208:                                ISGetLocalSize_Block,
    209:                                ISGetIndices_Block,
    210:                                ISRestoreIndices_Block,
    211:                                ISInvertPermutation_Block,
    212:                                ISSort_Block,
    213:                                ISSorted_Block,
    214:                                ISDuplicate_Block,
    215:                                ISDestroy_Block,
    216:                                ISView_Block,
    217:                                ISIdentity_Block };
    220: /*@C
    221:    ISCreateBlock - Creates a data structure for an index set containing
    222:    a list of integers. The indices are relative to entries, not blocks.

    224:    Collective on MPI_Comm

    226:    Input Parameters:
    227: +  n - the length of the index set (the number of blocks)
    228: .  bs - number of elements in each block
    229: .  idx - the list of integers
    230: -  comm - the MPI communicator

    232:    Output Parameter:
    233: .  is - the new index set

    235:    Notes:
    236:    When the communicator is not MPI_COMM_SELF, the operations on the
    237:    index sets, IS, are NOT conceptually the same as MPI_Group operations.
    238:    The index sets are then distributed sets of indices and thus certain operations
    239:    on them are collective.

    241:    Example:
    242:    If you wish to index the values {0,1,4,5}, then use
    243:    a block size of 2 and idx of {0,4}.

    245:    Level: beginner

    247:   Concepts: IS^block
    248:   Concepts: index sets^block
    249:   Concepts: block^index set

    251: .seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather()
    252: @*/
    253: PetscErrorCode PETSCVEC_DLLEXPORT ISCreateBlock(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscInt idx[],IS *is)
    254: {
    256:   PetscInt       i,min,max;
    257:   IS             Nindex;
    258:   IS_Block       *sub;
    259:   PetscTruth     sorted = PETSC_TRUE;

    263:   if (n < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
    265:   *is = PETSC_NULL;
    266: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
    267:   VecInitializePackage(PETSC_NULL);
    268: #endif

    270:   PetscHeaderCreate(Nindex,_p_IS,struct _ISOps,IS_COOKIE,IS_BLOCK,"IS",comm,ISDestroy,ISView);
    271:   PetscNew(IS_Block,&sub);
    272:   PetscLogObjectMemory(Nindex,sizeof(IS_Block)+n*sizeof(PetscInt)+sizeof(struct _p_IS));
    273:   PetscMalloc(n*sizeof(PetscInt),&sub->idx);
    274:   sub->n = n;
    275:   MPI_Allreduce(&n,&sub->N,1,MPIU_INT,MPI_SUM,comm);
    276:   for (i=1; i<n; i++) {
    277:     if (idx[i] < idx[i-1]) {sorted = PETSC_FALSE; break;}
    278:   }
    279:   if (n) {min = max = idx[0];} else {min = max = 0;}
    280:   for (i=1; i<n; i++) {
    281:     if (idx[i] < min) min = idx[i];
    282:     if (idx[i] > max) max = idx[i];
    283:   }
    284:   PetscMemcpy(sub->idx,idx,n*sizeof(PetscInt));
    285:   sub->sorted     = sorted;
    286:   sub->bs         = bs;
    287:   Nindex->min     = min;
    288:   Nindex->max     = max;
    289:   Nindex->data    = (void*)sub;
    290:   PetscMemcpy(Nindex->ops,&myops,sizeof(myops));
    291:   Nindex->isperm  = PETSC_FALSE;
    292:   *is = Nindex; return(0);
    293: }

    298: /*@C
    299:    ISBlockGetIndices - Gets the indices associated with each block.

    301:    Not Collective

    303:    Input Parameter:
    304: .  is - the index set

    306:    Output Parameter:
    307: .  idx - the integer indices

    309:    Level: intermediate

    311:    Concepts: IS^block
    312:    Concepts: index sets^getting indices
    313:    Concepts: index sets^block

    315: .seealso: ISGetIndices(), ISBlockRestoreIndices()
    316: @*/
    317: PetscErrorCode PETSCVEC_DLLEXPORT ISBlockGetIndices(IS in,PetscInt *idx[])
    318: {
    319:   IS_Block *sub;

    324:   if (in->type != IS_BLOCK) SETERRQ(PETSC_ERR_ARG_WRONG,"Not a block index set");

    326:   sub = (IS_Block*)in->data;
    327:   *idx = sub->idx;
    328:   return(0);
    329: }

    333: /*@C
    334:    ISBlockRestoreIndices - Restores the indices associated with each block.

    336:    Not Collective

    338:    Input Parameter:
    339: .  is - the index set

    341:    Output Parameter:
    342: .  idx - the integer indices

    344:    Level: intermediate

    346:    Concepts: IS^block
    347:    Concepts: index sets^getting indices
    348:    Concepts: index sets^block

    350: .seealso: ISRestoreIndices(), ISBlockGetIndices()
    351: @*/
    352: PetscErrorCode PETSCVEC_DLLEXPORT ISBlockRestoreIndices(IS is,PetscInt *idx[])
    353: {
    357:   if (is->type != IS_BLOCK) SETERRQ(PETSC_ERR_ARG_WRONG,"Not a block index set");
    358:   return(0);
    359: }

    363: /*@
    364:    ISBlockGetBlockSize - Returns the number of elements in a block.

    366:    Not Collective

    368:    Input Parameter:
    369: .  is - the index set

    371:    Output Parameter:
    372: .  size - the number of elements in a block

    374:    Level: intermediate

    376:    Concepts: IS^block size
    377:    Concepts: index sets^block size

    379: .seealso: ISBlockGetSize(), ISGetSize(), ISBlock(), ISCreateBlock()
    380: @*/
    381: PetscErrorCode PETSCVEC_DLLEXPORT ISBlockGetBlockSize(IS is,PetscInt *size)
    382: {
    383:   IS_Block *sub;

    388:   if (is->type != IS_BLOCK) SETERRQ(PETSC_ERR_ARG_WRONG,"Not a block index set");

    390:   sub = (IS_Block *)is->data;
    391:   *size = sub->bs;
    392:   return(0);
    393: }

    397: /*@C
    398:    ISBlock - Checks whether an index set is blocked.

    400:    Not Collective

    402:    Input Parameter:
    403: .  is - the index set

    405:    Output Parameter:
    406: .  flag - PETSC_TRUE if a block index set, else PETSC_FALSE

    408:    Level: intermediate

    410:    Concepts: IS^block
    411:    Concepts: index sets^block

    413: .seealso: ISBlockGetSize(), ISGetSize(), ISBlockGetBlockSize(), ISCreateBlock()
    414: @*/
    415: PetscErrorCode PETSCVEC_DLLEXPORT ISBlock(IS is,PetscTruth *flag)
    416: {
    420:   if (is->type != IS_BLOCK) *flag = PETSC_FALSE;
    421:   else                          *flag = PETSC_TRUE;
    422:   return(0);
    423: }

    427: /*@
    428:    ISBlockGetSize - Returns the number of blocks in the index set.

    430:    Not Collective

    432:    Input Parameter:
    433: .  is - the index set

    435:    Output Parameter:
    436: .  size - the number of blocks

    438:    Level: intermediate

    440:    Concepts: IS^block sizes
    441:    Concepts: index sets^block sizes

    443: .seealso: ISBlockGetBlockSize(), ISGetSize(), ISBlock(), ISCreateBlock()
    444: @*/
    445: PetscErrorCode PETSCVEC_DLLEXPORT ISBlockGetSize(IS is,PetscInt *size)
    446: {
    447:   IS_Block *sub;

    452:   if (is->type != IS_BLOCK) SETERRQ(PETSC_ERR_ARG_WRONG,"Not a block index set");

    454:   sub = (IS_Block *)is->data;
    455:   *size = sub->n;
    456:   return(0);
    457: }