impls/block/block.c

/*
     Provides the functions for index sets (IS) defined by a list of integers.
   These are for blocks of data, each block is indicated with a single integer.
*/
#include <petsc/private/isimpl.h>               /*I  "petscis.h"     I*/
#include <petscviewer.h>

typedef struct {
  PetscBool sorted;      /* are the blocks sorted? */
  PetscBool allocated;   /* did we allocate the index array ourselves? */
  PetscInt  *idx;
} IS_Block;

static PetscErrorCode ISDestroy_Block(IS is)
{
  IS_Block       *sub = (IS_Block*)is->data;

  PetscFunctionBegin;
  if (sub->allocated) PetscCall(PetscFree(sub->idx));
  PetscCall(PetscObjectComposeFunction((PetscObject)is,"ISBlockSetIndices_C",NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)is,"ISBlockGetIndices_C",NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)is,"ISBlockRestoreIndices_C",NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)is,"ISBlockGetSize_C",NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)is,"ISBlockGetLocalSize_C",NULL));
  PetscCall(PetscFree(is->data));
  PetscFunctionReturn(0);
}

static PetscErrorCode ISLocate_Block(IS is,PetscInt key,PetscInt *location)
{
  IS_Block       *sub = (IS_Block*)is->data;
  PetscInt       numIdx, i, bs, bkey, mkey;
  PetscBool      sorted;

  PetscFunctionBegin;
  PetscCall(PetscLayoutGetBlockSize(is->map,&bs));
  PetscCall(PetscLayoutGetSize(is->map,&numIdx));
  numIdx /= bs;
  bkey    = key / bs;
  mkey    = key % bs;
  if (mkey < 0) {
    bkey--;
    mkey += bs;
  }
  PetscCall(ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,&sorted));
  if (sorted) {
    PetscCall(PetscFindInt(bkey,numIdx,sub->idx,location));
  } else {
    const PetscInt *idx = sub->idx;

    *location = -1;
    for (i = 0; i < numIdx; i++) {
      if (idx[i] == bkey) {
        *location = i;
        break;
      }
    }
  }
  if (*location >= 0) {
    *location = *location * bs + mkey;
  }
  PetscFunctionReturn(0);
}

static PetscErrorCode ISGetIndices_Block(IS in,const PetscInt *idx[])
{
  IS_Block       *sub = (IS_Block*)in->data;
  PetscInt       i,j,k,bs,n,*ii,*jj;

  PetscFunctionBegin;
  PetscCall(PetscLayoutGetBlockSize(in->map, &bs));
  PetscCall(PetscLayoutGetLocalSize(in->map, &n));
  n   /= bs;
  if (bs == 1) *idx = sub->idx;
  else {
    if (n) {
      PetscCall(PetscMalloc1(bs*n,&jj));
      *idx = jj;
      k    = 0;
      ii   = sub->idx;
      for (i=0; i<n; i++)
        for (j=0; j<bs; j++)
          jj[k++] = bs*ii[i] + j;
    } else {
      /* do not malloc for zero size because F90Array1dCreate() inside ISRestoreArrayF90() does not keep array when zero length array */
      *idx = NULL;
    }
  }
  PetscFunctionReturn(0);
}

static PetscErrorCode ISRestoreIndices_Block(IS is,const PetscInt *idx[])
{
  IS_Block       *sub = (IS_Block*)is->data;
  PetscInt       bs;

  PetscFunctionBegin;
  PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
  if (bs != 1) {
    PetscCall(PetscFree(*(void**)idx));
  } else {
    /* F90Array1dCreate() inside ISRestoreArrayF90() does not keep array when zero length array */
    PetscCheck(is->map->n <= 0 || *idx == sub->idx,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must restore with value from ISGetIndices()");
  }
  PetscFunctionReturn(0);
}

static PetscErrorCode ISInvertPermutation_Block(IS is,PetscInt nlocal,IS *isout)
{
  IS_Block       *sub = (IS_Block*)is->data;
  PetscInt       i,*ii,bs,n,*idx = sub->idx;
  PetscMPIInt    size;

  PetscFunctionBegin;
  PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)is),&size));
  PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
  PetscCall(PetscLayoutGetLocalSize(is->map, &n));
  n   /= bs;
  if (size == 1) {
    PetscCall(PetscMalloc1(n,&ii));
    for (i=0; i<n; i++) ii[idx[i]] = i;
    PetscCall(ISCreateBlock(PETSC_COMM_SELF,bs,n,ii,PETSC_OWN_POINTER,isout));
    PetscCall(ISSetPermutation(*isout));
  } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No inversion written yet for block IS");
  PetscFunctionReturn(0);
}

