xref: /petsc/src/vec/is/utils/isltog.c (revision 633354d9cb03a0dcc5a871fca002dba4b1523513)

#include <petsc/private/isimpl.h> /*I "petscis.h"  I*/
#include <petsc/private/hashmapi.h>
#include <petscsf.h>
#include <petscviewer.h>
#include <petscbt.h>

PetscClassId          IS_LTOGM_CLASSID;
static PetscErrorCode ISLocalToGlobalMappingSetUpBlockInfo_Private(ISLocalToGlobalMapping);

typedef struct {
  PetscInt *globals;
} ISLocalToGlobalMapping_Basic;

typedef struct {
  PetscHMapI globalht;
} ISLocalToGlobalMapping_Hash;

/*@C
  ISGetPointRange - Returns a description of the points in an `IS` suitable for traversal

  Not Collective

  Input Parameter:
. pointIS - The `IS` object

  Output Parameters:
+ pStart - The first index, see notes
. pEnd   - One past the last index, see notes
- points - The indices, see notes

  Level: intermediate

  Notes:
  If the `IS` contains contiguous indices in an `ISSTRIDE`, then the indices are contained in [pStart, pEnd) and points = `NULL`.
  Otherwise, `pStart = 0`, `pEnd = numIndices`, and points is an array of the indices. This supports the following pattern
.vb
  ISGetPointRange(is, &pStart, &pEnd, &points);
  for (p = pStart; p < pEnd; ++p) {
    const PetscInt point = points ? points[p] : p;
    // use point
  }
  ISRestorePointRange(is, &pstart, &pEnd, &points);
.ve
  Hence the same code can be written for `pointIS` being a `ISSTRIDE` or `ISGENERAL`

.seealso: [](sec_scatter), `IS`, `ISRestorePointRange()`, `ISGetPointSubrange()`, `ISGetIndices()`, `ISCreateStride()`
@*/
PetscErrorCode ISGetPointRange(IS pointIS, PetscInt *pStart, PetscInt *pEnd, const PetscInt **points)
{
  PetscInt  numCells, step = 1;
  PetscBool isStride;

  PetscFunctionBeginHot;
  *pStart = 0;
  *points = NULL;
  PetscCall(ISGetLocalSize(pointIS, &numCells));
  PetscCall(PetscObjectTypeCompare((PetscObject)pointIS, ISSTRIDE, &isStride));
  if (isStride) PetscCall(ISStrideGetInfo(pointIS, pStart, &step));
  *pEnd = *pStart + numCells;
  if (!isStride || step != 1) PetscCall(ISGetIndices(pointIS, points));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  ISRestorePointRange - Destroys the traversal description created with `ISGetPointRange()`

  Not Collective

  Input Parameters:
+ pointIS - The `IS` object
. pStart  - The first index, from `ISGetPointRange()`
. pEnd    - One past the last index, from `ISGetPointRange()`
- points  - The indices, from `ISGetPointRange()`

  Level: intermediate

.seealso: [](sec_scatter), `IS`, `ISGetPointRange()`, `ISGetPointSubrange()`, `ISGetIndices()`, `ISCreateStride()`
@*/
PetscErrorCode ISRestorePointRange(IS pointIS, PetscInt *pStart, PetscInt *pEnd, const PetscInt **points)
{
  PetscInt  step = 1;
  PetscBool isStride;

  PetscFunctionBeginHot;
  PetscCall(PetscObjectTypeCompare((PetscObject)pointIS, ISSTRIDE, &isStride));
  if (isStride) PetscCall(ISStrideGetInfo(pointIS, pStart, &step));
  if (!isStride || step != 1) PetscCall(ISGetIndices(pointIS, points));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  ISGetPointSubrange - Configures the input `IS` to be a subrange for the traversal information given

  Not Collective

  Input Parameters:
+ subpointIS - The `IS` object to be configured
. pStart     - The first index of the subrange
. pEnd       - One past the last index for the subrange
- points     - The indices for the entire range, from `ISGetPointRange()`

  Output Parameters:
. subpointIS - The `IS` object now configured to be a subrange

  Level: intermediate

  Note:
  The input `IS` will now respond properly to calls to `ISGetPointRange()` and return the subrange.

.seealso: [](sec_scatter), `IS`, `ISGetPointRange()`, `ISRestorePointRange()`, `ISGetIndices()`, `ISCreateStride()`
@*/
PetscErrorCode ISGetPointSubrange(IS subpointIS, PetscInt pStart, PetscInt pEnd, const PetscInt *points)
{
  PetscFunctionBeginHot;
  if (points) {
    PetscCall(ISSetType(subpointIS, ISGENERAL));
    PetscCall(ISGeneralSetIndices(subpointIS, pEnd - pStart, &points[pStart], PETSC_USE_POINTER));
  } else {
    PetscCall(ISSetType(subpointIS, ISSTRIDE));
    PetscCall(ISStrideSetStride(subpointIS, pEnd - pStart, pStart, 1));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/* -----------------------------------------------------------------------------------------*/

/*
    Creates the global mapping information in the ISLocalToGlobalMapping structure

    If the user has not selected how to handle the global to local mapping then use HASH for "large" problems
*/
static PetscErrorCode ISGlobalToLocalMappingSetUp(ISLocalToGlobalMapping mapping)
{
  PetscInt i, *idx = mapping->indices, n = mapping->n, end, start;

  PetscFunctionBegin;
  if (mapping->data) PetscFunctionReturn(PETSC_SUCCESS);
  end   = 0;
  start = PETSC_MAX_INT;

  for (i = 0; i < n; i++) {
    if (idx[i] < 0) continue;
    if (idx[i] < start) start = idx[i];
    if (idx[i] > end) end = idx[i];
  }
  if (start > end) {
    start = 0;
    end   = -1;
  }
  mapping->globalstart = start;
  mapping->globalend   = end;
  if (!((PetscObject)mapping)->type_name) {
    if ((end - start) > PetscMax(4 * n, 1000000)) {
      PetscCall(ISLocalToGlobalMappingSetType(mapping, ISLOCALTOGLOBALMAPPINGHASH));
    } else {
      PetscCall(ISLocalToGlobalMappingSetType(mapping, ISLOCALTOGLOBALMAPPINGBASIC));
    }
  }
  PetscTryTypeMethod(mapping, globaltolocalmappingsetup);
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode ISGlobalToLocalMappingSetUp_Basic(ISLocalToGlobalMapping mapping)
{
  PetscInt                      i, *idx = mapping->indices, n = mapping->n, end, start, *globals;
  ISLocalToGlobalMapping_Basic *map;

  PetscFunctionBegin;
  start = mapping->globalstart;
  end   = mapping->globalend;
  PetscCall(PetscNew(&map));
  PetscCall(PetscMalloc1(end - start + 2, &globals));
  map->globals = globals;
  for (i = 0; i < end - start + 1; i++) globals[i] = -1;
  for (i = 0; i < n; i++) {
    if (idx[i] < 0) continue;
    globals[idx[i] - start] = i;
  }
  mapping->data = (void *)map;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode ISGlobalToLocalMappingSetUp_Hash(ISLocalToGlobalMapping mapping)
{
  PetscInt                     i, *idx = mapping->indices, n = mapping->n;
  ISLocalToGlobalMapping_Hash *map;

  PetscFunctionBegin;
  PetscCall(PetscNew(&map));
  PetscCall(PetscHMapICreate(&map->globalht));
  for (i = 0; i < n; i++) {
    if (idx[i] < 0) continue;
    PetscCall(PetscHMapISet(map->globalht, idx[i], i));
  }
  mapping->data = (void *)map;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode ISLocalToGlobalMappingDestroy_Basic(ISLocalToGlobalMapping mapping)
{
  ISLocalToGlobalMapping_Basic *map = (ISLocalToGlobalMapping_Basic *)mapping->data;

  PetscFunctionBegin;
  if (!map) PetscFunctionReturn(PETSC_SUCCESS);
  PetscCall(PetscFree(map->globals));
  PetscCall(PetscFree(mapping->data));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode ISLocalToGlobalMappingDestroy_Hash(ISLocalToGlobalMapping mapping)
{
  ISLocalToGlobalMapping_Hash *map = (ISLocalToGlobalMapping_Hash *)mapping->data;

  PetscFunctionBegin;
  if (!map) PetscFunctionReturn(PETSC_SUCCESS);
  PetscCall(PetscHMapIDestroy(&map->globalht));
  PetscCall(PetscFree(mapping->data));
  PetscFunctionReturn(PETSC_SUCCESS);
}

#define GTOLTYPE _Basic
#define GTOLNAME _Basic
#define GTOLBS   mapping->bs
#define GTOL(g, local) \
  do { \
    local = map->globals[g / bs - start]; \
    if (local >= 0) local = bs * local + (g % bs); \
  } while (0)

#include <../src/vec/is/utils/isltog.h>

#define GTOLTYPE _Basic
#define GTOLNAME Block_Basic
#define GTOLBS   1
#define GTOL(g, local) \
  do { \
    local = map->globals[g - start]; \
  } while (0)
#include <../src/vec/is/utils/isltog.h>

#define GTOLTYPE _Hash
#define GTOLNAME _Hash
#define GTOLBS   mapping->bs
#define GTOL(g, local) \
  do { \
    (void)PetscHMapIGet(map->globalht, g / bs, &local); \
    if (local >= 0) local = bs * local + (g % bs); \
  } while (0)
#include <../src/vec/is/utils/isltog.h>

#define GTOLTYPE _Hash
#define GTOLNAME Block_Hash
#define GTOLBS   1
#define GTOL(g, local) \
  do { \
    (void)PetscHMapIGet(map->globalht, g, &local); \
  } while (0)
#include <../src/vec/is/utils/isltog.h>

/*@
  ISLocalToGlobalMappingDuplicate - Duplicates the local to global mapping object

  Not Collective

  Input Parameter:
. ltog - local to global mapping

  Output Parameter:
. nltog - the duplicated local to global mapping

  Level: advanced

.seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreate()`
@*/
PetscErrorCode ISLocalToGlobalMappingDuplicate(ISLocalToGlobalMapping ltog, ISLocalToGlobalMapping *nltog)
{
  ISLocalToGlobalMappingType l2gtype;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ltog, IS_LTOGM_CLASSID, 1);
  PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)ltog), ltog->bs, ltog->n, ltog->indices, PETSC_COPY_VALUES, nltog));
  PetscCall(ISLocalToGlobalMappingGetType(ltog, &l2gtype));
  PetscCall(ISLocalToGlobalMappingSetType(*nltog, l2gtype));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  ISLocalToGlobalMappingGetSize - Gets the local size of a local to global mapping

  Not Collective

  Input Parameter:
. mapping - local to global mapping

  Output Parameter:
. n - the number of entries in the local mapping, `ISLocalToGlobalMappingGetIndices()` returns an array of this length

  Level: advanced

.seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreate()`
@*/
PetscErrorCode ISLocalToGlobalMappingGetSize(ISLocalToGlobalMapping mapping, PetscInt *n)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mapping, IS_LTOGM_CLASSID, 1);
  PetscAssertPointer(n, 2);
  *n = mapping->bs * mapping->n;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  ISLocalToGlobalMappingViewFromOptions - View an `ISLocalToGlobalMapping` based on values in the options database

  Collective

  Input Parameters:
+ A    - the local to global mapping object
. obj  - Optional object that provides the options prefix used for the options database query
- name - command line option

  Level: intermediate

  Note:
  See `PetscObjectViewFromOptions()` for the available `PetscViewer` and `PetscViewerFormat`

.seealso: [](sec_scatter), `PetscViewer`, ``ISLocalToGlobalMapping`, `ISLocalToGlobalMappingView`, `PetscObjectViewFromOptions()`, `ISLocalToGlobalMappingCreate()`
@*/
PetscErrorCode ISLocalToGlobalMappingViewFromOptions(ISLocalToGlobalMapping A, PetscObject obj, const char name[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(A, IS_LTOGM_CLASSID, 1);
  PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  ISLocalToGlobalMappingView - View a local to global mapping

  Not Collective

  Input Parameters:
+ mapping - local to global mapping
- viewer  - viewer

  Level: advanced

.seealso: [](sec_scatter), `PetscViewer`, `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreate()`
@*/
PetscErrorCode ISLocalToGlobalMappingView(ISLocalToGlobalMapping mapping, PetscViewer viewer)
{
  PetscInt    i;
  PetscMPIInt rank;
  PetscBool   iascii;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mapping, IS_LTOGM_CLASSID, 1);
  if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mapping), &viewer));
  PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);

  PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)mapping), &rank));
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
  if (iascii) {
    PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)mapping, viewer));
    PetscCall(PetscViewerASCIIPushSynchronized(viewer));
    for (i = 0; i < mapping->n; i++) {
      PetscInt bs = mapping->bs, g = mapping->indices[i];
      if (bs == 1) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] %" PetscInt_FMT " %" PetscInt_FMT "\n", rank, i, g));
      else PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] %" PetscInt_FMT ":%" PetscInt_FMT " %" PetscInt_FMT ":%" PetscInt_FMT "\n", rank, i * bs, (i + 1) * bs, g * bs, (g + 1) * bs));
    }
    PetscCall(PetscViewerFlush(viewer));
    PetscCall(PetscViewerASCIIPopSynchronized(viewer));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  ISLocalToGlobalMappingCreateIS - Creates a mapping between a local (0 to n)
  ordering and a global parallel ordering.

  Not Collective

  Input Parameter:
