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(PetscObjectComposeFunction((PetscObject)is, "ISShift_C", NULL));
  PetscCall(PetscFree(is->data));
  PetscFunctionReturn(PETSC_SUCCESS);
}

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(PetscLayoutGetLocalSize(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(PETSC_SUCCESS);
}

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(PETSC_SUCCESS);
}

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(PETSC_SUCCESS);
}

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(PETSC_SUCCESS);
}

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(PETSC_SUCCESS);
}

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(PETSC_SUCCESS);
}

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(PETSC_SUCCESS);
}

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

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(PETSC_SUCCESS);
}

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(PETSC_SUCCESS);
}

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(PETSC_SUCCESS);
}

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(PETSC_SUCCESS);
}

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(PETSC_SUCCESS);
}

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(PETSC_SUCCESS);
}

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(PETSC_SUCCESS);
}

static PetscErrorCode ISShift_Block(IS is, PetscInt shift, IS isy)
{
  IS_Block *isb  = (IS_Block *)is->data;
  IS_Block *isby = (IS_Block *)isy->data;
  PetscInt  i, n, bs;

  PetscFunctionBegin;
  PetscCall(PetscLayoutGetLocalSize(is->map, &n));
  PetscCall(PetscLayoutGetBlockSize(is->map, &bs));
  shift /= bs;
  for (i = 0; i < n / bs; i++) isby->idx[i] = isb->idx[i] + shift;
  PetscFunctionReturn(PETSC_SUCCESS);
}

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(PETSC_SUCCESS);
}

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(PETSC_SUCCESS);
}

// clang-format off
static const struct _ISOps myops = {
  PetscDesignatedInitializer(getindices, ISGetIndices_Block),
  PetscDesignatedInitializer(restoreindices, ISRestoreIndices_Block),
  PetscDesignatedInitializer(invertpermutation, ISInvertPermutation_Block),
  PetscDesignatedInitializer(sort, ISSort_Block),
  PetscDesignatedInitializer(sortremovedups, ISSortRemoveDups_Block),
  PetscDesignatedInitializer(sorted, ISSorted_Block),
  PetscDesignatedInitializer(duplicate, ISDuplicate_Block),
  PetscDesignatedInitializer(destroy, ISDestroy_Block),
  PetscDesignatedInitializer(view, ISView_Block),
  PetscDesignatedInitializer(load, ISLoad_Default),
  PetscDesignatedInitializer(copy, ISCopy_Block),
  PetscDesignatedInitializer(togeneral, ISToGeneral_Block),
  PetscDesignatedInitializer(oncomm, ISOnComm_Block),
  PetscDesignatedInitializer(setblocksize, ISSetBlockSize_Block),
  PetscDesignatedInitializer(contiguous, NULL),
  PetscDesignatedInitializer(locate, 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 */
  PetscDesignatedInitializer(sortedlocal, ISSortedLocal_Block),
  PetscDesignatedInitializer(sortedglobal, NULL),
  PetscDesignatedInitializer(uniquelocal, ISUniqueLocal_Block),
  PetscDesignatedInitializer(uniqueglobal, NULL),
  PetscDesignatedInitializer(permlocal, ISPermutationLocal_Block),
  PetscDesignatedInitializer(permglobal, NULL),
  PetscDesignatedInitializer(intervallocal, ISIntervalLocal_Block),
  PetscDesignatedInitializer(intervalglobal, NULL)
};
// clang-format on

/*@
  ISBlockSetIndices - Set integers representing blocks of indices in an index set of `ISType` `ISBLOCK`

  Collective

  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

  Level: beginner

  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.

  The convenience routine `ISCreateBlock()` allows one to create the `IS` and provide the blocks in a single function call.

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

.seealso: [](sec_scatter), `IS`, `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(PETSC_SUCCESS);
}

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) PetscAssertPointer(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(PetscArraycpy(sub->idx, idx, n));
    sub->allocated = PETSC_TRUE;
  } else if (mode == PETSC_OWN_POINTER) {
    sub->idx       = (PetscInt *)idx;
    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(PETSC_SUCCESS);
}

/*@
  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

  Level: beginner

  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.

  The routine `ISBlockSetIndices()` can be used to provide the indices to a preexisting block `IS`

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

.seealso: [](sec_scatter), `IS`, `ISCreateStride()`, `ISCreateGeneral()`, `ISAllGather()`, `ISBlockSetIndices()`, `ISBLOCK`, `ISGENERAL`
@*/
PetscErrorCode ISCreateBlock(MPI_Comm comm, PetscInt bs, PetscInt n, const PetscInt idx[], PetscCopyMode mode, IS *is)
{
  PetscFunctionBegin;
  PetscAssertPointer(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) PetscAssertPointer(idx, 4);

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

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

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

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

/*@C
  ISBlockGetIndices - Gets the indices associated with each block in an `ISBLOCK`

  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

  Note:
  Call `ISBlockRestoreIndices()` when you no longer need access to the indices

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

/*@C
  ISBlockRestoreIndices - Restores the indices associated with each block  in an `ISBLOCK` obtained with `ISBlockGetIndices()`

  Not Collective

  Input Parameter:
. is - the index set

  Output Parameter:
. idx - the integer indices

  Level: intermediate

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

/*@
  ISBlockGetLocalSize - Returns the local number of blocks in the index set of `ISType` `ISBLOCK`

  Not Collective

  Input Parameter:
. is - the index set

  Output Parameter:
. size - the local number of blocks

  Level: intermediate

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

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(PETSC_SUCCESS);
}

/*@
  ISBlockGetSize - Returns the global number of blocks in parallel in the index set of `ISType` `ISBLOCK`

  Not Collective

  Input Parameter:
. is - the index set

  Output Parameter:
. size - the global number of blocks

  Level: intermediate

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

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(PETSC_SUCCESS);
}

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

  PetscFunctionBegin;
  PetscCall(PetscNew(&sub));
  is->data   = (void *)sub;
  is->ops[0] = 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));
  PetscCall(PetscObjectComposeFunction((PetscObject)is, "ISShift_C", ISShift_Block));
  PetscFunctionReturn(PETSC_SUCCESS);
}