static PetscErrorCode ISView_Block(IS is, PetscViewer viewer)
{
  IS_Block       *sub = (IS_Block*)is->data;
  PetscInt       i,bs,n,*idx = sub->idx;
  PetscBool      iascii,ibinary;

  PetscFunctionBegin;
  PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
  PetscCall(PetscLayoutGetLocalSize(is->map, &n));
  n   /= bs;
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii));
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&ibinary));
  if (iascii) {
    PetscViewerFormat fmt;

    PetscCall(PetscViewerGetFormat(viewer,&fmt));
    if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
      IS             ist;
      const char     *name;
      const PetscInt *idx;
      PetscInt       n;

      PetscCall(PetscObjectGetName((PetscObject)is,&name));
      PetscCall(ISGetLocalSize(is,&n));
      PetscCall(ISGetIndices(is,&idx));
      PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is),n,idx,PETSC_USE_POINTER,&ist));
      PetscCall(PetscObjectSetName((PetscObject)ist,name));
      PetscCall(ISView(ist,viewer));
      PetscCall(ISDestroy(&ist));
      PetscCall(ISRestoreIndices(is,&idx));
    } else {
      PetscBool isperm;

      PetscCall(ISGetInfo(is,IS_PERMUTATION,IS_GLOBAL,PETSC_FALSE,&isperm));
      if (isperm) PetscCall(PetscViewerASCIIPrintf(viewer,"Block Index set is permutation\n"));
      PetscCall(PetscViewerASCIIPushSynchronized(viewer));
      PetscCall(PetscViewerASCIISynchronizedPrintf(viewer,"Block size %" PetscInt_FMT "\n",bs));
      PetscCall(PetscViewerASCIISynchronizedPrintf(viewer,"Number of block indices in set %" PetscInt_FMT "\n",n));
      PetscCall(PetscViewerASCIISynchronizedPrintf(viewer,"The first indices of each block are\n"));
      for (i=0; i<n; i++) {
        PetscCall(PetscViewerASCIISynchronizedPrintf(viewer,"Block %" PetscInt_FMT " Index %" PetscInt_FMT "\n",i,idx[i]));
      }
      PetscCall(PetscViewerFlush(viewer));
      PetscCall(PetscViewerASCIIPopSynchronized(viewer));
    }
  } else if (ibinary) {
    PetscCall(ISView_Binary(is,viewer));
  }
  PetscFunctionReturn(0);
}

static PetscErrorCode ISSort_Block(IS is)
{
  IS_Block       *sub = (IS_Block*)is->data;
  PetscInt       bs, n;

  PetscFunctionBegin;
  PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
  PetscCall(PetscLayoutGetLocalSize(is->map, &n));
  PetscCall(PetscIntSortSemiOrdered(n/bs,sub->idx));
  PetscFunctionReturn(0);
}

static PetscErrorCode ISSortRemoveDups_Block(IS is)
{
  IS_Block       *sub = (IS_Block*)is->data;
  PetscInt       bs, n, nb;
  PetscBool      sorted;

  PetscFunctionBegin;
  PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
  PetscCall(PetscLayoutGetLocalSize(is->map, &n));
  nb   = n/bs;
  PetscCall(ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,&sorted));
  if (sorted) {
    PetscCall(PetscSortedRemoveDupsInt(&nb,sub->idx));
  } else {
    PetscCall(PetscSortRemoveDupsInt(&nb,sub->idx));
  }
  PetscCall(PetscLayoutDestroy(&is->map));
  PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), nb*bs, PETSC_DECIDE, bs, &is->map));
  PetscFunctionReturn(0);
}