. is - index set containing the global numbers for each local number

  Output Parameter:
. mapping - new mapping data structure

  Level: advanced

  Note:
  the block size of the `IS` determines the block size of the mapping

.seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingSetFromOptions()`
@*/
PetscErrorCode ISLocalToGlobalMappingCreateIS(IS is, ISLocalToGlobalMapping *mapping)
{
  PetscInt        n, bs;
  const PetscInt *indices;
  MPI_Comm        comm;
  PetscBool       isblock;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(is, IS_CLASSID, 1);
  PetscAssertPointer(mapping, 2);

  PetscCall(PetscObjectGetComm((PetscObject)is, &comm));
  PetscCall(ISGetLocalSize(is, &n));
  PetscCall(PetscObjectTypeCompare((PetscObject)is, ISBLOCK, &isblock));
  if (!isblock) {
    PetscCall(ISGetIndices(is, &indices));
    PetscCall(ISLocalToGlobalMappingCreate(comm, 1, n, indices, PETSC_COPY_VALUES, mapping));
    PetscCall(ISRestoreIndices(is, &indices));
  } else {
    PetscCall(ISGetBlockSize(is, &bs));
    PetscCall(ISBlockGetIndices(is, &indices));
    PetscCall(ISLocalToGlobalMappingCreate(comm, bs, n / bs, indices, PETSC_COPY_VALUES, mapping));
    PetscCall(ISBlockRestoreIndices(is, &indices));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  ISLocalToGlobalMappingCreateSF - Creates a mapping between a local (0 to n)
  ordering and a global parallel ordering.

  Collective

  Input Parameters:
+ sf    - star forest mapping contiguous local indices to (rank, offset)
- start - first global index on this process, or `PETSC_DECIDE` to compute contiguous global numbering automatically

  Output Parameter:
. mapping - new mapping data structure

  Level: advanced

  Note:
  If a process calls this function with `start` = `PETSC_DECIDE` then all processes must, otherwise the program will hang.

.seealso: [](sec_scatter), `PetscSF`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingSetFromOptions()`
@*/
PetscErrorCode ISLocalToGlobalMappingCreateSF(PetscSF sf, PetscInt start, ISLocalToGlobalMapping *mapping)
{
  PetscInt i, maxlocal, nroots, nleaves, *globals, *ltog;
  MPI_Comm comm;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 1);
  PetscAssertPointer(mapping, 3);
  PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
  PetscCall(PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL));
  if (start == PETSC_DECIDE) {
    start = 0;
    PetscCallMPI(MPI_Exscan(&nroots, &start, 1, MPIU_INT, MPI_SUM, comm));
  } else PetscCheck(start >= 0, comm, PETSC_ERR_ARG_OUTOFRANGE, "start must be nonnegative or PETSC_DECIDE");
  PetscCall(PetscSFGetLeafRange(sf, NULL, &maxlocal));
  ++maxlocal;
  PetscCall(PetscMalloc1(nroots, &globals));
  PetscCall(PetscMalloc1(maxlocal, &ltog));
  for (i = 0; i < nroots; i++) globals[i] = start + i;
  for (i = 0; i < maxlocal; i++) ltog[i] = -1;
  PetscCall(PetscSFBcastBegin(sf, MPIU_INT, globals, ltog, MPI_REPLACE));
  PetscCall(PetscSFBcastEnd(sf, MPIU_INT, globals, ltog, MPI_REPLACE));
  PetscCall(ISLocalToGlobalMappingCreate(comm, 1, maxlocal, ltog, PETSC_OWN_POINTER, mapping));
  PetscCall(PetscFree(globals));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode ISLocalToGlobalMappingResetBlockInfo_Private(ISLocalToGlobalMapping mapping)
{
  PetscFunctionBegin;
  PetscCall(PetscFree(mapping->info_procs));
  PetscCall(PetscFree(mapping->info_numprocs));
  if (mapping->info_indices) {
    for (PetscInt i = 0; i < mapping->info_nproc; i++) PetscCall(PetscFree(mapping->info_indices[i]));
    PetscCall(PetscFree(mapping->info_indices));
  }
  if (mapping->info_nodei) PetscCall(PetscFree(mapping->info_nodei[0]));
  PetscCall(PetscFree2(mapping->info_nodec, mapping->info_nodei));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  ISLocalToGlobalMappingSetBlockSize - Sets the blocksize of the mapping

  Not Collective

  Input Parameters:
+ mapping - mapping data structure
- bs      - the blocksize

  Level: advanced

.seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`
@*/
PetscErrorCode ISLocalToGlobalMappingSetBlockSize(ISLocalToGlobalMapping mapping, PetscInt bs)
{
  PetscInt       *nid;
  const PetscInt *oid;
  PetscInt        i, cn, on, obs, nn;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mapping, IS_LTOGM_CLASSID, 1);
  PetscCheck(bs >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid block size %" PetscInt_FMT, bs);
  if (bs == mapping->bs) PetscFunctionReturn(PETSC_SUCCESS);
  on  = mapping->n;
  obs = mapping->bs;
  oid = mapping->indices;
  nn  = (on * obs) / bs;
  PetscCheck((on * obs) % bs == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Block size %" PetscInt_FMT " is inconsistent with block size %" PetscInt_FMT " and number of block indices %" PetscInt_FMT, bs, obs, on);

  PetscCall(PetscMalloc1(nn, &nid));
  PetscCall(ISLocalToGlobalMappingGetIndices(mapping, &oid));
  for (i = 0; i < nn; i++) {
    PetscInt j;
    for (j = 0, cn = 0; j < bs - 1; j++) {
      if (oid[i * bs + j] < 0) {
        cn++;
        continue;
      }
      PetscCheck(oid[i * bs + j] == oid[i * bs + j + 1] - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Block sizes %" PetscInt_FMT " and %" PetscInt_FMT " are incompatible with the block indices: non consecutive indices %" PetscInt_FMT " %" PetscInt_FMT, bs, obs, oid[i * bs + j], oid[i * bs + j + 1]);
    }
    if (oid[i * bs + j] < 0) cn++;
    if (cn) {
      PetscCheck(cn == bs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Block sizes %" PetscInt_FMT " and %" PetscInt_FMT " are incompatible with the block indices: invalid number of negative entries in block %" PetscInt_FMT, bs, obs, cn);
      nid[i] = -1;
    } else {
      nid[i] = oid[i * bs] / bs;
    }
  }
  PetscCall(ISLocalToGlobalMappingRestoreIndices(mapping, &oid));

  mapping->n  = nn;
  mapping->bs = bs;
  PetscCall(PetscFree(mapping->indices));
  mapping->indices     = nid;
  mapping->globalstart = 0;
  mapping->globalend   = 0;

  /* reset the cached information */
  PetscCall(ISLocalToGlobalMappingResetBlockInfo_Private(mapping));
  PetscTryTypeMethod(mapping, destroy);
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  ISLocalToGlobalMappingGetBlockSize - Gets the blocksize of the mapping
  ordering and a global parallel ordering.

  Not Collective

  Input Parameter:
. mapping - mapping data structure

  Output Parameter:
. bs - the blocksize

  Level: advanced

.seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`
@*/
PetscErrorCode ISLocalToGlobalMappingGetBlockSize(ISLocalToGlobalMapping mapping, PetscInt *bs)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mapping, IS_LTOGM_CLASSID, 1);
  *bs = mapping->bs;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  ISLocalToGlobalMappingCreate - Creates a mapping between a local (0 to n)
  ordering and a global parallel ordering.

  Not Collective, but communicator may have more than one process

  Input Parameters:
