impls/telescope/telescope.c

#include <petsc/private/petscimpl.h>
#include <petsc/private/matimpl.h>
#include <petsc/private/pcimpl.h>
#include <petscksp.h> /*I "petscksp.h" I*/
#include <petscdm.h>  /*I "petscdm.h" I*/
#include "../src/ksp/pc/impls/telescope/telescope.h"

static PetscBool  cited      = PETSC_FALSE;
static const char citation[] = "@inproceedings{MaySananRuppKnepleySmith2016,\n"
                               "  title     = {Extreme-Scale Multigrid Components within PETSc},\n"
                               "  author    = {Dave A. May and Patrick Sanan and Karl Rupp and Matthew G. Knepley and Barry F. Smith},\n"
                               "  booktitle = {Proceedings of the Platform for Advanced Scientific Computing Conference},\n"
                               "  series    = {PASC '16},\n"
                               "  isbn      = {978-1-4503-4126-4},\n"
                               "  location  = {Lausanne, Switzerland},\n"
                               "  pages     = {5:1--5:12},\n"
                               "  articleno = {5},\n"
                               "  numpages  = {12},\n"
                               "  url       = {https://doi.acm.org/10.1145/2929908.2929913},\n"
                               "  doi       = {10.1145/2929908.2929913},\n"
                               "  acmid     = {2929913},\n"
                               "  publisher = {ACM},\n"
                               "  address   = {New York, NY, USA},\n"
                               "  keywords  = {GPU, HPC, agglomeration, coarse-level solver, multigrid, parallel computing, preconditioning},\n"
                               "  year      = {2016}\n"
                               "}\n";

/*
 default setup mode

 [1a] scatter to (FORWARD)
 x(comm) -> xtmp(comm)
 [1b] local copy (to) ranks with color = 0
 xred(subcomm) <- xtmp

 [2] solve on sub KSP to obtain yred(subcomm)

 [3a] local copy (from) ranks with color = 0
 yred(subcomm) --> xtmp
 [2b] scatter from (REVERSE)
 xtmp(comm) -> y(comm)
*/

/*
  Collective[comm_f]
  Notes
   * Using comm_f = MPI_COMM_NULL will result in an error
   * Using comm_c = MPI_COMM_NULL is valid. If all instances of comm_c are NULL the subcomm is not valid.
   * If any non NULL comm_c communicator cannot map any of its ranks to comm_f, the subcomm is not valid.
*/
static PetscErrorCode PCTelescopeTestValidSubcomm(MPI_Comm comm_f, MPI_Comm comm_c, PetscBool *isvalid)
{
  PetscInt     valid = 1;
  MPI_Group    group_f, group_c;
  PetscMPIInt  count, k, size_f = 0, size_c = 0, size_c_sum = 0;
  PetscMPIInt *ranks_f, *ranks_c;

  PetscFunctionBegin;
  PetscCheck(comm_f != MPI_COMM_NULL, PETSC_COMM_SELF, PETSC_ERR_SUP, "comm_f cannot be MPI_COMM_NULL");

  PetscCallMPI(MPI_Comm_group(comm_f, &group_f));
  if (comm_c != MPI_COMM_NULL) PetscCallMPI(MPI_Comm_group(comm_c, &group_c));

  PetscCallMPI(MPI_Comm_size(comm_f, &size_f));
  if (comm_c != MPI_COMM_NULL) PetscCallMPI(MPI_Comm_size(comm_c, &size_c));

  /* check not all comm_c's are NULL */
  size_c_sum = size_c;
  PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &size_c_sum, 1, MPI_INT, MPI_SUM, comm_f));
  if (size_c_sum == 0) valid = 0;

  /* check we can map at least 1 rank in comm_c to comm_f */
  PetscCall(PetscMalloc1(size_f, &ranks_f));
  PetscCall(PetscMalloc1(size_c, &ranks_c));
  for (k = 0; k < size_f; k++) ranks_f[k] = MPI_UNDEFINED;
  for (k = 0; k < size_c; k++) ranks_c[k] = k;

  /*
   MPI_Group_translate_ranks() returns a non-zero exit code if any rank cannot be translated.
   I do not want the code to terminate immediately if this occurs, rather I want to throw
   the error later (during PCSetUp_Telescope()) via SETERRQ() with a message indicating
   that comm_c is not a valid sub-communicator.
   Hence I purposefully do not call PetscCall() after MPI_Group_translate_ranks().
  */
  count = 0;
  if (comm_c != MPI_COMM_NULL) {
    (void)MPI_Group_translate_ranks(group_c, size_c, ranks_c, group_f, ranks_f);
    for (k = 0; k < size_f; k++) {
      if (ranks_f[k] == MPI_UNDEFINED) count++;
    }
  }
  if (count == size_f) valid = 0;

  PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &valid, 1, MPIU_INT, MPI_MIN, comm_f));
  if (valid == 1) *isvalid = PETSC_TRUE;
  else *isvalid = PETSC_FALSE;

  PetscCall(PetscFree(ranks_f));
  PetscCall(PetscFree(ranks_c));
  PetscCallMPI(MPI_Group_free(&group_f));
  if (comm_c != MPI_COMM_NULL) PetscCallMPI(MPI_Group_free(&group_c));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static DM private_PCTelescopeGetSubDM(PC_Telescope sred)
{
  DM subdm = NULL;

  if (!PCTelescope_isActiveRank(sred)) {
    subdm = NULL;
  } else {
    switch (sred->sr_type) {
    case TELESCOPE_DEFAULT:
      subdm = NULL;
      break;
    case TELESCOPE_DMDA:
      subdm = ((PC_Telescope_DMDACtx *)sred->dm_ctx)->dmrepart;
      break;
    case TELESCOPE_DMPLEX:
      subdm = NULL;
      break;
    case TELESCOPE_COARSEDM:
      if (sred->ksp) PetscCallAbort(PETSC_COMM_SELF, KSPGetDM(sred->ksp, &subdm));
      break;
    }
  }
  return subdm;
}