static PetscErrorCode ISSorted_Block(IS is,PetscBool  *flg)
{
  PetscFunctionBegin;
  PetscCall(ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,flg));
  PetscFunctionReturn(0);
}

static PetscErrorCode ISSortedLocal_Block(IS is,PetscBool *flg)
{
  IS_Block       *sub = (IS_Block*)is->data;
  PetscInt       n, bs, i, *idx;

  PetscFunctionBegin;
  PetscCall(PetscLayoutGetLocalSize(is->map, &n));
  PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
  n   /= bs;
  idx  = sub->idx;
  for (i = 1; i < n; i++) if (idx[i] < idx[i - 1]) break;
  if (i < n) *flg = PETSC_FALSE;
  else       *flg = PETSC_TRUE;
  PetscFunctionReturn(0);
}

static PetscErrorCode ISUniqueLocal_Block(IS is,PetscBool *flg)
{
  IS_Block       *sub = (IS_Block*)is->data;
  PetscInt       n, bs, i, *idx, *idxcopy = NULL;
  PetscBool      sortedLocal;

  PetscFunctionBegin;
  PetscCall(PetscLayoutGetLocalSize(is->map, &n));
  PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
  n   /= bs;
  idx  = sub->idx;
  PetscCall(ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,&sortedLocal));
  if (!sortedLocal) {
    PetscCall(PetscMalloc1(n, &idxcopy));
    PetscCall(PetscArraycpy(idxcopy, idx, n));
    PetscCall(PetscIntSortSemiOrdered(n, idxcopy));
    idx = idxcopy;
  }
  for (i = 1; i < n; i++) if (idx[i] == idx[i - 1]) break;
  if (i < n) *flg = PETSC_FALSE;
  else       *flg = PETSC_TRUE;
  PetscCall(PetscFree(idxcopy));
  PetscFunctionReturn(0);
}

static PetscErrorCode ISPermutationLocal_Block(IS is,PetscBool *flg)
{
  IS_Block       *sub = (IS_Block*)is->data;
  PetscInt       n, bs, i, *idx, *idxcopy = NULL;
  PetscBool      sortedLocal;

  PetscFunctionBegin;
  PetscCall(PetscLayoutGetLocalSize(is->map, &n));
  PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
  n   /= bs;
  idx  = sub->idx;
  PetscCall(ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,&sortedLocal));
  if (!sortedLocal) {
    PetscCall(PetscMalloc1(n, &idxcopy));
    PetscCall(PetscArraycpy(idxcopy, idx, n));
    PetscCall(PetscIntSortSemiOrdered(n, idxcopy));
    idx = idxcopy;
  }
  for (i = 0; i < n; i++) if (idx[i] != i) break;
  if (i < n) *flg = PETSC_FALSE;
  else       *flg = PETSC_TRUE;
  PetscCall(PetscFree(idxcopy));
  PetscFunctionReturn(0);
}

static PetscErrorCode ISIntervalLocal_Block(IS is,PetscBool *flg)
{
  IS_Block       *sub = (IS_Block*)is->data;
  PetscInt       n, bs, i, *idx;

  PetscFunctionBegin;
  PetscCall(PetscLayoutGetLocalSize(is->map, &n));
  PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
  n   /= bs;
  idx  = sub->idx;
  for (i = 1; i < n; i++) if (idx[i] != idx[i - 1] + 1) break;
  if (i < n) *flg = PETSC_FALSE;
  else       *flg = PETSC_TRUE;
  PetscFunctionReturn(0);
}

static PetscErrorCode ISDuplicate_Block(IS is,IS *newIS)
{
  IS_Block       *sub = (IS_Block*)is->data;
  PetscInt        bs, n;

  PetscFunctionBegin;
  PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
  PetscCall(PetscLayoutGetLocalSize(is->map, &n));
  n   /= bs;
  PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)is),bs,n,sub->idx,PETSC_COPY_VALUES,newIS));
  PetscFunctionReturn(0);
}