+ comm    - MPI communicator
. bs      - the block size
. n       - the number of local elements divided by the block size, or equivalently the number of block indices
. indices - the global index for each local element, these do not need to be in increasing order (sorted), these values should not be scaled (i.e. multiplied) by the blocksize bs
- mode    - see PetscCopyMode

  Output Parameter:
. mapping - new mapping data structure

  Level: advanced

  Notes:
  There is one integer value in indices per block and it represents the actual indices bs*idx + j, where j=0,..,bs-1

  For "small" problems when using `ISGlobalToLocalMappingApply()` and `ISGlobalToLocalMappingApplyBlock()`, the `ISLocalToGlobalMappingType`
  of `ISLOCALTOGLOBALMAPPINGBASIC` will be used; this uses more memory but is faster; this approach is not scalable for extremely large mappings.

  For large problems `ISLOCALTOGLOBALMAPPINGHASH` is used, this is scalable.
  Use `ISLocalToGlobalMappingSetType()` or call `ISLocalToGlobalMappingSetFromOptions()` with the option
  `-islocaltoglobalmapping_type` <`basic`,`hash`> to control which is used.

.seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingSetFromOptions()`,
          `ISLOCALTOGLOBALMAPPINGBASIC`, `ISLOCALTOGLOBALMAPPINGHASH`
          `ISLocalToGlobalMappingSetType()`, `ISLocalToGlobalMappingType`
@*/
PetscErrorCode ISLocalToGlobalMappingCreate(MPI_Comm comm, PetscInt bs, PetscInt n, const PetscInt indices[], PetscCopyMode mode, ISLocalToGlobalMapping *mapping)
{
  PetscInt *in;

  PetscFunctionBegin;
  if (n) PetscAssertPointer(indices, 4);
  PetscAssertPointer(mapping, 6);

  *mapping = NULL;
  PetscCall(ISInitializePackage());

  PetscCall(PetscHeaderCreate(*mapping, IS_LTOGM_CLASSID, "ISLocalToGlobalMapping", "Local to global mapping", "IS", comm, ISLocalToGlobalMappingDestroy, ISLocalToGlobalMappingView));
  (*mapping)->n  = n;
  (*mapping)->bs = bs;
  if (mode == PETSC_COPY_VALUES) {
    PetscCall(PetscMalloc1(n, &in));
    PetscCall(PetscArraycpy(in, indices, n));
    (*mapping)->indices         = in;
    (*mapping)->dealloc_indices = PETSC_TRUE;
  } else if (mode == PETSC_OWN_POINTER) {
    (*mapping)->indices         = (PetscInt *)indices;
    (*mapping)->dealloc_indices = PETSC_TRUE;
  } else if (mode == PETSC_USE_POINTER) {
    (*mapping)->indices = (PetscInt *)indices;
  } else SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid mode %d", mode);
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscFunctionList ISLocalToGlobalMappingList = NULL;

/*@
  ISLocalToGlobalMappingSetFromOptions - Set mapping options from the options database.

  Not Collective

  Input Parameter:
. mapping - mapping data structure

  Options Database Key:
. -islocaltoglobalmapping_type - <basic,hash> nonscalable and scalable versions

  Level: advanced

.seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingDestroy()`,
`ISLocalToGlobalMappingCreateIS()`, `ISLOCALTOGLOBALMAPPINGBASIC`,
`ISLOCALTOGLOBALMAPPINGHASH`, `ISLocalToGlobalMappingSetType()`, `ISLocalToGlobalMappingType`
@*/
PetscErrorCode ISLocalToGlobalMappingSetFromOptions(ISLocalToGlobalMapping mapping)
{
  char                       type[256];
  ISLocalToGlobalMappingType defaulttype = "Not set";
  PetscBool                  flg;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mapping, IS_LTOGM_CLASSID, 1);
  PetscCall(ISLocalToGlobalMappingRegisterAll());
  PetscObjectOptionsBegin((PetscObject)mapping);
  PetscCall(PetscOptionsFList("-islocaltoglobalmapping_type", "ISLocalToGlobalMapping method", "ISLocalToGlobalMappingSetType", ISLocalToGlobalMappingList, (char *)(((PetscObject)mapping)->type_name) ? ((PetscObject)mapping)->type_name : defaulttype, type, 256, &flg));
  if (flg) PetscCall(ISLocalToGlobalMappingSetType(mapping, type));
  PetscOptionsEnd();
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  ISLocalToGlobalMappingDestroy - Destroys a mapping between a local (0 to n)
  ordering and a global parallel ordering.

  Not Collective

  Input Parameter:
. mapping - mapping data structure

  Level: advanced

.seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingCreate()`
@*/
PetscErrorCode ISLocalToGlobalMappingDestroy(ISLocalToGlobalMapping *mapping)
{
  PetscFunctionBegin;
  if (!*mapping) PetscFunctionReturn(PETSC_SUCCESS);
  PetscValidHeaderSpecific(*mapping, IS_LTOGM_CLASSID, 1);
  if (--((PetscObject)*mapping)->refct > 0) {
    *mapping = NULL;
    PetscFunctionReturn(PETSC_SUCCESS);
  }
  if ((*mapping)->dealloc_indices) PetscCall(PetscFree((*mapping)->indices));
  PetscCall(ISLocalToGlobalMappingResetBlockInfo_Private(*mapping));
  PetscTryTypeMethod(*mapping, destroy);
  PetscCall(PetscHeaderDestroy(mapping));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  ISLocalToGlobalMappingApplyIS - Creates from an `IS` in the local numbering
  a new index set using the global numbering defined in an `ISLocalToGlobalMapping`
  context.

  Collective

  Input Parameters:
+ mapping - mapping between local and global numbering
- is      - index set in local numbering

  Output Parameter:
. newis - index set in global numbering

  Level: advanced

  Note:
  The output `IS` will have the same communicator as the input `IS`.

.seealso: [](sec_scatter), `ISLocalToGlobalMappingApply()`, `ISLocalToGlobalMappingCreate()`,
          `ISLocalToGlobalMappingDestroy()`, `ISGlobalToLocalMappingApply()`
@*/
PetscErrorCode ISLocalToGlobalMappingApplyIS(ISLocalToGlobalMapping mapping, IS is, IS *newis)
{
  PetscInt        n, *idxout;
  const PetscInt *idxin;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mapping, IS_LTOGM_CLASSID, 1);
  PetscValidHeaderSpecific(is, IS_CLASSID, 2);
  PetscAssertPointer(newis, 3);

  PetscCall(ISGetLocalSize(is, &n));
  PetscCall(ISGetIndices(is, &idxin));
  PetscCall(PetscMalloc1(n, &idxout));
  PetscCall(ISLocalToGlobalMappingApply(mapping, n, idxin, idxout));
  PetscCall(ISRestoreIndices(is, &idxin));
  PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is), n, idxout, PETSC_OWN_POINTER, newis));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  ISLocalToGlobalMappingApply - Takes a list of integers in a local numbering
  and converts them to the global numbering.

  Not Collective

  Input Parameters:
+ mapping - the local to global mapping context
. N       - number of integers
- in      - input indices in local numbering

  Output Parameter:
. out - indices in global numbering

  Level: advanced

  Note:
  The `in` and `out` array parameters may be identical.

.seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingApplyBlock()`, `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingDestroy()`,
          `ISLocalToGlobalMappingApplyIS()`, `AOCreateBasic()`, `AOApplicationToPetsc()`,
          `AOPetscToApplication()`, `ISGlobalToLocalMappingApply()`
