xref: /petsc/src/vec/is/is/utils/isdiff.c (revision 5a884c48ab0c46bab83cd9bb8710f380fa6d8bcf)

#include <petsc/private/isimpl.h> /*I "petscis.h"  I*/
#include <petsc/private/sectionimpl.h>
#include <petscbt.h>

/*@
  ISDifference - Computes the difference between two index sets.

  Collective

  Input Parameters:
+ is1 - first index, to have items removed from it
- is2 - index values to be removed

  Output Parameter:
. isout - is1 - is2

  Level: intermediate

  Notes:
  Negative values are removed from the lists. `is2` may have values
  that are not in `is1`.

  This computation requires O(imax-imin) memory and O(imax-imin)
  work, where imin and imax are the bounds on the indices in is1.

  If `is2` is `NULL`, the result is the same as for an empty `IS`, i.e., a duplicate of `is1`.

  The difference is computed separately on each MPI rank

.seealso: [](sec_scatter), `IS`, `ISDestroy()`, `ISView()`, `ISSum()`, `ISExpand()`
@*/
PetscErrorCode ISDifference(IS is1, IS is2, IS *isout)
{
  PetscInt        i, n1, n2, imin, imax, nout, *iout;
  const PetscInt *i1, *i2;
  PetscBT         mask;
  MPI_Comm        comm;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(is1, IS_CLASSID, 1);
  PetscAssertPointer(isout, 3);
  if (!is2) {
    PetscCall(ISDuplicate(is1, isout));
    PetscFunctionReturn(PETSC_SUCCESS);
  }
  PetscValidHeaderSpecific(is2, IS_CLASSID, 2);

  PetscCall(ISGetIndices(is1, &i1));
  PetscCall(ISGetLocalSize(is1, &n1));

  /* Create a bit mask array to contain required values */
  if (n1) {
    imin = PETSC_INT_MAX;
    imax = 0;
    for (i = 0; i < n1; i++) {
      if (i1[i] < 0) continue;
      imin = PetscMin(imin, i1[i]);
      imax = PetscMax(imax, i1[i]);
    }
  } else imin = imax = 0;

  PetscCall(PetscBTCreate(imax - imin, &mask));
  /* Put the values from is1 */
  for (i = 0; i < n1; i++) {
    if (i1[i] < 0) continue;
    PetscCall(PetscBTSet(mask, i1[i] - imin));
  }
  PetscCall(ISRestoreIndices(is1, &i1));
  /* Remove the values from is2 */
  PetscCall(ISGetIndices(is2, &i2));
  PetscCall(ISGetLocalSize(is2, &n2));
  for (i = 0; i < n2; i++) {
    if (i2[i] < imin || i2[i] > imax) continue;
    PetscCall(PetscBTClear(mask, i2[i] - imin));
  }
  PetscCall(ISRestoreIndices(is2, &i2));

  /* Count the number in the difference */
  nout = 0;
  for (i = 0; i < imax - imin + 1; i++) {
    if (PetscBTLookup(mask, i)) nout++;
  }

  /* create the new IS containing the difference */
  PetscCall(PetscMalloc1(nout, &iout));
  nout = 0;
  for (i = 0; i < imax - imin + 1; i++) {
    if (PetscBTLookup(mask, i)) iout[nout++] = i + imin;
  }
  PetscCall(PetscObjectGetComm((PetscObject)is1, &comm));
  PetscCall(ISCreateGeneral(comm, nout, iout, PETSC_OWN_POINTER, isout));

  PetscCall(PetscBTDestroy(&mask));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  ISSum - Computes the sum (union) of two index sets.

  Only sequential version (at the moment)

  Input Parameters:
+ is1 - index set to be extended
- is2 - index values to be added

  Output Parameter:
. is3 - the sum; this can not be `is1` or `is2`

  Level: intermediate

  Notes:
  If n1 and n2 are the sizes of the sets, this takes O(n1+n2) time;

  Both index sets need to be sorted on input.

  The sum is computed separately on each MPI rank

.seealso: [](sec_scatter), `IS`, `ISDestroy()`, `ISView()`, `ISDifference()`, `ISExpand()`
@*/
PetscErrorCode ISSum(IS is1, IS is2, IS *is3)
{
  PetscBool       f;
  const PetscInt *i1, *i2;
  PetscInt        n1, n2, n3, p1, p2, *iout;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(is1, IS_CLASSID, 1);
  PetscValidHeaderSpecific(is2, IS_CLASSID, 2);
  PetscCheckSameComm(is1, 1, is2, 2);

  PetscCall(ISSorted(is1, &f));
  PetscCheck(f, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Arg 1 is not sorted");
  PetscCall(ISSorted(is2, &f));
  PetscCheck(f, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Arg 2 is not sorted");

  PetscCall(ISGetLocalSize(is1, &n1));
  PetscCall(ISGetLocalSize(is2, &n2));
  if (!n2) {
    PetscCall(ISDuplicate(is1, is3));
    PetscFunctionReturn(PETSC_SUCCESS);
  }
  PetscCall(ISGetIndices(is1, &i1));
  PetscCall(ISGetIndices(is2, &i2));

  p1 = 0;
  p2 = 0;
  n3 = 0;
  do {
    if (p1 == n1) { /* cleanup for is2 */
      n3 += n2 - p2;
      break;
    } else {
      while (p2 < n2 && i2[p2] < i1[p1]) {
        n3++;
        p2++;
      }
      if (p2 == n2) {
        /* cleanup for is1 */
        n3 += n1 - p1;
        break;
      } else {
        if (i2[p2] == i1[p1]) {
          n3++;
          p1++;
          p2++;
        }
      }
    }
    if (p2 == n2) {
      /* cleanup for is1 */
      n3 += n1 - p1;
      break;
    } else {
      while (p1 < n1 && i1[p1] < i2[p2]) {
        n3++;
        p1++;
      }
      if (p1 == n1) {
        /* clean up for is2 */
        n3 += n2 - p2;
        break;
      } else {
        if (i1[p1] == i2[p2]) {
          n3++;
          p1++;
          p2++;
        }
      }
    }
  } while (p1 < n1 || p2 < n2);

  if (n3 == n1) { /* no new elements to be added */
    PetscCall(ISRestoreIndices(is1, &i1));
    PetscCall(ISRestoreIndices(is2, &i2));
    PetscCall(ISDuplicate(is1, is3));
    PetscFunctionReturn(PETSC_SUCCESS);
  }
  PetscCall(PetscMalloc1(n3, &iout));

  p1 = 0;
  p2 = 0;
  n3 = 0;
  do {
    if (p1 == n1) { /* cleanup for is2 */
      while (p2 < n2) iout[n3++] = i2[p2++];
      break;
    } else {
      while (p2 < n2 && i2[p2] < i1[p1]) iout[n3++] = i2[p2++];
      if (p2 == n2) { /* cleanup for is1 */
        while (p1 < n1) iout[n3++] = i1[p1++];
        break;
      } else {
        if (i2[p2] == i1[p1]) {
          iout[n3++] = i1[p1++];
          p2++;
        }
      }
    }
    if (p2 == n2) { /* cleanup for is1 */
      while (p1 < n1) iout[n3++] = i1[p1++];
      break;
    } else {
      while (p1 < n1 && i1[p1] < i2[p2]) iout[n3++] = i1[p1++];
      if (p1 == n1) { /* clean up for is2 */
        while (p2 < n2) iout[n3++] = i2[p2++];
        break;
      } else {
        if (i1[p1] == i2[p2]) {
          iout[n3++] = i1[p1++];
          p2++;
        }
      }
    }
  } while (p1 < n1 || p2 < n2);

  PetscCall(ISRestoreIndices(is1, &i1));
  PetscCall(ISRestoreIndices(is2, &i2));
  PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is1), n3, iout, PETSC_OWN_POINTER, is3));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  ISExpand - Computes the union of two index sets, by concatenating 2 lists and
  removing duplicates.

  Collective

  Input Parameters:
+ is1 - first index set
- is2 - index values to be added

  Output Parameter:
. isout - `is1` + `is2` The index set `is2` is appended to `is1` removing duplicates

  Level: intermediate

  Notes:
  Negative values are removed from the lists. This requires O(imax-imin)
  memory and O(imax-imin) work, where imin and imax are the bounds on the
  indices in `is1` and `is2`.

  `is1` and `is2` do not need to be sorted.

  The operations are performed separately on each MPI rank

.seealso: [](sec_scatter), `IS`, `ISDestroy()`, `ISView()`, `ISDifference()`, `ISSum()`, `ISIntersect()`
@*/
PetscErrorCode ISExpand(IS is1, IS is2, IS *isout)
{
  PetscInt        i, n1, n2, imin, imax, nout, *iout;
  const PetscInt *i1, *i2;
  PetscBool       sorted1 = PETSC_TRUE, sorted2 = PETSC_TRUE;
  PetscBT         mask;
  MPI_Comm        comm;

  PetscFunctionBegin;
  if (is1) PetscValidHeaderSpecific(is1, IS_CLASSID, 1);
  if (is2) PetscValidHeaderSpecific(is2, IS_CLASSID, 2);
  PetscAssertPointer(isout, 3);

  PetscCheck(is1 || is2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Both arguments cannot be NULL");
  if (!is1) {
    PetscCall(ISDuplicate(is2, isout));
    PetscFunctionReturn(PETSC_SUCCESS);
  }
  if (!is2) {
    PetscCall(ISDuplicate(is1, isout));
    PetscFunctionReturn(PETSC_SUCCESS);
  }
  PetscCall(ISGetIndices(is1, &i1));
  PetscCall(ISGetLocalSize(is1, &n1));
  PetscCall(ISGetIndices(is2, &i2));
  PetscCall(ISGetLocalSize(is2, &n2));

  /* Create a bit mask array to contain required values */
  if (n1 || n2) {
    imin = PETSC_INT_MAX;
    imax = 0;
    if (n1) {
      PetscCall(ISSorted(is1, &sorted1));
      if (sorted1 && i1[0] >= 0) imin = i1[0], imax = i1[n1 - 1];
      else
        for (i = 0; i < n1; i++) {
          if (i1[i] < 0) continue;
          imin = PetscMin(imin, i1[i]);
          imax = PetscMax(imax, i1[i]);
        }
    }
    if (n2) {
      PetscCall(ISSorted(is2, &sorted2));
      if (sorted2 && i2[0] >= 0) imin = PetscMin(imin, i2[0]), imax = PetscMax(imax, i2[n2 - 1]);
      else
        for (i = 0; i < n2; i++) {
          if (i2[i] < 0) continue;
          imin = PetscMin(imin, i2[i]);
          imax = PetscMax(imax, i2[i]);
        }
    }
  } else imin = imax = 0;

  PetscCall(PetscMalloc1(n1 + n2, &iout));
  nout = 0;
  PetscCall(PetscBTCreate(imax - imin, &mask));
  /* Put the values from is1 */
  for (i = 0; i < n1; i++) {
    if (i1[i] < 0) continue;
    if (!PetscBTLookupSet(mask, i1[i] - imin)) iout[nout++] = i1[i];
  }
  n1 = -nout;
  PetscCall(ISRestoreIndices(is1, &i1));
  /* Put the values from is2 */
  for (i = 0; i < n2; i++) {
    if (i2[i] < 0) continue;
    if (!PetscBTLookupSet(mask, i2[i] - imin)) iout[nout++] = i2[i];
  }
  PetscCall(ISRestoreIndices(is2, &i2));

  /* create the new IS containing the sum */
  PetscCall(PetscObjectGetComm((PetscObject)is1, &comm));
  PetscCall(ISCreateGeneral(comm, nout, iout, PETSC_OWN_POINTER, isout));
  /* no entries of is2 (resp. is1) was inserted, so if is1 (resp. is2) is sorted, then so is isout */
  if ((-n1 == nout && sorted1) || (n1 == 0 && sorted2)) PetscCall(ISSetInfo(*isout, IS_SORTED, IS_LOCAL, PETSC_FALSE, PETSC_TRUE));
  PetscCall(PetscBTDestroy(&mask));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  ISIntersect - Computes the intersection of two index sets, by sorting and comparing.

  Collective

  Input Parameters:
+ is1 - first index set
- is2 - second index set

  Output Parameter:
. isout - the sorted intersection of `is1` and `is2`

  Level: intermediate

  Notes:
  Negative values are removed from the lists. This requires O(min(is1,is2))
  memory and O(max(is1,is2)log(max(is1,is2))) work

  `is1` and `is2` do not need to be sorted.

  The operations are performed separately on each MPI rank

.seealso: [](sec_scatter), `IS`, `ISDestroy()`, `ISView()`, `ISDifference()`, `ISSum()`, `ISExpand()`, `ISConcatenate()`
@*/
PetscErrorCode ISIntersect(IS is1, IS is2, IS *isout)
{
  PetscInt        i, n1, n2, nout, *iout;
  const PetscInt *i1, *i2;
  IS              is1sorted = NULL, is2sorted = NULL;
  PetscBool       sorted, lsorted;
  MPI_Comm        comm;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(is1, IS_CLASSID, 1);
  PetscValidHeaderSpecific(is2, IS_CLASSID, 2);
  PetscCheckSameComm(is1, 1, is2, 2);
  PetscAssertPointer(isout, 3);
  PetscCall(PetscObjectGetComm((PetscObject)is1, &comm));

  PetscCall(ISGetLocalSize(is1, &n1));
  PetscCall(ISGetLocalSize(is2, &n2));
  if (n1 < n2) {
    IS       tempis = is1;
    PetscInt ntemp  = n1;

    is1 = is2;
    is2 = tempis;
    n1  = n2;
    n2  = ntemp;
  }
  PetscCall(ISSorted(is1, &lsorted));
  PetscCallMPI(MPIU_Allreduce(&lsorted, &sorted, 1, MPI_C_BOOL, MPI_LAND, comm));
  if (!sorted) {
    PetscCall(ISDuplicate(is1, &is1sorted));
    PetscCall(ISSort(is1sorted));
    PetscCall(ISGetIndices(is1sorted, &i1));
  } else {
    is1sorted = is1;
    PetscCall(PetscObjectReference((PetscObject)is1));
    PetscCall(ISGetIndices(is1, &i1));
  }
  PetscCall(ISSorted(is2, &lsorted));
  PetscCallMPI(MPIU_Allreduce(&lsorted, &sorted, 1, MPI_C_BOOL, MPI_LAND, comm));
  if (!sorted) {
    PetscCall(ISDuplicate(is2, &is2sorted));
    PetscCall(ISSort(is2sorted));
    PetscCall(ISGetIndices(is2sorted, &i2));
  } else {
    is2sorted = is2;
    PetscCall(PetscObjectReference((PetscObject)is2));
    PetscCall(ISGetIndices(is2, &i2));
  }

  PetscCall(PetscMalloc1(n2, &iout));

  for (i = 0, nout = 0; i < n2; i++) {
    PetscInt key = i2[i];
    PetscInt loc;

    PetscCall(ISLocate(is1sorted, key, &loc));
    if (loc >= 0) {
      if (!nout || iout[nout - 1] < key) iout[nout++] = key;
    }
  }
  PetscCall(PetscRealloc(nout * sizeof(PetscInt), &iout));

  /* create the new IS containing the sum */
  PetscCall(ISCreateGeneral(comm, nout, iout, PETSC_OWN_POINTER, isout));
  PetscCall(ISSetInfo(*isout, IS_SORTED, IS_GLOBAL, PETSC_FALSE, PETSC_TRUE));
  PetscCall(ISRestoreIndices(is2sorted, &i2));
  PetscCall(ISDestroy(&is2sorted));
  PetscCall(ISRestoreIndices(is1sorted, &i1));
  PetscCall(ISDestroy(&is1sorted));
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode ISIntersect_Caching_Internal(IS is1, IS is2, IS *isect)
{
  PetscFunctionBegin;
  *isect = NULL;
  if (is2 && is1) {
    char          composeStr[33] = {0};
    PetscObjectId is2id;

    PetscCall(PetscObjectGetId((PetscObject)is2, &is2id));
    PetscCall(PetscSNPrintf(composeStr, 32, "ISIntersect_Caching_%" PetscInt64_FMT, is2id));
    PetscCall(PetscObjectQuery((PetscObject)is1, composeStr, (PetscObject *)isect));
    if (*isect == NULL) {
      PetscCall(ISIntersect(is1, is2, isect));
      PetscCall(PetscObjectCompose((PetscObject)is1, composeStr, (PetscObject)*isect));
    } else {
      PetscCall(PetscObjectReference((PetscObject)*isect));
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  ISConcatenate - Forms a new `IS` by locally concatenating the indices from an `IS` list without reordering.

  Collective

  Input Parameters:
+ comm   - communicator of the concatenated `IS`.
. len    - size of islist array (nonnegative)
- islist - array of index sets

  Output Parameter:
. isout - The concatenated index set; empty, if `len` == 0.

  Level: intermediate

  Notes:
  The semantics of calling this on comm imply that the comms of the members of `islist` also contain this rank.

.seealso: [](sec_scatter), `IS`, `ISDifference()`, `ISSum()`, `ISExpand()`, `ISIntersect()`
@*/
PetscErrorCode ISConcatenate(MPI_Comm comm, PetscInt len, const IS islist[], IS *isout)
{
  PetscInt        i, n, N;
  const PetscInt *iidx;
  PetscInt       *idx;

  PetscFunctionBegin;
  if (len) PetscAssertPointer(islist, 3);
  if (PetscDefined(USE_DEBUG)) {
    for (i = 0; i < len; ++i)
      if (islist[i]) PetscValidHeaderSpecific(islist[i], IS_CLASSID, 3);
  }
  PetscAssertPointer(isout, 4);
  if (!len) {
    PetscCall(ISCreateGeneral(comm, 0, NULL, PETSC_OWN_POINTER, isout));
    PetscFunctionReturn(PETSC_SUCCESS);
  }
  PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Negative array length: %" PetscInt_FMT, len);
  N = 0;
  for (i = 0; i < len; ++i) {
    if (islist[i]) {
      PetscCall(ISGetLocalSize(islist[i], &n));
      N += n;
    }
  }
  PetscCall(PetscMalloc1(N, &idx));
  N = 0;
  for (i = 0; i < len; ++i) {
    if (islist[i]) {
      PetscCall(ISGetLocalSize(islist[i], &n));
      if (n) {
        PetscCall(ISGetIndices(islist[i], &iidx));
        PetscCall(PetscArraycpy(idx + N, iidx, n));
        PetscCall(ISRestoreIndices(islist[i], &iidx));
        N += n;
      }
    }
  }
  PetscCall(ISCreateGeneral(comm, N, idx, PETSC_OWN_POINTER, isout));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  ISListToPair  - Convert an `IS` list to a pair of `IS` of equal length defining an equivalent integer multimap.
  Each `IS` in `islist` is assigned an integer j so that all of the indices of that `IS` are
  mapped to j.

  Collective

  Input Parameters:
+ comm    - `MPI_Comm`
. listlen - `IS` list length
- islist  - `IS` list

  Output Parameters:
+ xis - domain `IS`
- yis - range  `IS`

  Level: developer

  Notes:
  The global integers assigned to the `IS` of `islist` might not correspond to the
  local numbers of the `IS` on that list, but the two *orderings* are the same: the global
  integers assigned to the `IS` on the local list form a strictly increasing sequence.

  The `IS` in `islist` can belong to subcommunicators of `comm`, and the subcommunicators
  on the input `IS` list are assumed to be in a "deadlock-free" order.

  Local lists of `PetscObject`s (or their subcomms) on a comm are "deadlock-free" if subcomm1
  precedes subcomm2 on any local list, then it precedes subcomm2 on all ranks.
  Equivalently, the local numbers of the subcomms on each local list are drawn from some global
  numbering. This is ensured, for example, by `ISPairToList()`.

.seealso: `IS`, `ISPairToList()`
@*/
PetscErrorCode ISListToPair(MPI_Comm comm, PetscInt listlen, IS islist[], IS *xis, IS *yis)
{
  PetscInt        ncolors, *colors, leni, len, *xinds, *yinds, k;
  const PetscInt *indsi;

  PetscFunctionBegin;
  PetscCall(PetscMalloc1(listlen, &colors));
  PetscCall(PetscObjectsListGetGlobalNumbering(comm, listlen, (PetscObject *)islist, &ncolors, colors));
  len = 0;
  for (PetscInt i = 0; i < listlen; ++i) {
    PetscCall(ISGetLocalSize(islist[i], &leni));
    len += leni;
  }
  PetscCall(PetscMalloc1(len, &xinds));
  PetscCall(PetscMalloc1(len, &yinds));
  k = 0;
  for (PetscInt i = 0; i < listlen; ++i) {
    PetscCall(ISGetLocalSize(islist[i], &leni));
    PetscCall(ISGetIndices(islist[i], &indsi));
    for (PetscInt j = 0; j < leni; ++j) {
      xinds[k] = indsi[j];
      yinds[k] = colors[i];
      ++k;
    }
  }
  PetscCall(PetscFree(colors));
  PetscCall(ISCreateGeneral(comm, len, xinds, PETSC_OWN_POINTER, xis));
  PetscCall(ISCreateGeneral(comm, len, yinds, PETSC_OWN_POINTER, yis));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  ISPairToList - Convert an `IS` pair encoding an integer map to a list of `IS`.

  Collective

  Input Parameters:
+ xis - domain `IS`
- yis - range `IS`, the maximum value must be less than `PETSC_MPI_INT_MAX`

  Output Parameters:
+ listlen - length of `islist`
- islist  - list of `IS`s breaking up indices by color

  Level: developer

  Notes:
  Each `IS` in `islist` contains the preimage for each index on `yis`.
  The `IS` in `islist` are constructed on the subcommunicators of the input `IS`
  pair. Each subcommunicator corresponds to the preimage of some index j -- this subcomm
  contains exactly the MPI processes that assign some indices i to j.  This is essentially the inverse
  of `ISListToPair()`.

  `xis` and `yis` must be of the same length and have congruent communicators.

  The resulting `IS` have subcommunicators in a "deadlock-free" order (see `ISListToPair()`).

.seealso: `IS`, `ISListToPair()`
 @*/
PetscErrorCode ISPairToList(IS xis, IS yis, PetscInt *listlen, IS **islist)
{
  IS              indis = xis, coloris = yis;
  PetscInt       *inds, *colors, llen, ilen, lstart, lend, lcount, l;
  PetscMPIInt     rank, size, llow, lhigh, low, high, color, subsize;
  const PetscInt *ccolors, *cinds;
  MPI_Comm        comm, subcomm;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(xis, IS_CLASSID, 1);
  PetscValidHeaderSpecific(yis, IS_CLASSID, 2);
  PetscCheckSameComm(xis, 1, yis, 2);
  PetscAssertPointer(listlen, 3);
  PetscAssertPointer(islist, 4);
  PetscCall(PetscObjectGetComm((PetscObject)xis, &comm));
  PetscCallMPI(MPI_Comm_rank(comm, &rank));
  PetscCallMPI(MPI_Comm_size(comm, &size));
  /* Extract, copy and sort the local indices and colors on the color. */
  PetscCall(ISGetLocalSize(coloris, &llen));
  PetscCall(ISGetLocalSize(indis, &ilen));
  PetscCheck(llen == ilen, comm, PETSC_ERR_ARG_SIZ, "Incompatible IS sizes: %" PetscInt_FMT " and %" PetscInt_FMT, ilen, llen);
  PetscCall(ISGetIndices(coloris, &ccolors));
  PetscCall(ISGetIndices(indis, &cinds));
  PetscCall(PetscMalloc2(ilen, &inds, llen, &colors));
  PetscCall(PetscArraycpy(inds, cinds, ilen));
  PetscCall(PetscArraycpy(colors, ccolors, llen));
  PetscCall(PetscSortIntWithArray(llen, colors, inds));
  /* Determine the global extent of colors. */
  llow   = 0;
  lhigh  = -1;
  lstart = 0;
  lcount = 0;
  while (lstart < llen) {
    lend = lstart + 1;
    while (lend < llen && colors[lend] == colors[lstart]) ++lend;
    PetscCall(PetscMPIIntCast(PetscMin(llow, colors[lstart]), &llow));
    PetscCall(PetscMPIIntCast(PetscMax(lhigh, colors[lstart]), &lhigh));
    ++lcount;
  }
  PetscCallMPI(MPIU_Allreduce(&llow, &low, 1, MPI_INT, MPI_MIN, comm));
  PetscCallMPI(MPIU_Allreduce(&lhigh, &high, 1, MPI_INT, MPI_MAX, comm));
  *listlen = 0;
  if (low <= high) {
    if (lcount > 0) {
      *listlen = lcount;
      if (!*islist) PetscCall(PetscMalloc1(lcount, islist));
    }
    /*
     Traverse all possible global colors, and participate in the subcommunicators
     for the locally-supported colors.
     */
    lcount = 0;
    lstart = 0;
    lend   = 0;
    for (l = low; l <= high; ++l) {
      /*
       Find the range of indices with the same color, which is not smaller than l.
       Observe that, since colors is sorted, and is a subsequence of [low,high],
       as soon as we find a new color, it is >= l.
       */
      if (lstart < llen) {
        /* The start of the next locally-owned color is identified.  Now look for the end. */
        if (lstart == lend) {
          lend = lstart + 1;
          while (lend < llen && colors[lend] == colors[lstart]) ++lend;
        }
        /* Now check whether the identified color segment matches l. */
        PetscCheck(colors[lstart] >= l, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Locally owned color %" PetscInt_FMT " at location %" PetscInt_FMT " is < than the next global color %" PetscInt_FMT, colors[lstart], lcount, l);
      }
      color = (PetscMPIInt)(colors[lstart] == l);
      /* Check whether a proper subcommunicator exists. */
      PetscCallMPI(MPIU_Allreduce(&color, &subsize, 1, MPI_INT, MPI_SUM, comm));

      if (subsize == 1) subcomm = PETSC_COMM_SELF;
      else if (subsize == size) subcomm = comm;
      else {
        /* a proper communicator is necessary, so we create it. */
        PetscCallMPI(MPI_Comm_split(comm, color, rank, &subcomm));
      }
      if (colors[lstart] == l) {
        /* If we have l among the local colors, we create an IS to hold the corresponding indices. */
        PetscCall(ISCreateGeneral(subcomm, lend - lstart, inds + lstart, PETSC_COPY_VALUES, *islist + lcount));
        /* Position lstart at the beginning of the next local color. */
        lstart = lend;
        /* Increment the counter of the local colors split off into an IS. */
        ++lcount;
      }
      if (subsize > 0 && subsize < size) {
        /*
         Irrespective of color, destroy the split off subcomm:
         a subcomm used in the IS creation above is duplicated
         into a proper PETSc comm.
         */
        PetscCallMPI(MPI_Comm_free(&subcomm));
      }
    } /* for (l = low; l < high; ++l) */
  } /* if (low <= high) */
  PetscCall(PetscFree2(inds, colors));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  ISEmbed - Embed `IS` `a` into `IS` `b` by finding the locations in `b` that have the same indices as in `a`.
  If `c` is the `IS` of these locations, we have `a = b*c`, regarded as a composition of the
  corresponding `ISLocalToGlobalMapping`.

  Not Collective

  Input Parameters:
+ a    - `IS` to embed
. b    - `IS` to embed into
- drop - flag indicating whether to drop indices of `a` that are not in `b`.

  Output Parameter:
. c - local embedding indices

  Level: developer

  Notes:
  If some of the global indices of `a` are not among the indices of `b`, the embedding is impossible. The local indices of `a`
  corresponding to these global indices are either mapped to -1 (if `!drop`) or are omitted (if `drop`). In the former
  case the size of `c` is the same as that of `a`, in the latter case the size of `c` may be smaller.

  The resulting `IS` is sequential, since the index substitution it encodes is purely local.

.seealso: `IS`, `ISLocalToGlobalMapping`
 @*/
PetscErrorCode ISEmbed(IS a, IS b, PetscBool drop, IS *c)
{
  ISLocalToGlobalMapping     ltog;
  ISGlobalToLocalMappingMode gtoltype = IS_GTOLM_DROP;
  PetscInt                   alen, clen, *cindices, *cindices2;
  const PetscInt            *aindices;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(a, IS_CLASSID, 1);
  PetscValidHeaderSpecific(b, IS_CLASSID, 2);
  PetscAssertPointer(c, 4);
  PetscCall(ISLocalToGlobalMappingCreateIS(b, &ltog));
  PetscCall(ISGetLocalSize(a, &alen));
  PetscCall(ISGetIndices(a, &aindices));
  PetscCall(PetscMalloc1(alen, &cindices));
  if (!drop) gtoltype = IS_GTOLM_MASK;
  PetscCall(ISGlobalToLocalMappingApply(ltog, gtoltype, alen, aindices, &clen, cindices));
  PetscCall(ISRestoreIndices(a, &aindices));
  PetscCall(ISLocalToGlobalMappingDestroy(&ltog));
  if (clen != alen) {
    cindices2 = cindices;
    PetscCall(PetscMalloc1(clen, &cindices));
    PetscCall(PetscArraycpy(cindices, cindices2, clen));
    PetscCall(PetscFree(cindices2));
  }
  PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clen, cindices, PETSC_OWN_POINTER, c));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  ISSortPermutation  -  calculate the permutation of the indices into a nondecreasing order.

  Not Collective

  Input Parameters:
+ f      - `IS` to sort
- always - build the permutation even when `f`'s indices are nondecreasing.

  Output Parameter:
. h - permutation or `NULL`, if `f` is nondecreasing and `always` == `PETSC_FALSE`.

  Level: advanced

  Notes:
  Indices in `f` are unchanged. f[h[i]] is the i-th smallest f index.

  If always == `PETSC_FALSE`, an extra check is performed to see whether
  the `f` indices are nondecreasing. `h` is built on `PETSC_COMM_SELF`, since
  the permutation has a local meaning only.

.seealso: `IS`, `ISLocalToGlobalMapping`, `ISSort()`
 @*/
PetscErrorCode ISSortPermutation(IS f, PetscBool always, IS *h)
{
  const PetscInt *findices;
  PetscInt        fsize, *hindices, i;
  PetscBool       isincreasing;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(f, IS_CLASSID, 1);
  PetscAssertPointer(h, 3);
  PetscCall(ISGetLocalSize(f, &fsize));
  PetscCall(ISGetIndices(f, &findices));
  *h = NULL;
  if (!always) {
    isincreasing = PETSC_TRUE;
    for (i = 1; i < fsize; ++i) {
      if (findices[i] <= findices[i - 1]) {
        isincreasing = PETSC_FALSE;
        break;
      }
    }
    if (isincreasing) {
      PetscCall(ISRestoreIndices(f, &findices));
      PetscFunctionReturn(PETSC_SUCCESS);
    }
  }
  PetscCall(PetscMalloc1(fsize, &hindices));
  for (i = 0; i < fsize; ++i) hindices[i] = i;
  PetscCall(PetscSortIntWithPermutation(fsize, findices, hindices));
  PetscCall(ISRestoreIndices(f, &findices));
  PetscCall(ISCreateGeneral(PETSC_COMM_SELF, fsize, hindices, PETSC_OWN_POINTER, h));
  PetscCall(ISSetInfo(*h, IS_PERMUTATION, IS_LOCAL, PETSC_FALSE, PETSC_TRUE));
  PetscFunctionReturn(PETSC_SUCCESS);
}