static PetscErrorCode ISCopy_Block(IS is,IS isy)
{
  IS_Block       *is_block = (IS_Block*)is->data,*isy_block = (IS_Block*)isy->data;
  PetscInt       bs, n;

  PetscFunctionBegin;
  PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
  PetscCall(PetscLayoutGetLocalSize(is->map, &n));
  PetscCall(PetscArraycpy(isy_block->idx, is_block->idx, n/bs));
  PetscFunctionReturn(0);
}

static PetscErrorCode ISOnComm_Block(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
{
  IS_Block       *sub = (IS_Block*)is->data;
  PetscInt       bs, n;

  PetscFunctionBegin;
  PetscCheck(mode != PETSC_OWN_POINTER,comm,PETSC_ERR_ARG_WRONG,"Cannot use PETSC_OWN_POINTER");
  PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
  PetscCall(PetscLayoutGetLocalSize(is->map, &n));
  PetscCall(ISCreateBlock(comm,bs,n/bs,sub->idx,mode,newis));
  PetscFunctionReturn(0);
}

static PetscErrorCode ISSetBlockSize_Block(IS is,PetscInt bs)
{
  PetscFunctionBegin;
  PetscCheck(is->map->bs <= 0 || bs == is->map->bs,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot change blocksize %" PetscInt_FMT " (to %" PetscInt_FMT ") if ISType is ISBLOCK",is->map->bs,bs);
  PetscCall(PetscLayoutSetBlockSize(is->map, bs));
  PetscFunctionReturn(0);
}

static PetscErrorCode ISToGeneral_Block(IS inis)
{
  IS_Block       *sub   = (IS_Block*)inis->data;
  PetscInt       bs,n;
  const PetscInt *idx;

  PetscFunctionBegin;
  PetscCall(ISGetBlockSize(inis,&bs));
  PetscCall(ISGetLocalSize(inis,&n));
  PetscCall(ISGetIndices(inis,&idx));
  if (bs == 1) {
    PetscCopyMode mode = sub->allocated ? PETSC_OWN_POINTER : PETSC_USE_POINTER;
    sub->allocated = PETSC_FALSE; /* prevent deallocation when changing the subtype*/
    PetscCall(ISSetType(inis,ISGENERAL));
    PetscCall(ISGeneralSetIndices(inis,n,idx,mode));
  } else {
    PetscCall(ISSetType(inis,ISGENERAL));
    PetscCall(ISGeneralSetIndices(inis,n,idx,PETSC_OWN_POINTER));
  }
  PetscFunctionReturn(0);
}

static struct _ISOps myops = { ISGetIndices_Block,
                               ISRestoreIndices_Block,
                               ISInvertPermutation_Block,
                               ISSort_Block,
                               ISSortRemoveDups_Block,
                               ISSorted_Block,
                               ISDuplicate_Block,
                               ISDestroy_Block,
                               ISView_Block,
                               ISLoad_Default,
                               ISCopy_Block,
                               ISToGeneral_Block,
                               ISOnComm_Block,
                               ISSetBlockSize_Block,
                               NULL,
                               ISLocate_Block,
                               /* we can have specialized local routines for determining properties,
                                * but unless the block size is the same on each process (which is not guaranteed at
                                * the moment), then trying to do something specialized for global properties is too
                                * complicated */
                               ISSortedLocal_Block,
                               NULL,
                               ISUniqueLocal_Block,
                               NULL,
                               ISPermutationLocal_Block,
                               NULL,
                               ISIntervalLocal_Block,
                               NULL};

/*@
   ISBlockSetIndices - Set integers representing blocks of indices in an index set.

   Collective on IS

   Input Parameters:
+  is - the index set
.  bs - number of elements in each block
.   n - the length of the index set (the number of blocks)
.  idx - the list of integers, one for each block, the integers contain the index of the first index of each block divided by the block size
-  mode - see PetscCopyMode, only PETSC_COPY_VALUES and PETSC_OWN_POINTER are supported

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

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

   Level: beginner

.seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather(), ISCreateBlock(), ISBLOCK, ISGeneralSetIndices()
@*/
PetscErrorCode  ISBlockSetIndices(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
{
  PetscFunctionBegin;
  PetscCall(ISClearInfoCache(is,PETSC_FALSE));
  PetscUseMethod(is,"ISBlockSetIndices_C",(IS,PetscInt,PetscInt,const PetscInt[],PetscCopyMode),(is,bs,n,idx,mode));
  PetscFunctionReturn(0);
}

static PetscErrorCode  ISBlockSetIndices_Block(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
{
  PetscInt       i,min,max;
  IS_Block       *sub = (IS_Block*)is->data;
  PetscLayout    map;

  PetscFunctionBegin;
  PetscCheck(bs >= 1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"block size < 1");
  PetscCheck(n >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
  if (n) PetscValidIntPointer(idx,4);

  PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is),n*bs,is->map->N,bs,&map));
  PetscCall(PetscLayoutDestroy(&is->map));
  is->map = map;

  if (sub->allocated) PetscCall(PetscFree(sub->idx));
  if (mode == PETSC_COPY_VALUES) {
    PetscCall(PetscMalloc1(n,&sub->idx));
    PetscCall(PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt)));
    PetscCall(PetscArraycpy(sub->idx,idx,n));
    sub->allocated = PETSC_TRUE;
  } else if (mode == PETSC_OWN_POINTER) {
    sub->idx = (PetscInt*) idx;
    PetscCall(PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt)));
    sub->allocated = PETSC_TRUE;
  } else if (mode == PETSC_USE_POINTER) {
    sub->idx = (PetscInt*) idx;
    sub->allocated = PETSC_FALSE;
  }

  if (n) {
    min = max = idx[0];
    for (i=1; i<n; i++) {
      if (idx[i] < min) min = idx[i];
      if (idx[i] > max) max = idx[i];
    }
    is->min = bs*min;
    is->max = bs*max+bs-1;
  } else {
    is->min = PETSC_MAX_INT;
    is->max = PETSC_MIN_INT;
  }
  PetscFunctionReturn(0);
}