@*/
PetscErrorCode ISLocalToGlobalMappingApply(ISLocalToGlobalMapping mapping, PetscInt N, const PetscInt in[], PetscInt out[])
{
  PetscInt i, bs, Nmax;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mapping, IS_LTOGM_CLASSID, 1);
  bs   = mapping->bs;
  Nmax = bs * mapping->n;
  if (bs == 1) {
    const PetscInt *idx = mapping->indices;
    for (i = 0; i < N; i++) {
      if (in[i] < 0) {
        out[i] = in[i];
        continue;
      }
      PetscCheck(in[i] < Nmax, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local index %" PetscInt_FMT " too large %" PetscInt_FMT " (max) at %" PetscInt_FMT, in[i], Nmax - 1, i);
      out[i] = idx[in[i]];
    }
  } else {
    const PetscInt *idx = mapping->indices;
    for (i = 0; i < N; i++) {
      if (in[i] < 0) {
        out[i] = in[i];
        continue;
      }
      PetscCheck(in[i] < Nmax, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local index %" PetscInt_FMT " too large %" PetscInt_FMT " (max) at %" PetscInt_FMT, in[i], Nmax - 1, i);
      out[i] = idx[in[i] / bs] * bs + (in[i] % bs);
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  ISLocalToGlobalMappingApplyBlock - Takes a list of integers in a local block numbering and converts them to the global block numbering

  Not Collective

  Input Parameters:
+ mapping - the local to global mapping context
. N       - number of integers
- in      - input indices in local block numbering

  Output Parameter:
. out - indices in global block numbering

  Example:
  If the index values are {0,1,6,7} set with a call to `ISLocalToGlobalMappingCreate`(`PETSC_COMM_SELF`,2,2,{0,3}) then the mapping applied to 0
  (the first block) would produce 0 and the mapping applied to 1 (the second block) would produce 3.

  Level: advanced

  Note:
  The `in` and `out` array parameters may be identical.

.seealso: [](sec_scatter), `ISLocalToGlobalMappingApply()`, `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingDestroy()`,
          `ISLocalToGlobalMappingApplyIS()`, `AOCreateBasic()`, `AOApplicationToPetsc()`,
          `AOPetscToApplication()`, `ISGlobalToLocalMappingApply()`
@*/
PetscErrorCode ISLocalToGlobalMappingApplyBlock(ISLocalToGlobalMapping mapping, PetscInt N, const PetscInt in[], PetscInt out[])
{
  PetscInt        i, Nmax;
  const PetscInt *idx;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mapping, IS_LTOGM_CLASSID, 1);
  Nmax = mapping->n;
  idx  = mapping->indices;
  for (i = 0; i < N; i++) {
    if (in[i] < 0) {
      out[i] = in[i];
      continue;
    }
    PetscCheck(in[i] < Nmax, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local block index %" PetscInt_FMT " too large %" PetscInt_FMT " (max) at %" PetscInt_FMT, in[i], Nmax - 1, i);
    out[i] = idx[in[i]];
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  ISGlobalToLocalMappingApply - Provides the local numbering for a list of integers
  specified with a global numbering.

  Not Collective

  Input Parameters:
+ mapping - mapping between local and global numbering
. type    - `IS_GTOLM_MASK` - maps global indices with no local value to -1 in the output list (i.e., mask them)
           `IS_GTOLM_DROP` - drops the indices with no local value from the output list
. n       - number of global indices to map
- idx     - global indices to map

  Output Parameters:
+ nout   - number of indices in output array (if type == `IS_GTOLM_MASK` then nout = n)
- idxout - local index of each global index, one must pass in an array long enough
             to hold all the indices. You can call `ISGlobalToLocalMappingApply()` with
             idxout == NULL to determine the required length (returned in nout)
             and then allocate the required space and call `ISGlobalToLocalMappingApply()`
             a second time to set the values.

  Level: advanced

  Notes:
  Either `nout` or `idxout` may be `NULL`. `idx` and `idxout` may be identical.

  For "small" problems when using `ISGlobalToLocalMappingApply()` and `ISGlobalToLocalMappingApplyBlock()`, the `ISLocalToGlobalMappingType` of
  `ISLOCALTOGLOBALMAPPINGBASIC` will be used;
  this uses more memory but is faster; this approach is not scalable for extremely large mappings. For large problems `ISLOCALTOGLOBALMAPPINGHASH` is used, this is scalable.
  Use `ISLocalToGlobalMappingSetType()` or call `ISLocalToGlobalMappingSetFromOptions()` with the option -islocaltoglobalmapping_type <basic,hash> to control which is used.

  Developer Notes:
  The manual page states that `idx` and `idxout` may be identical but the calling
  sequence declares `idx` as const so it cannot be the same as `idxout`.

.seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingApply()`, `ISGlobalToLocalMappingApplyBlock()`, `ISLocalToGlobalMappingCreate()`,
          `ISLocalToGlobalMappingDestroy()`
@*/
PetscErrorCode ISGlobalToLocalMappingApply(ISLocalToGlobalMapping mapping, ISGlobalToLocalMappingMode type, PetscInt n, const PetscInt idx[], PetscInt *nout, PetscInt idxout[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mapping, IS_LTOGM_CLASSID, 1);
  if (!mapping->data) PetscCall(ISGlobalToLocalMappingSetUp(mapping));
  PetscUseTypeMethod(mapping, globaltolocalmappingapply, type, n, idx, nout, idxout);
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  ISGlobalToLocalMappingApplyIS - Creates from an `IS` in the global numbering
  a new index set using the local numbering defined in an `ISLocalToGlobalMapping`
  context.

  Not Collective

  Input Parameters:
+ mapping - mapping between local and global numbering
. type    - `IS_GTOLM_MASK` - maps global indices with no local value to -1 in the output list (i.e., mask them)
           `IS_GTOLM_DROP` - drops the indices with no local value from the output list
- is      - index set in global numbering

  Output Parameter:
. newis - index set in local numbering

  Level: advanced

  Note:
  The output `IS` will be sequential, as it encodes a purely local operation

.seealso: [](sec_scatter), `ISGlobalToLocalMapping`, `ISGlobalToLocalMappingApply()`, `ISLocalToGlobalMappingCreate()`,
          `ISLocalToGlobalMappingDestroy()`
@*/
PetscErrorCode ISGlobalToLocalMappingApplyIS(ISLocalToGlobalMapping mapping, ISGlobalToLocalMappingMode type, IS is, IS *newis)
{
  PetscInt        n, nout, *idxout;
  const PetscInt *idxin;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mapping, IS_LTOGM_CLASSID, 1);
  PetscValidHeaderSpecific(is, IS_CLASSID, 3);
  PetscAssertPointer(newis, 4);

  PetscCall(ISGetLocalSize(is, &n));
  PetscCall(ISGetIndices(is, &idxin));
  if (type == IS_GTOLM_MASK) {
    PetscCall(PetscMalloc1(n, &idxout));
  } else {
    PetscCall(ISGlobalToLocalMappingApply(mapping, type, n, idxin, &nout, NULL));
    PetscCall(PetscMalloc1(nout, &idxout));
  }
  PetscCall(ISGlobalToLocalMappingApply(mapping, type, n, idxin, &nout, idxout));
  PetscCall(ISRestoreIndices(is, &idxin));
  PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nout, idxout, PETSC_OWN_POINTER, newis));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  ISGlobalToLocalMappingApplyBlock - Provides the local block numbering for a list of integers
  specified with a block global numbering.

  Not Collective

  Input Parameters:
+ mapping - mapping between local and global numbering
. type    - `IS_GTOLM_MASK` - maps global indices with no local value to -1 in the output list (i.e., mask them)
           `IS_GTOLM_DROP` - drops the indices with no local value from the output list
. n       - number of global indices to map
- idx     - global indices to map

  Output Parameters:
+ nout   - number of indices in output array (if type == `IS_GTOLM_MASK` then nout = n)
- idxout - local index of each global index, one must pass in an array long enough
             to hold all the indices. You can call `ISGlobalToLocalMappingApplyBlock()` with
             idxout == NULL to determine the required length (returned in nout)
             and then allocate the required space and call `ISGlobalToLocalMappingApplyBlock()`
             a second time to set the values.

  Level: advanced

  Notes:
  Either `nout` or `idxout` may be `NULL`. `idx` and `idxout` may be identical.

  For "small" problems when using `ISGlobalToLocalMappingApply()` and `ISGlobalToLocalMappingApplyBlock()`, the `ISLocalToGlobalMappingType` of
  `ISLOCALTOGLOBALMAPPINGBASIC` will be used;
  this uses more memory but is faster; this approach is not scalable for extremely large mappings. For large problems `ISLOCALTOGLOBALMAPPINGHASH` is used, this is scalable.
  Use `ISLocalToGlobalMappingSetType()` or call `ISLocalToGlobalMappingSetFromOptions()` with the option -islocaltoglobalmapping_type <basic,hash> to control which is used.

  Developer Notes:
  The manual page states that `idx` and `idxout` may be identical but the calling
  sequence declares `idx` as const so it cannot be the same as `idxout`.

.seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingApply()`, `ISGlobalToLocalMappingApply()`, `ISLocalToGlobalMappingCreate()`,
          `ISLocalToGlobalMappingDestroy()`
@*/
PetscErrorCode ISGlobalToLocalMappingApplyBlock(ISLocalToGlobalMapping mapping, ISGlobalToLocalMappingMode type, PetscInt n, const PetscInt idx[], PetscInt *nout, PetscInt idxout[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mapping, IS_LTOGM_CLASSID, 1);
  if (!mapping->data) PetscCall(ISGlobalToLocalMappingSetUp(mapping));
  PetscUseTypeMethod(mapping, globaltolocalmappingapplyblock, type, n, idx, nout, idxout);
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  ISLocalToGlobalMappingGetBlockInfo - Gets the neighbor information

  Collective the first time it is called

  Input Parameter:
. mapping - the mapping from local to global indexing

  Output Parameters:
+ nproc    - number of processes that are connected to the calling process
. procs    - neighboring processes
. numprocs - number of block indices for each process
- indices  - block indices (in local numbering) shared with neighbors (sorted by global numbering)

  Level: advanced

.seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
          `ISLocalToGlobalMappingRestoreBlockInfo()`
@*/
PetscErrorCode ISLocalToGlobalMappingGetBlockInfo(ISLocalToGlobalMapping mapping, PetscInt *nproc, PetscInt *procs[], PetscInt *numprocs[], PetscInt **indices[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mapping, IS_LTOGM_CLASSID, 1);
  PetscCall(ISLocalToGlobalMappingSetUpBlockInfo_Private(mapping));
  if (nproc) *nproc = mapping->info_nproc;
  if (procs) *procs = mapping->info_procs;
  if (numprocs) *numprocs = mapping->info_numprocs;
  if (indices) *indices = mapping->info_indices;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  ISLocalToGlobalMappingGetBlockNodeInfo - Gets the neighbor information for each local block index

  Collective the first time it is called

  Input Parameter:
. mapping - the mapping from local to global indexing

  Output Parameter:
+ n       - number of local block nodes
. n_procs - an array storing the number of processes for each local block node (including self)
- procs   - the processes' rank for each local block node (sorted, self is first)

  Level: advanced

  Notes:
  The user needs to call `ISLocalToGlobalMappingRestoreBlockNodeInfo()` when the data is no longer needed.
  The information returned by this function complements that of `ISLocalToGlobalMappingGetBlockInfo()`.
  The latter only provides local information, and the neighboring information
  cannot be inferred in the general case, unless the mapping is locally one-to-one on each process.

.seealso: `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
          `ISLocalToGlobalMappingGetBlockInfo()`, `ISLocalToGlobalMappingRestoreBlockNodeInfo()`, `ISLocalToGlobalMappingGetNodeInfo()`
@*/
PetscErrorCode ISLocalToGlobalMappingGetBlockNodeInfo(ISLocalToGlobalMapping mapping, PetscInt *n, PetscInt *n_procs[], PetscInt **procs[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mapping, IS_LTOGM_CLASSID, 1);
  PetscCall(ISLocalToGlobalMappingSetUpBlockInfo_Private(mapping));
  if (n) *n = mapping->n;
  if (n_procs) *n_procs = mapping->info_nodec;
  if (procs) *procs = mapping->info_nodei;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  ISLocalToGlobalMappingRestoreBlockNodeInfo - Frees the memory allocated by `ISLocalToGlobalMappingGetBlockNodeInfo()`

  Not Collective

  Input Parameters:
+ mapping - the mapping from local to global indexing
. n       - number of local block nodes
. n_procs - an array storing the number of processes for each local block nodes (including self)
- procs   - the processes' rank for each local block node (sorted, self is first)

  Level: advanced

.seealso: `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
          `ISLocalToGlobalMappingGetBlockNodeInfo()`
@*/
PetscErrorCode ISLocalToGlobalMappingRestoreBlockNodeInfo(ISLocalToGlobalMapping mapping, PetscInt *n, PetscInt *n_procs[], PetscInt **procs[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mapping, IS_LTOGM_CLASSID, 1);
  if (n) *n = 0;
  if (n_procs) *n_procs = NULL;
  if (procs) *procs = NULL;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode ISLocalToGlobalMappingSetUpBlockInfo_Private(ISLocalToGlobalMapping mapping)
{
  PetscSF            sf;
  MPI_Comm           comm;
  const PetscSFNode *sfnode;
  PetscSFNode       *newsfnode;
  PetscLayout        layout;
  PetscHMapI         neighs;
  PetscHashIter      iter;
  PetscBool          missing;
  const PetscInt    *gidxs, *rootdegree;
  PetscInt          *mask, *mrootdata, *leafdata, *newleafdata, *leafrd, *tmpg;
  PetscInt           nroots, nleaves, newnleaves, bs, i, j, m, mnroots, p;
  PetscMPIInt        rank, size;

  PetscFunctionBegin;
  if (mapping->info_numprocs) PetscFunctionReturn(PETSC_SUCCESS);
  PetscCall(PetscObjectGetComm((PetscObject)mapping, &comm));
  PetscCallMPI(MPI_Comm_size(comm, &size));
  PetscCallMPI(MPI_Comm_rank(comm, &rank));

  /* Get mapping indices */
  PetscCall(ISLocalToGlobalMappingGetBlockSize(mapping, &bs));
  PetscCall(ISLocalToGlobalMappingGetBlockIndices(mapping, &gidxs));
  PetscCall(ISLocalToGlobalMappingGetSize(mapping, &nleaves));
  nleaves /= bs;

  /* Create layout for global indices */
  for (i = 0, m = 0; i < nleaves; i++) m = PetscMax(m, gidxs[i]);
  PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &m, 1, MPIU_INT, MPI_MAX, comm));
  PetscCall(PetscLayoutCreate(comm, &layout));
  PetscCall(PetscLayoutSetSize(layout, m + 1));
  PetscCall(PetscLayoutSetUp(layout));

  /* Create SF to share global indices */
  PetscCall(PetscSFCreate(comm, &sf));
  PetscCall(PetscSFSetGraphLayout(sf, layout, nleaves, NULL, PETSC_OWN_POINTER, gidxs));
  PetscCall(PetscSFSetUp(sf));
  PetscCall(PetscLayoutDestroy(&layout));

  /* communicate root degree to leaves */
  PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, &sfnode));
  PetscCall(PetscSFComputeDegreeBegin(sf, &rootdegree));
  PetscCall(PetscSFComputeDegreeEnd(sf, &rootdegree));
  for (i = 0, mnroots = 0; i < nroots; i++) mnroots += rootdegree[i];
  PetscCall(PetscMalloc3(2 * PetscMax(mnroots, nroots), &mrootdata, 2 * nleaves, &leafdata, nleaves, &leafrd));
  for (i = 0, m = 0; i < nroots; i++) {
    mrootdata[2 * i + 0] = rootdegree[i];
    mrootdata[2 * i + 1] = m;
    m += rootdegree[i];
  }
  PetscCall(PetscSFBcastBegin(sf, MPIU_2INT, mrootdata, leafdata, MPI_REPLACE));
  PetscCall(PetscSFBcastEnd(sf, MPIU_2INT, mrootdata, leafdata, MPI_REPLACE));

  /* allocate enough space to store ranks */
  for (i = 0, newnleaves = 0; i < nleaves; i++) {
    newnleaves += leafdata[2 * i];
    leafrd[i] = leafdata[2 * i];
  }

  /* create new SF nodes to collect multi-root data at leaves */
  PetscCall(PetscMalloc1(newnleaves, &newsfnode));
  for (i = 0, m = 0; i < nleaves; i++) {
    for (j = 0; j < leafrd[i]; j++) {
      newsfnode[m].rank  = sfnode[i].rank;
      newsfnode[m].index = leafdata[2 * i + 1] + j;
      m++;
    }
  }

  /* gather ranks at multi roots */
  for (i = 0; i < mnroots; i++) mrootdata[i] = -1;
  for (i = 0; i < nleaves; i++) leafdata[i] = (PetscInt)rank;

  PetscCall(PetscSFGatherBegin(sf, MPIU_INT, leafdata, mrootdata));
  PetscCall(PetscSFGatherEnd(sf, MPIU_INT, leafdata, mrootdata));

  /* set new multi-leaves graph into the SF */
  PetscCall(PetscSFSetGraph(sf, mnroots, newnleaves, NULL, PETSC_OWN_POINTER, newsfnode, PETSC_OWN_POINTER));
  PetscCall(PetscSFSetUp(sf));

  /* broadcast multi-root data to multi-leaves */
  PetscCall(PetscMalloc1(newnleaves, &newleafdata));
  PetscCall(PetscSFBcastBegin(sf, MPIU_INT, mrootdata, newleafdata, MPI_REPLACE));
  PetscCall(PetscSFBcastEnd(sf, MPIU_INT, mrootdata, newleafdata, MPI_REPLACE));

  /* sort sharing ranks */
  for (i = 0, m = 0; i < nleaves; i++) {
    PetscCall(PetscSortInt(leafrd[i], newleafdata + m));
    m += leafrd[i];
  }

  /* Number of neighbors and their ranks */
  PetscCall(PetscHMapICreate(&neighs));
  for (i = 0; i < newnleaves; i++) PetscCall(PetscHMapIPut(neighs, newleafdata[i], &iter, &missing));
  PetscCall(PetscHMapIGetSize(neighs, &mapping->info_nproc));
  PetscCall(PetscMalloc1(mapping->info_nproc + 1, &mapping->info_procs));
  PetscCall(PetscHMapIGetKeys(neighs, (i = 0, &i), mapping->info_procs));
  for (i = 0; i < mapping->info_nproc; i++) { /* put info for self first */
    if (mapping->info_procs[i] == rank) {
      PetscInt newr = mapping->info_procs[0];

      mapping->info_procs[0] = rank;
      mapping->info_procs[i] = newr;
      break;
    }
  }
  if (mapping->info_nproc) PetscCall(PetscSortInt(mapping->info_nproc - 1, mapping->info_procs + 1));
  PetscCall(PetscHMapIDestroy(&neighs));

  /* collect info data */
  PetscCall(PetscMalloc1(mapping->info_nproc + 1, &mapping->info_numprocs));
  PetscCall(PetscMalloc1(mapping->info_nproc + 1, &mapping->info_indices));
  for (i = 0; i < mapping->info_nproc + 1; i++) mapping->info_indices[i] = NULL;

  PetscCall(PetscMalloc1(nleaves, &mask));
  PetscCall(PetscMalloc1(nleaves, &tmpg));
  for (p = 0; p < mapping->info_nproc; p++) {
    PetscInt *tmp, trank = mapping->info_procs[p];

    PetscCall(PetscMemzero(mask, nleaves * sizeof(*mask)));
    for (i = 0, m = 0; i < nleaves; i++) {
      for (j = 0; j < leafrd[i]; j++) {
        if (newleafdata[m] == trank) mask[i]++;
        if (!p && newleafdata[m] != rank) mask[i]++;
        m++;
      }
    }
    for (i = 0, m = 0; i < nleaves; i++)
      if (mask[i] > (!p ? 1 : 0)) m++;

    PetscCall(PetscMalloc1(m, &tmp));
    for (i = 0, m = 0; i < nleaves; i++)
      if (mask[i] > (!p ? 1 : 0)) {
        tmp[m]  = i;
        tmpg[m] = gidxs[i];
        m++;
      }
    PetscCall(PetscSortIntWithArray(m, tmpg, tmp));
    mapping->info_indices[p]  = tmp;
    mapping->info_numprocs[p] = m;
  }

  /* Node info */
  PetscCall(PetscMalloc2(nleaves, &mapping->info_nodec, nleaves + 1, &mapping->info_nodei));
  PetscCall(PetscArraycpy(mapping->info_nodec, leafrd, nleaves));
  PetscCall(PetscMalloc1(newnleaves, &mapping->info_nodei[0]));
  for (i = 0; i < nleaves - 1; i++) mapping->info_nodei[i + 1] = mapping->info_nodei[i] + mapping->info_nodec[i];
  PetscCall(PetscArraycpy(mapping->info_nodei[0], newleafdata, newnleaves));

  PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(mapping, &gidxs));
  PetscCall(PetscFree(tmpg));
  PetscCall(PetscFree(mask));
  PetscCall(PetscSFDestroy(&sf));
  PetscCall(PetscFree3(mrootdata, leafdata, leafrd));
  PetscCall(PetscFree(newleafdata));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  ISLocalToGlobalMappingRestoreBlockInfo - Frees the memory allocated by `ISLocalToGlobalMappingGetBlockInfo()`

  Not Collective

  Input Parameters:
+ mapping  - the mapping from local to global indexing
. nproc    - number of processes that are connected to the calling process
. procs    - neighboring processes
. numprocs - number of block indices for each process
- indices  - block indices (in local numbering) shared with neighbors (sorted by global numbering)

  Level: advanced

.seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
          `ISLocalToGlobalMappingGetInfo()`
@*/
PetscErrorCode ISLocalToGlobalMappingRestoreBlockInfo(ISLocalToGlobalMapping mapping, PetscInt *nproc, PetscInt *procs[], PetscInt *numprocs[], PetscInt **indices[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mapping, IS_LTOGM_CLASSID, 1);
  if (nproc) *nproc = 0;
  if (procs) *procs = NULL;
  if (numprocs) *numprocs = NULL;
  if (indices) *indices = NULL;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  ISLocalToGlobalMappingGetInfo - Gets the neighbor information for each process

  Collective the first time it is called

  Input Parameter:
. mapping - the mapping from local to global indexing

  Output Parameters:
+ nproc    - number of processes that are connected to the calling process
. procs    - neighboring processes
. numprocs - number of indices for each process
- indices  - indices (in local numbering) shared with neighbors (sorted by global numbering)

  Level: advanced

  Note:
  The user needs to call `ISLocalToGlobalMappingRestoreInfo()` when the data is no longer needed.

  Fortran Notes:
  There is no `ISLocalToGlobalMappingRestoreInfo()` in Fortran. You must make sure that
  `procs`[], `numprocs`[] and `indices`[][] are large enough arrays, either by allocating them
  dynamically or defining static ones large enough.

.seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
          `ISLocalToGlobalMappingRestoreInfo()`, `ISLocalToGlobalMappingGetNodeInfo()`
@*/
PetscErrorCode ISLocalToGlobalMappingGetInfo(ISLocalToGlobalMapping mapping, PetscInt *nproc, PetscInt *procs[], PetscInt *numprocs[], PetscInt **indices[])
{
  PetscInt **bindices = NULL, *bnumprocs = NULL, bs, i, j, k, n, *bprocs;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mapping, IS_LTOGM_CLASSID, 1);
  bs = mapping->bs;
  PetscCall(ISLocalToGlobalMappingGetBlockInfo(mapping, &n, &bprocs, &bnumprocs, &bindices));
  if (bs > 1) { /* we need to expand the cached info */
    if (indices) PetscCall(PetscCalloc1(n, indices));
    if (numprocs) PetscCall(PetscCalloc1(n, numprocs));
    if (indices || numprocs) {
      for (i = 0; i < n; i++) {
        if (indices) {
          PetscCall(PetscMalloc1(bs * bnumprocs[i], &(*indices)[i]));
          for (j = 0; j < bnumprocs[i]; j++) {
            for (k = 0; k < bs; k++) (*indices)[i][j * bs + k] = bs * bindices[i][j] + k;
          }
        }
        if (numprocs) (*numprocs)[i] = bnumprocs[i] * bs;
      }
    }
  } else {
    if (numprocs) *numprocs = bnumprocs;
    if (indices) *indices = bindices;
  }
  if (nproc) *nproc = n;
  if (procs) *procs = bprocs;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  ISLocalToGlobalMappingRestoreInfo - Frees the memory allocated by `ISLocalToGlobalMappingGetInfo()`

  Not Collective

  Input Parameters:
+ mapping  - the mapping from local to global indexing
. nproc    - number of processes that are connected to the calling process
. procs    - neighboring processes
. numprocs - number of indices for each process
- indices  - indices (in local numbering) shared with neighbors (sorted by global numbering)

  Level: advanced

.seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
          `ISLocalToGlobalMappingGetInfo()`
@*/
PetscErrorCode ISLocalToGlobalMappingRestoreInfo(ISLocalToGlobalMapping mapping, PetscInt *nproc, PetscInt *procs[], PetscInt *numprocs[], PetscInt **indices[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mapping, IS_LTOGM_CLASSID, 1);
  if (mapping->bs > 1) {
    if (numprocs) PetscCall(PetscFree(*numprocs));
    if (indices) {
      if (*indices)
        for (PetscInt i = 0; i < *nproc; i++) PetscCall(PetscFree((*indices)[i]));
      PetscCall(PetscFree(*indices));
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  ISLocalToGlobalMappingGetNodeInfo - Gets the neighbor information of local nodes

  Collective the first time it is called

  Input Parameter:
. mapping - the mapping from local to global indexing

  Output Parameters:
+ n       - number of local nodes
. n_procs - an array storing the number of processes for each local node (including self)
- procs   - the processes' rank for each local node (sorted, self is first)

  Level: advanced

  Note:
  The user needs to call `ISLocalToGlobalMappingRestoreNodeInfo()` when the data is no longer needed.

.seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
          `ISLocalToGlobalMappingGetInfo()`, `ISLocalToGlobalMappingRestoreNodeInfo()`, `ISLocalToGlobalMappingGetBlockNodeInfo()`
@*/
PetscErrorCode ISLocalToGlobalMappingGetNodeInfo(ISLocalToGlobalMapping mapping, PetscInt *n, PetscInt *n_procs[], PetscInt **procs[])
{
  PetscInt **bprocs = NULL, *bn_procs = NULL, bs, i, j, k, bn;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mapping, IS_LTOGM_CLASSID, 1);
  bs = mapping->bs;
  PetscCall(ISLocalToGlobalMappingGetBlockNodeInfo(mapping, &bn, &bn_procs, &bprocs));
  if (bs > 1) { /* we need to expand the cached info */
    PetscInt *tn_procs;
    PetscInt  c;

    PetscCall(PetscMalloc1(bn * bs, &tn_procs));
    for (i = 0, c = 0; i < bn; i++) {
      for (k = 0; k < bs; k++) tn_procs[i * bs + k] = bn_procs[i];
      c += bs * bn_procs[i];
    }
    if (n) *n = bn * bs;
    if (procs) {
      PetscInt **tprocs;
      PetscInt   tn = bn * bs;

      PetscCall(PetscMalloc1(tn, &tprocs));
      if (tn) PetscCall(PetscMalloc1(c, &tprocs[0]));
      for (i = 0; i < tn - 1; i++) tprocs[i + 1] = tprocs[i] + tn_procs[i];
      for (i = 0; i < bn; i++) {
        for (k = 0; k < bs; k++) {
          for (j = 0; j < bn_procs[i]; j++) tprocs[i * bs + k][j] = bprocs[i][j];
        }
      }
      *procs = tprocs;
    }
    if (n_procs) *n_procs = tn_procs;
    else PetscCall(PetscFree(tn_procs));
  } else {
    if (n) *n = bn;
    if (n_procs) *n_procs = bn_procs;
    if (procs) *procs = bprocs;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  ISLocalToGlobalMappingRestoreNodeInfo - Frees the memory allocated by `ISLocalToGlobalMappingGetNodeInfo()`

  Not Collective

  Input Parameters:
+ mapping - the mapping from local to global indexing
. n       - number of local nodes
. n_procs - an array storing the number of processes for each local node (including self)
- procs   - the processes' rank for each local node (sorted, self is first)

  Level: advanced

.seealso: [](sec_scatter), `ISLocalToGlobalMappingDestroy()`, `ISLocalToGlobalMappingCreateIS()`, `ISLocalToGlobalMappingCreate()`,
          `ISLocalToGlobalMappingGetInfo()`
@*/
PetscErrorCode ISLocalToGlobalMappingRestoreNodeInfo(ISLocalToGlobalMapping mapping, PetscInt *n, PetscInt *n_procs[], PetscInt **procs[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(mapping, IS_LTOGM_CLASSID, 1);
  if (mapping->bs > 1) {
    if (n_procs) PetscCall(PetscFree(*n_procs));
    if (procs) {
      if (*procs) PetscCall(PetscFree((*procs)[0]));
      PetscCall(PetscFree(*procs));
    }
  }
  PetscCall(ISLocalToGlobalMappingRestoreBlockNodeInfo(mapping, n, n_procs, procs));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*MC
    ISLocalToGlobalMappingGetIndicesF90 - Get global indices for every local point that is mapped

    Synopsis:
    ISLocalToGlobalMappingGetIndicesF90(ISLocalToGlobalMapping ltog, PetscInt, pointer :: array(:)}, integer ierr)

    Not Collective

    Input Parameter:
.   A - the matrix

    Output Parameter:
.   array - array of indices, the length of this array may be obtained with `ISLocalToGlobalMappingGetSize()`

    Level: advanced

    Note:
    Use  `ISLocalToGlobalMappingGetIndicesF90()` when you no longer need access to the data

.seealso: [](sec_fortranarrays), [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingGetIndices()`, `ISLocalToGlobalMappingRestoreIndices()`,
          `ISLocalToGlobalMappingRestoreIndicesF90()`
M*/

/*MC
    ISLocalToGlobalMappingRestoreIndicesF90 - restores the global indices for every local point that is mapped obtained with `ISLocalToGlobalMappingGetIndicesF90()`

    Synopsis:
    ISLocalToGlobalMappingRestoreIndicesF90(ISLocalToGlobalMapping ltog, PetscInt, pointer :: array(:)}, integer ierr)

    Not Collective

    Input Parameters:
+   A - the matrix
-   array - array of indices, the length of this array may be obtained with `ISLocalToGlobalMappingGetSize()`

    Level: advanced

.seealso: [](sec_fortranarrays), [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingGetIndices()`, `ISLocalToGlobalMappingRestoreIndices()`,
          `ISLocalToGlobalMappingGetIndicesF90()`
M*/

/*@C
  ISLocalToGlobalMappingGetIndices - Get global indices for every local point that is mapped

  Not Collective

  Input Parameter:
. ltog - local to global mapping

  Output Parameter:
. array - array of indices, the length of this array may be obtained with `ISLocalToGlobalMappingGetSize()`

  Level: advanced

  Note:
  `ISLocalToGlobalMappingGetSize()` returns the length the this array

.seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingApply()`, `ISLocalToGlobalMappingRestoreIndices()`,
          `ISLocalToGlobalMappingGetBlockIndices()`, `ISLocalToGlobalMappingRestoreBlockIndices()`
@*/
PetscErrorCode ISLocalToGlobalMappingGetIndices(ISLocalToGlobalMapping ltog, const PetscInt **array)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ltog, IS_LTOGM_CLASSID, 1);
  PetscAssertPointer(array, 2);
  if (ltog->bs == 1) {
    *array = ltog->indices;
  } else {
    PetscInt       *jj, k, i, j, n = ltog->n, bs = ltog->bs;
    const PetscInt *ii;

    PetscCall(PetscMalloc1(bs * n, &jj));
    *array = jj;
    k      = 0;
    ii     = ltog->indices;
    for (i = 0; i < n; i++)
      for (j = 0; j < bs; j++) jj[k++] = bs * ii[i] + j;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  ISLocalToGlobalMappingRestoreIndices - Restore indices obtained with `ISLocalToGlobalMappingGetIndices()`

  Not Collective

  Input Parameters:
+ ltog  - local to global mapping
- array - array of indices

  Level: advanced

.seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingApply()`, `ISLocalToGlobalMappingGetIndices()`
@*/
PetscErrorCode ISLocalToGlobalMappingRestoreIndices(ISLocalToGlobalMapping ltog, const PetscInt **array)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ltog, IS_LTOGM_CLASSID, 1);
  PetscAssertPointer(array, 2);
  PetscCheck(ltog->bs != 1 || *array == ltog->indices, PETSC_COMM_SELF, PETSC_ERR_ARG_BADPTR, "Trying to return mismatched pointer");
  if (ltog->bs > 1) PetscCall(PetscFree(*(void **)array));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  ISLocalToGlobalMappingGetBlockIndices - Get global indices for every local block

  Not Collective

  Input Parameter:
. ltog - local to global mapping

  Output Parameter:
. array - array of indices

  Level: advanced

.seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingApply()`, `ISLocalToGlobalMappingRestoreBlockIndices()`
@*/
PetscErrorCode ISLocalToGlobalMappingGetBlockIndices(ISLocalToGlobalMapping ltog, const PetscInt **array)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ltog, IS_LTOGM_CLASSID, 1);
  PetscAssertPointer(array, 2);
  *array = ltog->indices;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  ISLocalToGlobalMappingRestoreBlockIndices - Restore indices obtained with `ISLocalToGlobalMappingGetBlockIndices()`

  Not Collective

  Input Parameters:
+ ltog  - local to global mapping
- array - array of indices

  Level: advanced

.seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingApply()`, `ISLocalToGlobalMappingGetIndices()`
@*/
PetscErrorCode ISLocalToGlobalMappingRestoreBlockIndices(ISLocalToGlobalMapping ltog, const PetscInt **array)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ltog, IS_LTOGM_CLASSID, 1);
  PetscAssertPointer(array, 2);
  PetscCheck(*array == ltog->indices, PETSC_COMM_SELF, PETSC_ERR_ARG_BADPTR, "Trying to return mismatched pointer");
  *array = NULL;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  ISLocalToGlobalMappingConcatenate - Create a new mapping that concatenates a list of mappings

  Not Collective

  Input Parameters:
+ comm  - communicator for the new mapping, must contain the communicator of every mapping to concatenate
. n     - number of mappings to concatenate
- ltogs - local to global mappings

  Output Parameter:
. ltogcat - new mapping

  Level: advanced

  Note:
  This currently always returns a mapping with block size of 1

  Developer Notes:
  If all the input mapping have the same block size we could easily handle that as a special case

.seealso: [](sec_scatter), `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingCreate()`
@*/
PetscErrorCode ISLocalToGlobalMappingConcatenate(MPI_Comm comm, PetscInt n, const ISLocalToGlobalMapping ltogs[], ISLocalToGlobalMapping *ltogcat)
{
  PetscInt i, cnt, m, *idx;

  PetscFunctionBegin;
  PetscCheck(n >= 0, comm, PETSC_ERR_ARG_OUTOFRANGE, "Must have a non-negative number of mappings, given %" PetscInt_FMT, n);
  if (n > 0) PetscAssertPointer(ltogs, 3);
  for (i = 0; i < n; i++) PetscValidHeaderSpecific(ltogs[i], IS_LTOGM_CLASSID, 3);
  PetscAssertPointer(ltogcat, 4);
  for (cnt = 0, i = 0; i < n; i++) {
    PetscCall(ISLocalToGlobalMappingGetSize(ltogs[i], &m));
    cnt += m;
  }
  PetscCall(PetscMalloc1(cnt, &idx));
  for (cnt = 0, i = 0; i < n; i++) {
    const PetscInt *subidx;
    PetscCall(ISLocalToGlobalMappingGetSize(ltogs[i], &m));
    PetscCall(ISLocalToGlobalMappingGetIndices(ltogs[i], &subidx));
    PetscCall(PetscArraycpy(&idx[cnt], subidx, m));
    PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogs[i], &subidx));
    cnt += m;
  }
  PetscCall(ISLocalToGlobalMappingCreate(comm, 1, cnt, idx, PETSC_OWN_POINTER, ltogcat));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*MC
      ISLOCALTOGLOBALMAPPINGBASIC - basic implementation of the `ISLocalToGlobalMapping` object. When `ISGlobalToLocalMappingApply()` is
                                    used this is good for only small and moderate size problems.

   Options Database Key:
.   -islocaltoglobalmapping_type basic - select this method

   Level: beginner

   Developer Note:
   This stores all the mapping information on each MPI rank.

.seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingSetType()`, `ISLOCALTOGLOBALMAPPINGHASH`
M*/
PETSC_EXTERN PetscErrorCode ISLocalToGlobalMappingCreate_Basic(ISLocalToGlobalMapping ltog)
{
  PetscFunctionBegin;
  ltog->ops->globaltolocalmappingapply      = ISGlobalToLocalMappingApply_Basic;
  ltog->ops->globaltolocalmappingsetup      = ISGlobalToLocalMappingSetUp_Basic;
  ltog->ops->globaltolocalmappingapplyblock = ISGlobalToLocalMappingApplyBlock_Basic;
  ltog->ops->destroy                        = ISLocalToGlobalMappingDestroy_Basic;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*MC
      ISLOCALTOGLOBALMAPPINGHASH - hash implementation of the `ISLocalToGlobalMapping` object. When `ISGlobalToLocalMappingApply()` is
                                    used this is good for large memory problems.

   Options Database Key:
.   -islocaltoglobalmapping_type hash - select this method

   Level: beginner

   Note:
    This is selected automatically for large problems if the user does not set the type.

.seealso: [](sec_scatter), `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingSetType()`, `ISLOCALTOGLOBALMAPPINGBASIC`
M*/
PETSC_EXTERN PetscErrorCode ISLocalToGlobalMappingCreate_Hash(ISLocalToGlobalMapping ltog)
{
  PetscFunctionBegin;
  ltog->ops->globaltolocalmappingapply      = ISGlobalToLocalMappingApply_Hash;
  ltog->ops->globaltolocalmappingsetup      = ISGlobalToLocalMappingSetUp_Hash;
  ltog->ops->globaltolocalmappingapplyblock = ISGlobalToLocalMappingApplyBlock_Hash;
  ltog->ops->destroy                        = ISLocalToGlobalMappingDestroy_Hash;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  ISLocalToGlobalMappingRegister -  Registers a method for applying a global to local mapping with an `ISLocalToGlobalMapping`

  Not Collective

  Input Parameters:
+ sname    - name of a new method
- function - routine to create method context

  Example Usage:
.vb
   ISLocalToGlobalMappingRegister("my_mapper", MyCreate);
.ve

  Then, your mapping can be chosen with the procedural interface via
$     ISLocalToGlobalMappingSetType(ltog, "my_mapper")
  or at runtime via the option
$     -islocaltoglobalmapping_type my_mapper

  Level: advanced

  Note:
  `ISLocalToGlobalMappingRegister()` may be called multiple times to add several user-defined mappings.

.seealso: [](sec_scatter), `ISLocalToGlobalMappingRegisterAll()`, `ISLocalToGlobalMappingRegisterDestroy()`, `ISLOCALTOGLOBALMAPPINGBASIC`,
          `ISLOCALTOGLOBALMAPPINGHASH`, `ISLocalToGlobalMapping`, `ISLocalToGlobalMappingApply()`
@*/
PetscErrorCode ISLocalToGlobalMappingRegister(const char sname[], PetscErrorCode (*function)(ISLocalToGlobalMapping))
{
  PetscFunctionBegin;
  PetscCall(ISInitializePackage());
  PetscCall(PetscFunctionListAdd(&ISLocalToGlobalMappingList, sname, function));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  ISLocalToGlobalMappingSetType - Sets the implementation type `ISLocalToGlobalMapping` will use

  Logically Collective

  Input Parameters:
+ ltog - the `ISLocalToGlobalMapping` object
- type - a known method

  Options Database Key:
. -islocaltoglobalmapping_type  <method> - Sets the method; use -help for a list of available methods (for instance, basic or hash)

  Level: intermediate

  Notes:
  See `ISLocalToGlobalMappingType` for available methods

  Normally, it is best to use the `ISLocalToGlobalMappingSetFromOptions()` command and
  then set the `ISLocalToGlobalMappingType` from the options database rather than by using
  this routine.

  Developer Notes:
  `ISLocalToGlobalMappingRegister()` is used to add new types to `ISLocalToGlobalMappingList` from which they
  are accessed by `ISLocalToGlobalMappingSetType()`.

.seealso: [](sec_scatter), `ISLocalToGlobalMappingType`, `ISLocalToGlobalMappingRegister()`, `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingGetType()`
@*/
PetscErrorCode ISLocalToGlobalMappingSetType(ISLocalToGlobalMapping ltog, ISLocalToGlobalMappingType type)
{
  PetscBool match;
  PetscErrorCode (*r)(ISLocalToGlobalMapping) = NULL;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ltog, IS_LTOGM_CLASSID, 1);
  if (type) PetscAssertPointer(type, 2);

  PetscCall(PetscObjectTypeCompare((PetscObject)ltog, type, &match));
  if (match) PetscFunctionReturn(PETSC_SUCCESS);

  /* L2G maps defer type setup at globaltolocal calls, allow passing NULL here */
  if (type) {
    PetscCall(PetscFunctionListFind(ISLocalToGlobalMappingList, type, &r));
    PetscCheck(r, PetscObjectComm((PetscObject)ltog), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unable to find requested ISLocalToGlobalMapping type %s", type);
  }
  /* Destroy the previous private LTOG context */
  PetscTryTypeMethod(ltog, destroy);
  ltog->ops->destroy = NULL;

  PetscCall(PetscObjectChangeTypeName((PetscObject)ltog, type));
  if (r) PetscCall((*r)(ltog));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  ISLocalToGlobalMappingGetType - Get the type of the `ISLocalToGlobalMapping`

  Not Collective

  Input Parameter:
. ltog - the `ISLocalToGlobalMapping` object

  Output Parameter:
. type - the type

  Level: intermediate

.seealso: [](sec_scatter), `ISLocalToGlobalMappingType`, `ISLocalToGlobalMappingRegister()`, `ISLocalToGlobalMappingCreate()`, `ISLocalToGlobalMappingSetType()`
@*/
PetscErrorCode ISLocalToGlobalMappingGetType(ISLocalToGlobalMapping ltog, ISLocalToGlobalMappingType *type)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ltog, IS_LTOGM_CLASSID, 1);
  PetscAssertPointer(type, 2);
  *type = ((PetscObject)ltog)->type_name;
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscBool ISLocalToGlobalMappingRegisterAllCalled = PETSC_FALSE;

/*@C
  ISLocalToGlobalMappingRegisterAll - Registers all of the local to global mapping components in the `IS` package.

  Not Collective

  Level: advanced

.seealso: [](sec_scatter), `ISRegister()`, `ISLocalToGlobalRegister()`
@*/
PetscErrorCode ISLocalToGlobalMappingRegisterAll(void)
{
  PetscFunctionBegin;
  if (ISLocalToGlobalMappingRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
  ISLocalToGlobalMappingRegisterAllCalled = PETSC_TRUE;
  PetscCall(ISLocalToGlobalMappingRegister(ISLOCALTOGLOBALMAPPINGBASIC, ISLocalToGlobalMappingCreate_Basic));
  PetscCall(ISLocalToGlobalMappingRegister(ISLOCALTOGLOBALMAPPINGHASH, ISLocalToGlobalMappingCreate_Hash));
  PetscFunctionReturn(PETSC_SUCCESS);
}