static PetscErrorCode PCTelescopeSetUp_default(PC pc, PC_Telescope sred)
{
  PetscInt   m, M, bs, st, ed;
  Vec        x, xred, yred, xtmp;
  Mat        B;
  MPI_Comm   comm, subcomm;
  VecScatter scatter;
  IS         isin;
  VecType    vectype;

  PetscFunctionBegin;
  PetscCall(PetscInfo(pc, "PCTelescope: setup (default)\n"));
  comm    = PetscSubcommParent(sred->psubcomm);
  subcomm = PetscSubcommChild(sred->psubcomm);

  PetscCall(PCGetOperators(pc, NULL, &B));
  PetscCall(MatGetSize(B, &M, NULL));
  PetscCall(MatGetBlockSize(B, &bs));
  PetscCall(MatCreateVecs(B, &x, NULL));
  PetscCall(MatGetVecType(B, &vectype)); /* Use the vectype of the matrix used to construct the preconditioner by default */

  xred = NULL;
  m    = 0;
  if (PCTelescope_isActiveRank(sred)) {
    PetscCall(VecCreate(subcomm, &xred));
    PetscCall(VecSetSizes(xred, PETSC_DECIDE, M));
    PetscCall(VecSetBlockSize(xred, bs));
    PetscCall(VecSetType(xred, vectype));
    PetscCall(VecSetFromOptions(xred));
    PetscCall(VecGetLocalSize(xred, &m));
  }

  yred = NULL;
  if (PCTelescope_isActiveRank(sred)) PetscCall(VecDuplicate(xred, &yred));

  PetscCall(VecCreate(comm, &xtmp));
  PetscCall(VecSetSizes(xtmp, m, PETSC_DECIDE));
  PetscCall(VecSetBlockSize(xtmp, bs));
  PetscCall(VecSetType(xtmp, vectype));

  if (PCTelescope_isActiveRank(sred)) {
    PetscCall(VecGetOwnershipRange(xred, &st, &ed));
    PetscCall(ISCreateStride(comm, ed - st, st, 1, &isin));
  } else {
    PetscCall(VecGetOwnershipRange(x, &st, &ed));
    PetscCall(ISCreateStride(comm, 0, st, 1, &isin));
  }
  PetscCall(ISSetBlockSize(isin, bs));

  PetscCall(VecScatterCreate(x, isin, xtmp, NULL, &scatter));

  sred->isin    = isin;
  sred->scatter = scatter;
  sred->xred    = xred;
  sred->yred    = yred;
  sred->xtmp    = xtmp;
  PetscCall(VecDestroy(&x));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCTelescopeMatCreate_default(PC pc, PC_Telescope sred, MatReuse reuse, Mat *A)
{
  MPI_Comm comm, subcomm;
  Mat      Bred, B;
  PetscInt nr, nc, bs;
  IS       isrow, iscol;
  Mat      Blocal, *_Blocal;

  PetscFunctionBegin;
  PetscCall(PetscInfo(pc, "PCTelescope: updating the redundant preconditioned operator (default)\n"));
  PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
  subcomm = PetscSubcommChild(sred->psubcomm);
  PetscCall(PCGetOperators(pc, NULL, &B));
  PetscCall(MatGetSize(B, &nr, &nc));
  isrow = sred->isin;
  PetscCall(ISCreateStride(PETSC_COMM_SELF, nc, 0, 1, &iscol));
  PetscCall(ISSetIdentity(iscol));
  PetscCall(MatGetBlockSizes(B, NULL, &bs));
  PetscCall(ISSetBlockSize(iscol, bs));
  PetscCall(MatSetOption(B, MAT_SUBMAT_SINGLEIS, PETSC_TRUE));
  PetscCall(MatCreateSubMatrices(B, 1, &isrow, &iscol, MAT_INITIAL_MATRIX, &_Blocal));
  Blocal = *_Blocal;
  PetscCall(PetscFree(_Blocal));
  Bred = NULL;
  if (PCTelescope_isActiveRank(sred)) {
    PetscInt mm;

    if (reuse != MAT_INITIAL_MATRIX) Bred = *A;

    PetscCall(MatGetSize(Blocal, &mm, NULL));
    PetscCall(MatCreateMPIMatConcatenateSeqMat(subcomm, Blocal, mm, reuse, &Bred));
  }
  *A = Bred;
  PetscCall(ISDestroy(&iscol));
  PetscCall(MatDestroy(&Blocal));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCTelescopeSubNullSpaceCreate_Telescope(PC pc, PC_Telescope sred, MatNullSpace nullspace, MatNullSpace *sub_nullspace)
{
  PetscBool  has_const;
  const Vec *vecs;
  Vec       *sub_vecs = NULL;
  PetscInt   i, k, n = 0;
  MPI_Comm   subcomm;

  PetscFunctionBegin;
  subcomm = PetscSubcommChild(sred->psubcomm);
  PetscCall(MatNullSpaceGetVecs(nullspace, &has_const, &n, &vecs));

  if (PCTelescope_isActiveRank(sred)) {
    if (n) PetscCall(VecDuplicateVecs(sred->xred, n, &sub_vecs));
  }

  /* copy entries */
  for (k = 0; k < n; k++) {
    const PetscScalar *x_array;
    PetscScalar       *LA_sub_vec;
    PetscInt           st, ed;

    /* pull in vector x->xtmp */
    PetscCall(VecScatterBegin(sred->scatter, vecs[k], sred->xtmp, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(VecScatterEnd(sred->scatter, vecs[k], sred->xtmp, INSERT_VALUES, SCATTER_FORWARD));
    if (sub_vecs) {
      /* copy vector entries into xred */
      PetscCall(VecGetArrayRead(sred->xtmp, &x_array));
      if (sub_vecs[k]) {
        PetscCall(VecGetOwnershipRange(sub_vecs[k], &st, &ed));
        PetscCall(VecGetArray(sub_vecs[k], &LA_sub_vec));
        for (i = 0; i < ed - st; i++) LA_sub_vec[i] = x_array[i];
        PetscCall(VecRestoreArray(sub_vecs[k], &LA_sub_vec));
      }
      PetscCall(VecRestoreArrayRead(sred->xtmp, &x_array));
    }
  }

  if (PCTelescope_isActiveRank(sred)) {
    /* create new (near) nullspace for redundant object */
    PetscCall(MatNullSpaceCreate(subcomm, has_const, n, sub_vecs, sub_nullspace));
    PetscCall(VecDestroyVecs(n, &sub_vecs));
    PetscCheck(!nullspace->remove, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Propagation of custom remove callbacks not supported when propagating (near) nullspaces with PCTelescope");
    PetscCheck(!nullspace->rmctx, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Propagation of custom remove callback context not supported when propagating (near) nullspaces with PCTelescope");
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCTelescopeMatNullSpaceCreate_default(PC pc, PC_Telescope sred, Mat sub_mat)
{
  Mat B;

  PetscFunctionBegin;
  PetscCall(PCGetOperators(pc, NULL, &B));
  /* Propagate the nullspace if it exists */
  {
    MatNullSpace nullspace, sub_nullspace;
    PetscCall(MatGetNullSpace(B, &nullspace));
    if (nullspace) {
      PetscCall(PetscInfo(pc, "PCTelescope: generating nullspace (default)\n"));
      PetscCall(PCTelescopeSubNullSpaceCreate_Telescope(pc, sred, nullspace, &sub_nullspace));
      if (PCTelescope_isActiveRank(sred)) {
        PetscCall(MatSetNullSpace(sub_mat, sub_nullspace));
        PetscCall(MatNullSpaceDestroy(&sub_nullspace));
      }
    }
  }
  /* Propagate the near nullspace if it exists */
  {
    MatNullSpace nearnullspace, sub_nearnullspace;
    PetscCall(MatGetNearNullSpace(B, &nearnullspace));
    if (nearnullspace) {
      PetscCall(PetscInfo(pc, "PCTelescope: generating near nullspace (default)\n"));
      PetscCall(PCTelescopeSubNullSpaceCreate_Telescope(pc, sred, nearnullspace, &sub_nearnullspace));
      if (PCTelescope_isActiveRank(sred)) {
        PetscCall(MatSetNearNullSpace(sub_mat, sub_nearnullspace));
        PetscCall(MatNullSpaceDestroy(&sub_nearnullspace));
      }
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCView_Telescope(PC pc, PetscViewer viewer)
{
  PC_Telescope sred = (PC_Telescope)pc->data;
  PetscBool    isascii, isstring;
  PetscViewer  subviewer;

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
  if (isascii) {
    {
      MPI_Comm    comm, subcomm;
      PetscMPIInt comm_size, subcomm_size;
      DM          dm = NULL, subdm = NULL;

      PetscCall(PCGetDM(pc, &dm));
      subdm = private_PCTelescopeGetSubDM(sred);

      if (sred->psubcomm) {
        comm    = PetscSubcommParent(sred->psubcomm);
        subcomm = PetscSubcommChild(sred->psubcomm);
        PetscCallMPI(MPI_Comm_size(comm, &comm_size));
        PetscCallMPI(MPI_Comm_size(subcomm, &subcomm_size));

        PetscCall(PetscViewerASCIIPushTab(viewer));
        PetscCall(PetscViewerASCIIPrintf(viewer, "PETSc subcomm: parent comm size reduction factor = %" PetscInt_FMT "\n", sred->redfactor));
        PetscCall(PetscViewerASCIIPrintf(viewer, "PETSc subcomm: parent_size = %d , subcomm_size = %d\n", comm_size, subcomm_size));
        switch (sred->subcommtype) {
        case PETSC_SUBCOMM_INTERLACED:
          PetscCall(PetscViewerASCIIPrintf(viewer, "PETSc subcomm: type = %s\n", PetscSubcommTypes[sred->subcommtype]));
          break;
        case PETSC_SUBCOMM_CONTIGUOUS:
          PetscCall(PetscViewerASCIIPrintf(viewer, "PETSc subcomm type = %s\n", PetscSubcommTypes[sred->subcommtype]));
          break;
        default:
          SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "General subcomm type not supported by PCTelescope");
        }
        PetscCall(PetscViewerASCIIPopTab(viewer));
      } else {
        PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
        subcomm = sred->subcomm;
        if (!PCTelescope_isActiveRank(sred)) subcomm = PETSC_COMM_SELF;

        PetscCall(PetscViewerASCIIPushTab(viewer));
        PetscCall(PetscViewerASCIIPrintf(viewer, "subcomm: using user provided sub-communicator\n"));
        PetscCall(PetscViewerASCIIPopTab(viewer));
      }

      PetscCall(PetscViewerGetSubViewer(viewer, subcomm, &subviewer));
      if (PCTelescope_isActiveRank(sred)) {
        PetscCall(PetscViewerASCIIPushTab(subviewer));

        if (dm && sred->ignore_dm) PetscCall(PetscViewerASCIIPrintf(subviewer, "ignoring DM\n"));
        if (sred->ignore_kspcomputeoperators) PetscCall(PetscViewerASCIIPrintf(subviewer, "ignoring KSPComputeOperators\n"));
        switch (sred->sr_type) {
        case TELESCOPE_DEFAULT:
          PetscCall(PetscViewerASCIIPrintf(subviewer, "setup type: default\n"));
          break;
        case TELESCOPE_DMDA:
          PetscCall(PetscViewerASCIIPrintf(subviewer, "setup type: DMDA auto-repartitioning\n"));
          PetscCall(DMView_DA_Short(subdm, subviewer));
          break;
        case TELESCOPE_DMPLEX:
          PetscCall(PetscViewerASCIIPrintf(subviewer, "setup type: DMPLEX auto-repartitioning\n"));
          break;
        case TELESCOPE_COARSEDM:
          PetscCall(PetscViewerASCIIPrintf(subviewer, "setup type: coarse DM\n"));
          break;
        }

        if (dm) {
          PetscObject obj = (PetscObject)dm;
          PetscCall(PetscViewerASCIIPrintf(subviewer, "Parent DM object:"));
          PetscCall(PetscViewerASCIIUseTabs(subviewer, PETSC_FALSE));
          if (obj->type_name) PetscCall(PetscViewerASCIIPrintf(subviewer, " type = %s;", obj->type_name));
          if (obj->name) PetscCall(PetscViewerASCIIPrintf(subviewer, " name = %s;", obj->name));
          if (obj->prefix) PetscCall(PetscViewerASCIIPrintf(subviewer, " prefix = %s", obj->prefix));
          PetscCall(PetscViewerASCIIPrintf(subviewer, "\n"));
          PetscCall(PetscViewerASCIIUseTabs(subviewer, PETSC_TRUE));
        } else {
          PetscCall(PetscViewerASCIIPrintf(subviewer, "Parent DM object: NULL\n"));
        }
        if (subdm) {
          PetscObject obj = (PetscObject)subdm;
          PetscCall(PetscViewerASCIIPrintf(subviewer, "Sub DM object:"));
          PetscCall(PetscViewerASCIIUseTabs(subviewer, PETSC_FALSE));
          if (obj->type_name) PetscCall(PetscViewerASCIIPrintf(subviewer, " type = %s;", obj->type_name));
          if (obj->name) PetscCall(PetscViewerASCIIPrintf(subviewer, " name = %s;", obj->name));
          if (obj->prefix) PetscCall(PetscViewerASCIIPrintf(subviewer, " prefix = %s", obj->prefix));
          PetscCall(PetscViewerASCIIPrintf(subviewer, "\n"));
          PetscCall(PetscViewerASCIIUseTabs(subviewer, PETSC_TRUE));
        } else {
          PetscCall(PetscViewerASCIIPrintf(subviewer, "Sub DM object: NULL\n"));
        }

        PetscCall(KSPView(sred->ksp, subviewer));
        PetscCall(PetscViewerASCIIPopTab(subviewer));
      }
      PetscCall(PetscViewerRestoreSubViewer(viewer, subcomm, &subviewer));
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCSetUp_Telescope(PC pc)
{
  PC_Telescope    sred = (PC_Telescope)pc->data;
  MPI_Comm        comm, subcomm = 0;
  PCTelescopeType sr_type;

  PetscFunctionBegin;
  PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));

  /* Determine type of setup/update */
  if (!pc->setupcalled) {
    PetscBool has_dm, same;
    DM        dm;

    sr_type = TELESCOPE_DEFAULT;
    has_dm  = PETSC_FALSE;
    PetscCall(PCGetDM(pc, &dm));
    if (dm) has_dm = PETSC_TRUE;
    if (has_dm) {
      /* check for dmda */
      PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMDA, &same));
      if (same) {
        PetscCall(PetscInfo(pc, "PCTelescope: found DMDA\n"));
        sr_type = TELESCOPE_DMDA;
      }
      /* check for dmplex */
      PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &same));
      if (same) {
        PetscCall(PetscInfo(pc, "PCTelescope: found DMPLEX\n"));
        sr_type = TELESCOPE_DMPLEX;
      }

      if (sred->use_coarse_dm) {
        PetscCall(PetscInfo(pc, "PCTelescope: using coarse DM\n"));
        sr_type = TELESCOPE_COARSEDM;
      }

      if (sred->ignore_dm) {
        PetscCall(PetscInfo(pc, "PCTelescope: ignoring DM\n"));
        sr_type = TELESCOPE_DEFAULT;
      }
    }
    sred->sr_type = sr_type;
  } else {
    sr_type = sred->sr_type;
  }

  /* set function pointers for repartition setup, matrix creation/update, matrix (near) nullspace, and reset functionality */
  switch (sr_type) {
  case TELESCOPE_DEFAULT:
    sred->pctelescope_setup_type              = PCTelescopeSetUp_default;
    sred->pctelescope_matcreate_type          = PCTelescopeMatCreate_default;
    sred->pctelescope_matnullspacecreate_type = PCTelescopeMatNullSpaceCreate_default;
    sred->pctelescope_reset_type              = NULL;
    break;
  case TELESCOPE_DMDA:
    pc->ops->apply                            = PCApply_Telescope_dmda;
    pc->ops->applyrichardson                  = PCApplyRichardson_Telescope_dmda;
    sred->pctelescope_setup_type              = PCTelescopeSetUp_dmda;
    sred->pctelescope_matcreate_type          = PCTelescopeMatCreate_dmda;
    sred->pctelescope_matnullspacecreate_type = PCTelescopeMatNullSpaceCreate_dmda;
    sred->pctelescope_reset_type              = PCReset_Telescope_dmda;
    break;
  case TELESCOPE_DMPLEX:
    SETERRQ(comm, PETSC_ERR_SUP, "Support for DMPLEX is currently not available");
  case TELESCOPE_COARSEDM:
    pc->ops->apply                            = PCApply_Telescope_CoarseDM;
    pc->ops->applyrichardson                  = PCApplyRichardson_Telescope_CoarseDM;
    sred->pctelescope_setup_type              = PCTelescopeSetUp_CoarseDM;
    sred->pctelescope_matcreate_type          = NULL;
    sred->pctelescope_matnullspacecreate_type = NULL; /* PCTelescopeMatNullSpaceCreate_CoarseDM; */
    sred->pctelescope_reset_type              = PCReset_Telescope_CoarseDM;
    break;
  default:
    SETERRQ(comm, PETSC_ERR_SUP, "Support only provided for: repartitioning an operator; repartitioning a DMDA; or using a coarse DM");
  }

  /* subcomm definition */
  if (!pc->setupcalled) {
    if ((sr_type == TELESCOPE_DEFAULT) || (sr_type == TELESCOPE_DMDA)) {
      if (!sred->psubcomm) {
        PetscCall(PetscSubcommCreate(comm, &sred->psubcomm));
        PetscCall(PetscSubcommSetNumber(sred->psubcomm, sred->redfactor));
        PetscCall(PetscSubcommSetType(sred->psubcomm, sred->subcommtype));
        sred->subcomm = PetscSubcommChild(sred->psubcomm);
      }
    } else { /* query PC for DM, check communicators */
      DM          dm, dm_coarse_partition          = NULL;
      MPI_Comm    comm_fine, comm_coarse_partition = MPI_COMM_NULL;
      PetscMPIInt csize_fine = 0, csize_coarse_partition = 0, cs[2], csg[2], cnt = 0;
      PetscBool   isvalidsubcomm = PETSC_TRUE;

      PetscCall(PCGetDM(pc, &dm));
      comm_fine = PetscObjectComm((PetscObject)dm);
      PetscCall(DMGetCoarseDM(dm, &dm_coarse_partition));
      if (dm_coarse_partition) cnt = 1;
      PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &cnt, 1, MPI_INT, MPI_SUM, comm_fine));
      PetscCheck(cnt != 0, comm_fine, PETSC_ERR_SUP, "Zero instances of a coarse DM were found");

      PetscCallMPI(MPI_Comm_size(comm_fine, &csize_fine));
      if (dm_coarse_partition) {
        comm_coarse_partition = PetscObjectComm((PetscObject)dm_coarse_partition);
        PetscCallMPI(MPI_Comm_size(comm_coarse_partition, &csize_coarse_partition));
      }

      cs[0] = csize_fine;
      cs[1] = csize_coarse_partition;
      PetscCallMPI(MPIU_Allreduce(cs, csg, 2, MPI_INT, MPI_MAX, comm_fine));
      PetscCheck(csg[0] != csg[1], comm_fine, PETSC_ERR_SUP, "Coarse DM uses the same size communicator as the parent DM attached to the PC");

      PetscCall(PCTelescopeTestValidSubcomm(comm_fine, comm_coarse_partition, &isvalidsubcomm));
      PetscCheck(isvalidsubcomm, comm_fine, PETSC_ERR_SUP, "Coarse DM communicator is not a sub-communicator of parentDM->comm");
      sred->subcomm = comm_coarse_partition;
    }
  }
  subcomm = sred->subcomm;

  /* internal KSP */
  if (!pc->setupcalled) {
    const char *prefix;

    if (PCTelescope_isActiveRank(sred)) {
      PetscCall(KSPCreate(subcomm, &sred->ksp));
      PetscCall(KSPSetNestLevel(sred->ksp, pc->kspnestlevel));
      PetscCall(KSPSetErrorIfNotConverged(sred->ksp, pc->erroriffailure));
      PetscCall(PetscObjectIncrementTabLevel((PetscObject)sred->ksp, (PetscObject)pc, 1));
      PetscCall(KSPSetType(sred->ksp, KSPPREONLY));
      PetscCall(PCGetOptionsPrefix(pc, &prefix));
      PetscCall(KSPSetOptionsPrefix(sred->ksp, prefix));
      PetscCall(KSPAppendOptionsPrefix(sred->ksp, "telescope_"));
    }
  }

  /* setup */
  if (!pc->setupcalled && sred->pctelescope_setup_type) PetscCall(sred->pctelescope_setup_type(pc, sred));
  /* update */
  if (!pc->setupcalled) {
    if (sred->pctelescope_matcreate_type) PetscCall(sred->pctelescope_matcreate_type(pc, sred, MAT_INITIAL_MATRIX, &sred->Bred));
    if (sred->pctelescope_matnullspacecreate_type) PetscCall(sred->pctelescope_matnullspacecreate_type(pc, sred, sred->Bred));
  } else {
    if (sred->pctelescope_matcreate_type) PetscCall(sred->pctelescope_matcreate_type(pc, sred, MAT_REUSE_MATRIX, &sred->Bred));
  }

  /* common - no construction */
  if (PCTelescope_isActiveRank(sred)) {
    PetscCall(KSPSetOperators(sred->ksp, sred->Bred, sred->Bred));
    if (pc->setfromoptionscalled && !pc->setupcalled) PetscCall(KSPSetFromOptions(sred->ksp));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCApply_Telescope(PC pc, Vec x, Vec y)
{
  PC_Telescope       sred = (PC_Telescope)pc->data;
  Vec                xtmp, xred, yred;
  PetscInt           i, st, ed;
  VecScatter         scatter;
  PetscScalar       *array;
  const PetscScalar *x_array;

  PetscFunctionBegin;
  PetscCall(PetscCitationsRegister(citation, &cited));

  xtmp    = sred->xtmp;
  scatter = sred->scatter;
  xred    = sred->xred;
  yred    = sred->yred;

  /* pull in vector x->xtmp */
  PetscCall(VecScatterBegin(scatter, x, xtmp, INSERT_VALUES, SCATTER_FORWARD));
  PetscCall(VecScatterEnd(scatter, x, xtmp, INSERT_VALUES, SCATTER_FORWARD));

  /* copy vector entries into xred */
  PetscCall(VecGetArrayRead(xtmp, &x_array));
  if (xred) {
    PetscScalar *LA_xred;
    PetscCall(VecGetOwnershipRange(xred, &st, &ed));
    PetscCall(VecGetArray(xred, &LA_xred));
    for (i = 0; i < ed - st; i++) LA_xred[i] = x_array[i];
    PetscCall(VecRestoreArray(xred, &LA_xred));
  }
  PetscCall(VecRestoreArrayRead(xtmp, &x_array));
  /* solve */
  if (PCTelescope_isActiveRank(sred)) {
    PetscCall(KSPSolve(sred->ksp, xred, yred));
    PetscCall(KSPCheckSolve(sred->ksp, pc, yred));
  }
  /* return vector */
  PetscCall(VecGetArray(xtmp, &array));
  if (yred) {
    const PetscScalar *LA_yred;
    PetscCall(VecGetOwnershipRange(yred, &st, &ed));
    PetscCall(VecGetArrayRead(yred, &LA_yred));
    for (i = 0; i < ed - st; i++) array[i] = LA_yred[i];
    PetscCall(VecRestoreArrayRead(yred, &LA_yred));
  }
  PetscCall(VecRestoreArray(xtmp, &array));
  PetscCall(VecScatterBegin(scatter, xtmp, y, INSERT_VALUES, SCATTER_REVERSE));
  PetscCall(VecScatterEnd(scatter, xtmp, y, INSERT_VALUES, SCATTER_REVERSE));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCApplyRichardson_Telescope(PC pc, Vec x, Vec y, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt its, PetscBool zeroguess, PetscInt *outits, PCRichardsonConvergedReason *reason)
{
  PC_Telescope       sred = (PC_Telescope)pc->data;
  Vec                xtmp, yred;
  PetscInt           i, st, ed;
  VecScatter         scatter;
  const PetscScalar *x_array;
  PetscBool          default_init_guess_value;

  PetscFunctionBegin;
  xtmp    = sred->xtmp;
  scatter = sred->scatter;
  yred    = sred->yred;

  PetscCheck(its <= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PCApplyRichardson_Telescope only supports max_it = 1");
  *reason = (PCRichardsonConvergedReason)0;

  if (!zeroguess) {
    PetscCall(PetscInfo(pc, "PCTelescope: Scattering y for non-zero initial guess\n"));
    /* pull in vector y->xtmp */
    PetscCall(VecScatterBegin(scatter, y, xtmp, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(VecScatterEnd(scatter, y, xtmp, INSERT_VALUES, SCATTER_FORWARD));

    /* copy vector entries into xred */
    PetscCall(VecGetArrayRead(xtmp, &x_array));
    if (yred) {
      PetscScalar *LA_yred;
      PetscCall(VecGetOwnershipRange(yred, &st, &ed));
      PetscCall(VecGetArray(yred, &LA_yred));
      for (i = 0; i < ed - st; i++) LA_yred[i] = x_array[i];
      PetscCall(VecRestoreArray(yred, &LA_yred));
    }
    PetscCall(VecRestoreArrayRead(xtmp, &x_array));
  }

  if (PCTelescope_isActiveRank(sred)) {
    PetscCall(KSPGetInitialGuessNonzero(sred->ksp, &default_init_guess_value));
    if (!zeroguess) PetscCall(KSPSetInitialGuessNonzero(sred->ksp, PETSC_TRUE));
  }

  PetscCall(PCApply_Telescope(pc, x, y));

  if (PCTelescope_isActiveRank(sred)) PetscCall(KSPSetInitialGuessNonzero(sred->ksp, default_init_guess_value));

  if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS;
  *outits = 1;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCReset_Telescope(PC pc)
{
  PC_Telescope sred = (PC_Telescope)pc->data;

  PetscFunctionBegin;
  PetscCall(ISDestroy(&sred->isin));
  PetscCall(VecScatterDestroy(&sred->scatter));
  PetscCall(VecDestroy(&sred->xred));
  PetscCall(VecDestroy(&sred->yred));
  PetscCall(VecDestroy(&sred->xtmp));
  PetscCall(MatDestroy(&sred->Bred));
  PetscCall(KSPReset(sred->ksp));
  if (sred->pctelescope_reset_type) PetscCall(sred->pctelescope_reset_type(pc));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCDestroy_Telescope(PC pc)
{
  PC_Telescope sred = (PC_Telescope)pc->data;

  PetscFunctionBegin;
  PetscCall(PCReset_Telescope(pc));
  PetscCall(KSPDestroy(&sred->ksp));
  PetscCall(PetscSubcommDestroy(&sred->psubcomm));
  PetscCall(PetscFree(sred->dm_ctx));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeGetKSP_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeGetSubcommType_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeSetSubcommType_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeGetReductionFactor_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeSetReductionFactor_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeGetIgnoreDM_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeSetIgnoreDM_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeGetIgnoreKSPComputeOperators_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeSetIgnoreKSPComputeOperators_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeGetDM_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeGetUseCoarseDM_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeSetUseCoarseDM_C", NULL));
  PetscCall(PetscFree(pc->data));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCSetFromOptions_Telescope(PC pc, PetscOptionItems PetscOptionsObject)
{
  PC_Telescope     sred = (PC_Telescope)pc->data;
  MPI_Comm         comm;
  PetscMPIInt      size;
  PetscBool        flg;
  PetscSubcommType subcommtype;

  PetscFunctionBegin;
  PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
  PetscCallMPI(MPI_Comm_size(comm, &size));
  PetscOptionsHeadBegin(PetscOptionsObject, "Telescope options");
  PetscCall(PetscOptionsEnum("-pc_telescope_subcomm_type", "Subcomm type (interlaced or contiguous)", "PCTelescopeSetSubcommType", PetscSubcommTypes, (PetscEnum)sred->subcommtype, (PetscEnum *)&subcommtype, &flg));
  if (flg) PetscCall(PCTelescopeSetSubcommType(pc, subcommtype));
  PetscCall(PetscOptionsInt("-pc_telescope_reduction_factor", "Factor to reduce comm size by", "PCTelescopeSetReductionFactor", sred->redfactor, &sred->redfactor, NULL));
  PetscCheck(sred->redfactor <= size, comm, PETSC_ERR_ARG_WRONG, "-pc_telescope_reduction_factor <= comm size");
  PetscCall(PetscOptionsBool("-pc_telescope_ignore_dm", "Ignore any DM attached to the PC", "PCTelescopeSetIgnoreDM", sred->ignore_dm, &sred->ignore_dm, NULL));
  PetscCall(PetscOptionsBool("-pc_telescope_ignore_kspcomputeoperators", "Ignore method used to compute A", "PCTelescopeSetIgnoreKSPComputeOperators", sred->ignore_kspcomputeoperators, &sred->ignore_kspcomputeoperators, NULL));
  PetscCall(PetscOptionsBool("-pc_telescope_use_coarse_dm", "Define sub-communicator from the coarse DM", "PCTelescopeSetUseCoarseDM", sred->use_coarse_dm, &sred->use_coarse_dm, NULL));
  PetscOptionsHeadEnd();
  PetscFunctionReturn(PETSC_SUCCESS);
}

/* PC simplementation specific API's */

static PetscErrorCode PCTelescopeGetKSP_Telescope(PC pc, KSP *ksp)
{
  PC_Telescope red = (PC_Telescope)pc->data;

  PetscFunctionBegin;
  if (ksp) *ksp = red->ksp;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCTelescopeGetSubcommType_Telescope(PC pc, PetscSubcommType *subcommtype)
{
  PC_Telescope red = (PC_Telescope)pc->data;

  PetscFunctionBegin;
  if (subcommtype) *subcommtype = red->subcommtype;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCTelescopeSetSubcommType_Telescope(PC pc, PetscSubcommType subcommtype)
{
  PC_Telescope red = (PC_Telescope)pc->data;

  PetscFunctionBegin;
  PetscCheck(!pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "You cannot change the subcommunicator type for PCTelescope after it has been set up.");
  red->subcommtype = subcommtype;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCTelescopeGetReductionFactor_Telescope(PC pc, PetscInt *fact)
{
  PC_Telescope red = (PC_Telescope)pc->data;

  PetscFunctionBegin;
  if (fact) *fact = red->redfactor;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCTelescopeSetReductionFactor_Telescope(PC pc, PetscInt fact)
{
  PC_Telescope red = (PC_Telescope)pc->data;
  PetscMPIInt  size;

  PetscFunctionBegin;
  PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
  PetscCheck(fact > 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Reduction factor of telescoping PC %" PetscInt_FMT " must be positive", fact);
  PetscCheck(fact <= size, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Reduction factor of telescoping PC %" PetscInt_FMT " must be <= comm.size", fact);
  red->redfactor = fact;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCTelescopeGetIgnoreDM_Telescope(PC pc, PetscBool *v)
{
  PC_Telescope red = (PC_Telescope)pc->data;

  PetscFunctionBegin;
  if (v) *v = red->ignore_dm;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCTelescopeSetIgnoreDM_Telescope(PC pc, PetscBool v)
{
  PC_Telescope red = (PC_Telescope)pc->data;

  PetscFunctionBegin;
  red->ignore_dm = v;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCTelescopeGetUseCoarseDM_Telescope(PC pc, PetscBool *v)
{
  PC_Telescope red = (PC_Telescope)pc->data;

  PetscFunctionBegin;
  if (v) *v = red->use_coarse_dm;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCTelescopeSetUseCoarseDM_Telescope(PC pc, PetscBool v)
{
  PC_Telescope red = (PC_Telescope)pc->data;

  PetscFunctionBegin;
  red->use_coarse_dm = v;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCTelescopeGetIgnoreKSPComputeOperators_Telescope(PC pc, PetscBool *v)
{
  PC_Telescope red = (PC_Telescope)pc->data;

  PetscFunctionBegin;
  if (v) *v = red->ignore_kspcomputeoperators;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCTelescopeSetIgnoreKSPComputeOperators_Telescope(PC pc, PetscBool v)
{
  PC_Telescope red = (PC_Telescope)pc->data;

  PetscFunctionBegin;
  red->ignore_kspcomputeoperators = v;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCTelescopeGetDM_Telescope(PC pc, DM *dm)
{
  PC_Telescope red = (PC_Telescope)pc->data;

  PetscFunctionBegin;
  *dm = private_PCTelescopeGetSubDM(red);
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCTelescopeGetKSP - Gets the `KSP` created by the telescoping `PC`.

  Not Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameter:
. subksp - the `KSP` defined on the smaller set of processes

  Level: advanced

.seealso: [](ch_ksp), `PC`, `KSP`, `PCTELESCOPE`
@*/
PetscErrorCode PCTelescopeGetKSP(PC pc, KSP *subksp)
{
  PetscFunctionBegin;
  PetscUseMethod(pc, "PCTelescopeGetKSP_C", (PC, KSP *), (pc, subksp));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCTelescopeGetReductionFactor - Gets the factor by which the original number of MPI processes has been reduced by that was set by
  `PCTelescopeSetReductionFactor()`

  Not Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameter:
. fact - the reduction factor

  Level: advanced

.seealso: [](ch_ksp), `PC`, `PCTELESCOPE`, `PCTelescopeSetReductionFactor()`
@*/
PetscErrorCode PCTelescopeGetReductionFactor(PC pc, PetscInt *fact)
{
  PetscFunctionBegin;
  PetscUseMethod(pc, "PCTelescopeGetReductionFactor_C", (PC, PetscInt *), (pc, fact));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCTelescopeSetReductionFactor - Sets the factor by which the original number of MPI processes will been reduced by when
  constructing the subcommunicator to be used with the `PCTELESCOPE`.

  Not Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameter:
. fact - the reduction factor

  Level: advanced

.seealso: [](ch_ksp), `PCTELESCOPE`, `PCTelescopeGetReductionFactor()`
@*/
PetscErrorCode PCTelescopeSetReductionFactor(PC pc, PetscInt fact)
{
  PetscFunctionBegin;
  PetscTryMethod(pc, "PCTelescopeSetReductionFactor_C", (PC, PetscInt), (pc, fact));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCTelescopeGetIgnoreDM - Get the flag indicating if any `DM` attached to the `PC` will be used in constructing the `PC` on the
  reduced number of MPI processes

  Not Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameter:
. v - the flag

  Level: advanced

.seealso: [](ch_ksp), `DM`, `PCTELESCOPE`, `PCTelescopeSetIgnoreDM()`
@*/
PetscErrorCode PCTelescopeGetIgnoreDM(PC pc, PetscBool *v)
{
  PetscFunctionBegin;
  PetscUseMethod(pc, "PCTelescopeGetIgnoreDM_C", (PC, PetscBool *), (pc, v));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCTelescopeSetIgnoreDM - Set a flag to ignore any `DM` attached to the `PC` when constructing the `PC` on the
  reduced number of MPI processes

  Not Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameter:
. v - Use `PETSC_TRUE` to ignore any `DM`

  Level: advanced

.seealso: [](ch_ksp), `DM`, `PCTELESCOPE`, `PCTelescopeGetIgnoreDM()`
@*/
PetscErrorCode PCTelescopeSetIgnoreDM(PC pc, PetscBool v)
{
  PetscFunctionBegin;
  PetscTryMethod(pc, "PCTelescopeSetIgnoreDM_C", (PC, PetscBool), (pc, v));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCTelescopeGetUseCoarseDM - Get the flag indicating if the coarse `DM` attached to `DM` associated with the `PC` will be used in constructing
  the `PC` on the reduced number of MPI processes

  Not Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameter:
. v - the flag

  Level: advanced

.seealso: [](ch_ksp), `DM`, `PCTELESCOPE`, `PCTelescopeSetIgnoreDM()`, `PCTelescopeSetUseCoarseDM()`
@*/
PetscErrorCode PCTelescopeGetUseCoarseDM(PC pc, PetscBool *v)
{
  PetscFunctionBegin;
  PetscUseMethod(pc, "PCTelescopeGetUseCoarseDM_C", (PC, PetscBool *), (pc, v));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCTelescopeSetUseCoarseDM - Set a flag to query the `DM` attached to the `PC` if it also has a coarse `DM` and utilize that `DM`
  in constructing the `PC` on the reduced number of MPI processes

  Not Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameter:
. v - Use `PETSC_FALSE` to ignore any coarse `DM`

  Level: advanced

  Notes:
  When you have specified to use a coarse `DM`, the communicator used to create the sub-`KSP` within `PCTELESCOPE`
  will be that of the coarse `DM`. Hence the flags `-pc_telescope_reduction_factor` and
  `-pc_telescope_subcomm_type` will not be used.

  It is required that the communicator associated with the parent (fine) and the coarse `DM` are of different sizes.
  An error will occur of the size if the communicator associated with the coarse `DM` is the same as that of the parent `DM`.
  Furthermore, it is required that the communicator on the coarse `DM` is a sub-communicator of the parent.
  This will be checked at the time the preconditioner is setup and an error will occur if
  the coarse `DM` does not define a sub-communicator of that used by the parent `DM`.

  The particular Telescope setup invoked when using a coarse `DM` is agnostic with respect to the type of
  the `DM` used (e.g. it supports `DMSHELL`, `DMPLEX`, etc).

  Support is currently only provided for the case when you are using `KSPSetComputeOperators()`

  The user is required to compose a function with the parent `DM` to facilitate the transfer of fields (`Vec`)
  between the different decompositions defined by the fine and coarse `DM`s.
  In the user code, this is achieved via
.vb
   {
     DM dm_fine;
     PetscObjectCompose((PetscObject)dm_fine,"PCTelescopeFieldScatter",your_field_scatter_method);
   }
.ve
  The signature of the user provided field scatter method is
.vb
   PetscErrorCode your_field_scatter_method(DM dm_fine,Vec x_fine,ScatterMode mode,DM dm_coarse,Vec x_coarse);
.ve
  The user must provide support for both mode `SCATTER_FORWARD` and mode `SCATTER_REVERSE`.
  `SCATTER_FORWARD` implies the direction of transfer is from the parent (fine) `DM` to the coarse `DM`.

  Optionally, the user may also compose a function with the parent `DM` to facilitate the transfer
  of state variables between the fine and coarse `DM`s.
  In the context of a finite element discretization, an example state variable might be
  values associated with quadrature points within each element.
  A user provided state scatter method is composed via
.vb
   {
     DM dm_fine;
     PetscObjectCompose((PetscObject)dm_fine,"PCTelescopeStateScatter",your_state_scatter_method);
   }
.ve
  The signature of the user provided state scatter method is
.vb
   PetscErrorCode your_state_scatter_method(DM dm_fine,ScatterMode mode,DM dm_coarse);
.ve
  `SCATTER_FORWARD` implies the direction of transfer is from the fine `DM` to the coarse `DM`.
  The user is only required to support mode = `SCATTER_FORWARD`.
  No assumption is made about the data type of the state variables.
  These must be managed by the user and must be accessible from the `DM`.

  Care must be taken in defining the user context passed to `KSPSetComputeOperators()` which is to be
  associated with the sub-`KSP` residing within `PCTELESCOPE`.
  In general, `PCTELESCOPE` assumes that the context on the fine and coarse `DM` used with
  `KSPSetComputeOperators()` should be "similar" in type or origin.
  Specifically the following rules are used to infer what context on the sub-`KSP` should be.

  First the contexts from the `KSP` and the fine and coarse `DM`s are retrieved.
  Note that the special case of a `DMSHELL` context is queried.

.vb
   DMKSPGetComputeOperators(dm_fine,&dmfine_kspfunc,&dmfine_kspctx);
   DMGetApplicationContext(dm_fine,&dmfine_appctx);
   DMShellGetContext(dm_fine,&dmfine_shellctx);

   DMGetApplicationContext(dm_coarse,&dmcoarse_appctx);
   DMShellGetContext(dm_coarse,&dmcoarse_shellctx);
.ve

  The following rules are then enforced\:

  1. If `dmfine_kspctx` = `NULL`, then we provide a `NULL` pointer as the context for the sub-`KSP`\:
  `KSPSetComputeOperators`(`sub_ksp`,`dmfine_kspfunc`,`NULL`);

  2. If `dmfine_kspctx` != `NULL` and `dmfine_kspctx` == `dmfine_appctx`,

  check that `dmcoarse_appctx` is also non-`NULL`. If this is true, then\:
  `KSPSetComputeOperators`(`sub_ksp`,`dmfine_kspfunc`,`dmcoarse_appctx`);

  3. If `dmfine_kspctx` != `NULL` and `dmfine_kspctx` == `dmfine_shellctx`,

  check that `dmcoarse_shellctx` is also non-`NULL`. If this is true, then\:
  `KSPSetComputeOperators`(`sub_ksp`,`dmfine_kspfunc`,`dmcoarse_shellctx`);

  If neither of the above three tests passed, then `PCTELESCOPE` cannot safely determine what
  context should be provided to `KSPSetComputeOperators()` for use with the sub-`KSP`.
  In this case, an additional mechanism is provided via a composed function which will return
  the actual context to be used. To use this feature you must compose the "getter" function
  with the coarse `DM`, e.g.
.vb
   {
     DM dm_coarse;
     PetscObjectCompose((PetscObject)dm_coarse,"PCTelescopeGetCoarseDMKSPContext",your_coarse_context_getter);
   }
.ve
  The signature of the user provided method is
.vb
   PetscErrorCode your_coarse_context_getter(DM dm_coarse,void **your_kspcontext);
.ve

.seealso: [](ch_ksp), `DM`, `PCTELESCOPE`, `PCTelescopeSetIgnoreDM()`
@*/
PetscErrorCode PCTelescopeSetUseCoarseDM(PC pc, PetscBool v)
{
  PetscFunctionBegin;
  PetscTryMethod(pc, "PCTelescopeSetUseCoarseDM_C", (PC, PetscBool), (pc, v));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCTelescopeGetIgnoreKSPComputeOperators - Get the flag indicating if `KSPComputeOperators()` will be used to construct
  the matrix on the reduced number of MPI processes

  Not Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameter:
. v - the flag

  Level: advanced

.seealso: [](ch_ksp), `PCTELESCOPE`, `PCTelescopeSetIgnoreDM()`, `PCTelescopeSetUseCoarseDM()`, `PCTelescopeSetIgnoreKSPComputeOperators()`
@*/
PetscErrorCode PCTelescopeGetIgnoreKSPComputeOperators(PC pc, PetscBool *v)
{
  PetscFunctionBegin;
  PetscUseMethod(pc, "PCTelescopeGetIgnoreKSPComputeOperators_C", (PC, PetscBool *), (pc, v));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCTelescopeSetIgnoreKSPComputeOperators - Set a flag to have `PCTELESCOPE` ignore the function provided to `KSPComputeOperators()` in
  constructint the matrix on the reduced number of MPI processes

  Not Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameter:
. v - Use `PETSC_TRUE` to ignore the function (if defined) set via `KSPSetComputeOperators()` on `pc`

  Level: advanced

.seealso: [](ch_ksp), `PCTELESCOPE`, `PCTelescopeSetIgnoreDM()`, `PCTelescopeSetUseCoarseDM()`, `PCTelescopeGetIgnoreKSPComputeOperators()`
@*/
PetscErrorCode PCTelescopeSetIgnoreKSPComputeOperators(PC pc, PetscBool v)
{
  PetscFunctionBegin;
  PetscTryMethod(pc, "PCTelescopeSetIgnoreKSPComputeOperators_C", (PC, PetscBool), (pc, v));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCTelescopeGetDM - Get the re-partitioned `DM` attached to the sub-`KSP`.

  Not Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameter:
. subdm - The re-partitioned `DM`

  Level: advanced

.seealso: [](ch_ksp), `DM`, `PCTELESCOPE`, `PCTelescopeSetIgnoreDM()`, `PCTelescopeSetUseCoarseDM()`, `PCTelescopeGetIgnoreKSPComputeOperators()`
@*/
PetscErrorCode PCTelescopeGetDM(PC pc, DM *subdm)
{
  PetscFunctionBegin;
  PetscUseMethod(pc, "PCTelescopeGetDM_C", (PC, DM *), (pc, subdm));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCTelescopeSetSubcommType - set subcommunicator type `PetscSubcommType` (interlaced or contiguous) to be used when
  the subcommunicator is generated from the given `PC`

  Logically Collective

  Input Parameters:
+ pc          - the preconditioner context
- subcommtype - the subcommunicator type (see `PetscSubcommType`)

  Level: advanced

.seealso: [](ch_ksp), `PetscSubcommType`, `PetscSubcomm`, `PCTELESCOPE`, `PCTelescopeGetSubcommType()`
@*/
PetscErrorCode PCTelescopeSetSubcommType(PC pc, PetscSubcommType subcommtype)
{
  PetscFunctionBegin;
  PetscTryMethod(pc, "PCTelescopeSetSubcommType_C", (PC, PetscSubcommType), (pc, subcommtype));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCTelescopeGetSubcommType - Get the subcommunicator type `PetscSubcommType` (interlaced or contiguous) set with `PCTelescopeSetSubcommType()`

  Not Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameter:
. subcommtype - the subcommunicator type (see `PetscSubcommType`)

  Level: advanced

.seealso: [](ch_ksp), `PetscSubcomm`, `PetscSubcommType`, `PCTELESCOPE`, `PCTelescopeSetSubcommType()`
@*/
PetscErrorCode PCTelescopeGetSubcommType(PC pc, PetscSubcommType *subcommtype)
{
  PetscFunctionBegin;
  PetscUseMethod(pc, "PCTelescopeGetSubcommType_C", (PC, PetscSubcommType *), (pc, subcommtype));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*MC
   PCTELESCOPE - Runs a `KSP` solver on a sub-communicator {cite}`maysananruppknepleysmith2016` of the communicator used by the original `KSP`.
                 MPI processes not in the sub-communicator are idle during the solve. Usually used to solve the smaller coarser grid problems in multigrid
                 (`PCMG`) that could not be efficiently solved on the entire communication

   Options Database Keys:
+  -pc_telescope_reduction_factor <r>                 - factor to reduce the communicator size by. e.g. with 64 MPI ranks and r=4, the new sub-communicator will have 64/4 = 16 ranks.
.  -pc_telescope_ignore_dm                            - flag to indicate whether an attached `DM` should be ignored in constructing the new `PC`
.  -pc_telescope_subcomm_type <interlaced,contiguous> - defines the selection of MPI processes on the sub-communicator. see `PetscSubcomm` for more information.
.  -pc_telescope_ignore_kspcomputeoperators           - flag to indicate whether `KSPSetComputeOperators()` should be used on the sub-`KSP`.
-  -pc_telescope_use_coarse_dm                        - flag to indicate whether the coarse `DM` should be used to define the sub-communicator.

   Level: advanced

   Notes:
   Assuming that the parent preconditioner `PC` is defined on a communicator c, this implementation
   creates a child sub-communicator (c') containing fewer MPI processes than the original parent preconditioner `PC`.
   The preconditioner is deemed telescopic as it only calls `KSPSolve()` on a single
   sub-communicator, in contrast with `PCREDUNDANT` which calls `KSPSolve()` on N sub-communicators.
   This means there will be MPI processes which will be idle during the application of this preconditioner.
   Additionally, in comparison with `PCREDUNDANT`, `PCTELESCOPE` can utilize an attached `DM` to construct `DM` dependent preconditioner, such as `PCMG`

   The default type `KSPType` of the sub `KSP` (the `KSP` defined on c') is `KSPPREONLY`.

   There are three setup mechanisms for `PCTELESCOPE`. Features support by each type are described below.
   In the following, we will refer to the operators B and B', these are the `Bmat` provided to the `KSP` on the
   communicators c and c' respectively.

   [1] Default setup
   The sub-communicator c' is created via `PetscSubcommCreate()`.
   Any explicitly defined nullspace and near nullspace vectors attached to B with `MatSetNullSpace()` and `MatSetNearNullSpace()` are transferred to B'.
   Currently there is no support for nullspaces provided with `MatNullSpaceSetFunction()`).
   No support is provided for `KSPSetComputeOperators()`.
   Currently there is no support for the flag `-pc_use_amat`.

   [2] `DM` aware setup
   The sub-communicator c' is created via `PetscSubcommCreate()`.
   If a `DM` is attached to the `PC`, it is re-partitioned on the sub-communicator c'.
   Both the `Bmat` operator and the right-hand side vector are permuted into the new DOF ordering defined by the re-partitioned `DM`.
   Currently only support for re-partitioning a `DMDA` is provided.
   Any explicitly defined nullspace or near nullspace vectors attached to the original B with `MatSetNullSpace()`
   and `MatSetNearNullSpace()` are extracted, re-partitioned and set on B'
   (currently there is no support for nullspaces provided with `MatNullSpaceSetFunction()`).
   Support is provided for `KSPSetComputeOperators()`. The user provided function and context is propagated to the sub `KSP`.
   This is fragile since the user must ensure that their user context is valid for use on c'.
   Currently there is no support for the flag `-pc_use_amat`.

   [3] Coarse `DM` setup
   If a `DM` (dmfine) is attached to the `PC`, dmfine is queried for a "coarse" `DM` (call this dmcoarse) via `DMGetCoarseDM()`.
   `PCTELESCOPE` will interpret the coarse `DM` as being defined on a sub-communicator of c.
   The communicator associated with dmcoarse will define the c' to be used within `PCTELESCOPE`.
   `PCTELESCOPE` will check that c' is in fact a sub-communicator of c. If it is not, an error will be reported.
   The intention of this setup type is that `PCTELESCOPE` will use an existing (e.g. user defined) communicator hierarchy, say as would be
   available with using multi-grid on unstructured meshes.
   This setup will not use the command line options `-pc_telescope_reduction_factor` or `-pc_telescope_subcomm_type`.
   Any explicitly defined nullspace or near nullspace vectors attached to the B are extracted, scattered into the correct ordering consistent
   with dmcoarse and set on B'
   (currently there is no support for nullspaces provided with `MatNullSpaceSetFunction()`).
   There is no general method to permute field orderings, hence only `KSPSetComputeOperators()` is supported.
   The user must use `PetscObjectComposeFunction()` with dmfine to define the method to scatter fields from dmfine to dmcoarse.
   Propagation of the user context for `KSPSetComputeOperators()` on the sub `KSP` is attempted by querying the `DM` contexts associated with
   dmfine and dmcoarse. Alternatively, the user may use `PetscObjectComposeFunction()` with dmcoarse to define a method which will return the appropriate user context for `KSPSetComputeOperators()`.
   Currently there is no support for the flag `-pc_use_amat`.
   This setup can be invoked by the option `-pc_telescope_use_coarse_dm` or by calling `PCTelescopeSetUseCoarseDM`(pc,`PETSC_TRUE`);
   Further information about the user-provided methods required by this setup type are described here `PCTelescopeSetUseCoarseDM()`.

   Developer Notes:
   During `PCSetUp()`, the B operator is scattered onto c'.
   Within `PCApply()`, the RHS vector (x) is scattered into a redundant vector, xred (defined on c').
   Then, `KSPSolve()` is executed on the c' communicator.

   The communicator used within the telescoping preconditioner is defined by a `PetscSubcomm` using the INTERLACED
   creation routine by default (this can be changed with `-pc_telescope_subcomm_type`). We run the sub `KSP` on only
   the ranks within the communicator which have a color equal to zero.

   The telescoping preconditioner is aware of nullspaces and near nullspaces which are attached to the B operator.
   In the case where B has a (near) nullspace attached, the (near) nullspace vectors are extracted from B and mapped into
   a new (near) nullspace, defined on the sub-communicator, which is attached to B' (the B operator which was scattered to c')

   The telescoping preconditioner can re-partition an attached `DM` if it is a `DMDA` (2D or 3D -
   support for 1D `DMDA`s is not provided). If a `DMDA` is found, a topologically equivalent `DMDA` is created on c'
   and this new `DM` is attached the sub `KSP`. The design of telescope is such that it should be possible to extend support
   for re-partitioning other to `DM`'s (e.g. `DMPLEX`). The user can supply a flag to ignore attached DMs.
   Alternatively, user-provided re-partitioned `DM`s can be used via `-pc_telescope_use_coarse_dm`.

   With the default setup mode, B' is defined by fusing rows (in order) associated with MPI processes common to c and c'.

   When a `DMDA` is attached to the parent preconditioner, B' is defined by: (i) performing a symmetric permutation of B
   into the ordering defined by the `DMDA` on c', (ii) extracting the local chunks via `MatCreateSubMatrices()`, (iii) fusing the
   locally (sequential) matrices defined on the ranks common to c and c' into B' using `MatCreateMPIMatConcatenateSeqMat()`

   Limitations/improvements include the following.
   `VecPlaceArray()` could be used within `PCApply()` to improve efficiency and reduce memory usage.
   A unified mechanism to query for user contexts as required by `KSPSetComputeOperators()` and `MatNullSpaceSetFunction()`.

   The symmetric permutation used when a `DMDA` is encountered is performed via explicitly assembling a permutation matrix P,
   and performing P^T.A.P. Possibly it might be more efficient to use `MatPermute()`. We opted to use P^T.A.P as it appears
   `VecPermute()` does not support the use case required here. By computing P, one can permute both the operator and RHS in a
   consistent manner.

   Mapping of vectors (default setup mode) is performed in the following way.
   Suppose the parent communicator size was 4, and we set a reduction factor of 2; this would give a comm size on c' of 2.
   Using the interlaced creation routine, the ranks in c with color = 0 will be rank 0 and 2.
   We perform the scatter to the sub-communicator in the following way.
   [1] Given a vector x defined on communicator c

.vb
   rank(c)  local values of x
   ------- ----------------------------------------
        0   [  0.0,  1.0,  2.0,  3.0,  4.0,  5.0 ]
        1   [  6.0,  7.0,  8.0,  9.0, 10.0, 11.0 ]
        2   [ 12.0, 13.0, 14.0, 15.0, 16.0, 17.0 ]
        3   [ 18.0, 19.0, 20.0, 21.0, 22.0, 23.0 ]
.ve

   scatter into xtmp defined also on comm c, so that we have the following values

.vb
   rank(c)  local values of xtmp
   ------- ----------------------------------------------------------------------------
        0   [  0.0,  1.0,  2.0,  3.0,  4.0,  5.0,  6.0,  7.0,  8.0,  9.0, 10.0, 11.0 ]
        1   [ ]
        2   [ 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0 ]
        3   [ ]
.ve

   The entries on rank 1 and 3 (ranks which do not have a color = 0 in c') have no values

   [2] Copy the values from ranks 0, 2 (indices with respect to comm c) into the vector xred which is defined on communicator c'.
   Ranks 0 and 2 are the only ranks in the subcomm which have a color = 0.

.vb
   rank(c')  local values of xred
   -------- ----------------------------------------------------------------------------
         0   [  0.0,  1.0,  2.0,  3.0,  4.0,  5.0,  6.0,  7.0,  8.0,  9.0, 10.0, 11.0 ]
         1   [ 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0 ]
.ve

  Contributed by:
  Dave May

.seealso: [](ch_ksp), `PCTelescopeGetKSP()`, `PCTelescopeGetDM()`, `PCTelescopeGetReductionFactor()`, `PCTelescopeSetReductionFactor()`, `PCTelescopeGetIgnoreDM()`, `PCTelescopeSetIgnoreDM()`, `PCREDUNDANT`
M*/
PETSC_EXTERN PetscErrorCode PCCreate_Telescope(PC pc)
{
  struct _PC_Telescope *sred;

  PetscFunctionBegin;
  PetscCall(PetscNew(&sred));
  sred->psubcomm                   = NULL;
  sred->subcommtype                = PETSC_SUBCOMM_INTERLACED;
  sred->subcomm                    = MPI_COMM_NULL;
  sred->redfactor                  = 1;
  sred->ignore_dm                  = PETSC_FALSE;
  sred->ignore_kspcomputeoperators = PETSC_FALSE;
  sred->use_coarse_dm              = PETSC_FALSE;
  pc->data                         = (void *)sred;

  pc->ops->apply           = PCApply_Telescope;
  pc->ops->applytranspose  = NULL;
  pc->ops->applyrichardson = PCApplyRichardson_Telescope;
  pc->ops->setup           = PCSetUp_Telescope;
  pc->ops->destroy         = PCDestroy_Telescope;
  pc->ops->reset           = PCReset_Telescope;
  pc->ops->setfromoptions  = PCSetFromOptions_Telescope;
  pc->ops->view            = PCView_Telescope;

  sred->pctelescope_setup_type              = PCTelescopeSetUp_default;
  sred->pctelescope_matcreate_type          = PCTelescopeMatCreate_default;
  sred->pctelescope_matnullspacecreate_type = PCTelescopeMatNullSpaceCreate_default;
  sred->pctelescope_reset_type              = NULL;

  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeGetKSP_C", PCTelescopeGetKSP_Telescope));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeGetSubcommType_C", PCTelescopeGetSubcommType_Telescope));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeSetSubcommType_C", PCTelescopeSetSubcommType_Telescope));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeGetReductionFactor_C", PCTelescopeGetReductionFactor_Telescope));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeSetReductionFactor_C", PCTelescopeSetReductionFactor_Telescope));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeGetIgnoreDM_C", PCTelescopeGetIgnoreDM_Telescope));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeSetIgnoreDM_C", PCTelescopeSetIgnoreDM_Telescope));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeGetIgnoreKSPComputeOperators_C", PCTelescopeGetIgnoreKSPComputeOperators_Telescope));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeSetIgnoreKSPComputeOperators_C", PCTelescopeSetIgnoreKSPComputeOperators_Telescope));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeGetDM_C", PCTelescopeGetDM_Telescope));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeGetUseCoarseDM_C", PCTelescopeGetUseCoarseDM_Telescope));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCTelescopeSetUseCoarseDM_C", PCTelescopeSetUseCoarseDM_Telescope));
  PetscFunctionReturn(PETSC_SUCCESS);
}