/*@
   ISCreateBlock - Creates a data structure for an index set containing
   a list of integers. Each integer represents a fixed block size set of indices.

   Collective

   Input Parameters:
+  comm - the MPI communicator
.  bs - number of elements in each block
.  n - the length of the index set (the number of blocks)
.  idx - the list of integers, one for each block, the integers contain the index of the first entry of each block divided by the block size
-  mode - see PetscCopyMode, only PETSC_COPY_VALUES and PETSC_OWN_POINTER are supported in this routine

   Output Parameter:
.  is - the new index set

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

   Example:
   If you wish to index the values {0,1,6,7}, then use
   a block size of 2 and idx of {0,3}.

   Level: beginner

.seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather(), ISBlockSetIndices(), ISBLOCK, ISGENERAL
@*/
PetscErrorCode  ISCreateBlock(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode,IS *is)
{
  PetscFunctionBegin;
  PetscValidPointer(is,6);
  PetscCheck(bs >= 1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"block size < 1");
  PetscCheck(n >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
  if (n) PetscValidIntPointer(idx,4);

  PetscCall(ISCreate(comm,is));
  PetscCall(ISSetType(*is,ISBLOCK));
  PetscCall(ISBlockSetIndices(*is,bs,n,idx,mode));
  PetscFunctionReturn(0);
}

static PetscErrorCode  ISBlockGetIndices_Block(IS is,const PetscInt *idx[])
{
  IS_Block *sub = (IS_Block*)is->data;

  PetscFunctionBegin;
  *idx = sub->idx;
  PetscFunctionReturn(0);
}

static PetscErrorCode  ISBlockRestoreIndices_Block(IS is,const PetscInt *idx[])
{
  PetscFunctionBegin;
  PetscFunctionReturn(0);
}

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

   Not Collective

   Input Parameter:
.  is - the index set

   Output Parameter:
.  idx - the integer indices, one for each block and count of block not indices

   Level: intermediate

.seealso: ISGetIndices(), ISBlockRestoreIndices(), ISBLOCK, ISBlockSetIndices(), ISCreateBlock()
@*/
PetscErrorCode  ISBlockGetIndices(IS is,const PetscInt *idx[])
{
  PetscFunctionBegin;
  PetscUseMethod(is,"ISBlockGetIndices_C",(IS,const PetscInt*[]),(is,idx));
  PetscFunctionReturn(0);
}

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

   Not Collective

   Input Parameter:
.  is - the index set

   Output Parameter:
.  idx - the integer indices

   Level: intermediate

.seealso: ISRestoreIndices(), ISBlockGetIndices()
@*/
PetscErrorCode  ISBlockRestoreIndices(IS is,const PetscInt *idx[])
{
  PetscFunctionBegin;
  PetscUseMethod(is,"ISBlockRestoreIndices_C",(IS,const PetscInt*[]),(is,idx));
  PetscFunctionReturn(0);
}

/*@
   ISBlockGetLocalSize - Returns the local number of blocks in the index set.

   Not Collective

   Input Parameter:
.  is - the index set

   Output Parameter:
.  size - the local number of blocks

   Level: intermediate

.seealso: ISGetBlockSize(), ISBlockGetSize(), ISGetSize(), ISCreateBlock(), ISBLOCK
@*/
PetscErrorCode  ISBlockGetLocalSize(IS is,PetscInt *size)
{
  PetscFunctionBegin;
  PetscUseMethod(is,"ISBlockGetLocalSize_C",(IS,PetscInt*),(is,size));
  PetscFunctionReturn(0);
}

static PetscErrorCode  ISBlockGetLocalSize_Block(IS is,PetscInt *size)
{
  PetscInt       bs, n;

  PetscFunctionBegin;
  PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
  PetscCall(PetscLayoutGetLocalSize(is->map, &n));
  *size = n/bs;
  PetscFunctionReturn(0);
}

/*@
   ISBlockGetSize - Returns the global number of blocks in the index set.

   Not Collective

   Input Parameter:
.  is - the index set

   Output Parameter:
.  size - the global number of blocks

   Level: intermediate

.seealso: ISGetBlockSize(), ISBlockGetLocalSize(), ISGetSize(), ISCreateBlock(), ISBLOCK
@*/
PetscErrorCode  ISBlockGetSize(IS is,PetscInt *size)
{
  PetscFunctionBegin;
  PetscUseMethod(is,"ISBlockGetSize_C",(IS,PetscInt*),(is,size));
  PetscFunctionReturn(0);
}

static PetscErrorCode  ISBlockGetSize_Block(IS is,PetscInt *size)
{
  PetscInt       bs, N;

  PetscFunctionBegin;
  PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
  PetscCall(PetscLayoutGetSize(is->map, &N));
  *size = N/bs;
  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode ISCreate_Block(IS is)
{
  IS_Block       *sub;

  PetscFunctionBegin;
  PetscCall(PetscNewLog(is,&sub));
  is->data = (void *) sub;
  PetscCall(PetscMemcpy(is->ops,&myops,sizeof(myops)));
  PetscCall(PetscObjectComposeFunction((PetscObject)is,"ISBlockSetIndices_C",ISBlockSetIndices_Block));
  PetscCall(PetscObjectComposeFunction((PetscObject)is,"ISBlockGetIndices_C",ISBlockGetIndices_Block));
  PetscCall(PetscObjectComposeFunction((PetscObject)is,"ISBlockRestoreIndices_C",ISBlockRestoreIndices_Block));
  PetscCall(PetscObjectComposeFunction((PetscObject)is,"ISBlockGetSize_C",ISBlockGetSize_Block));
  PetscCall(PetscObjectComposeFunction((PetscObject)is,"ISBlockGetLocalSize_C",ISBlockGetLocalSize_Block));
  PetscFunctionReturn(0);
}