xref: /petsc/src/ksp/pc/impls/fieldsplit/fieldsplit.c (revision d484b38499953fbf56fc48e9f43f331c10de329f)

#include <petsc/private/pcimpl.h>  /*I "petscpc.h" I*/
#include <petsc/private/kspimpl.h> /*  This is needed to provide the appropriate PETSC_EXTERN for KSP_Solve_FS ....*/
#include <petsc/private/matimpl.h> /* MatScatterDense_Private() for PCMatApply() */
#include <petscdm.h>
#include <petscdevice.h>
#if PetscDefined(HAVE_CUDA)
  #include <petscdevice_cuda.h>
#endif
#if PetscDefined(HAVE_HIP)
  #include <petscdevice_hip.h>
#endif

const char *const PCFieldSplitSchurPreTypes[]  = {"SELF", "SELFP", "A11", "USER", "FULL", "PCFieldSplitSchurPreType", "PC_FIELDSPLIT_SCHUR_PRE_", NULL};
const char *const PCFieldSplitSchurFactTypes[] = {"DIAG", "LOWER", "UPPER", "FULL", "PCFieldSplitSchurFactType", "PC_FIELDSPLIT_SCHUR_FACT_", NULL};

PetscLogEvent KSP_Solve_FS_0, KSP_Solve_FS_1, KSP_Solve_FS_S, KSP_Solve_FS_U, KSP_Solve_FS_L, KSP_Solve_FS_2, KSP_Solve_FS_3, KSP_Solve_FS_4;

typedef struct _PC_FieldSplitLink *PC_FieldSplitLink;
struct _PC_FieldSplitLink {
  KSP               ksp;
  Vec               x, y, z;
  Mat               X, Y;
  char             *splitname;
  PetscInt          nfields;
  PetscInt         *fields, *fields_col;
  VecScatter        sctx;
  IS                is, is_col;
  PC_FieldSplitLink next, previous;
  PetscLogEvent     event;

  /* Used only when setting coordinates with PCSetCoordinates */
  PetscInt   dim;
  PetscInt   ndofs;
  PetscReal *coords;
};

typedef struct {
  PCCompositeType type;
  PetscBool       defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */
  PetscBool       splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */
  PetscInt        bs;           /* Block size for IS and Mat structures */
  PetscInt        nsplits;      /* Number of field divisions defined */
  Vec            *x, *y, w1, w2;
  Mat            *mat;    /* The diagonal block for each split */
  Mat            *pmat;   /* The preconditioning diagonal block for each split */
  Mat            *Afield; /* The rows of the matrix associated with each split */
  PetscBool       issetup;

  /* Only used when Schur complement preconditioning is used */
  Mat                       B;          /* The (0,1) block */
  Mat                       C;          /* The (1,0) block */
  Mat                       schur;      /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */
  Mat                       schurp;     /* Assembled approximation to S built by MatSchurComplement to be used as a matrix for constructing the preconditioner when solving with S */
  Mat                       schur_user; /* User-provided matrix for constructing the preconditioner for the Schur complement */
  PCFieldSplitSchurPreType  schurpre;   /* Determines which matrix is used for the Schur complement */
  PCFieldSplitSchurFactType schurfactorization;
  KSP                       kspschur;   /* The solver for S */
  KSP                       kspupper;   /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */
  PetscScalar               schurscale; /* Scaling factor for the Schur complement solution with DIAG factorization */

  /* Only used when Golub-Kahan bidiagonalization preconditioning is used */
  Mat          H;           /* The modified matrix H = A00 + nu*A01*A01'              */
  PetscReal    gkbtol;      /* Stopping tolerance for lower bound estimate            */
  PetscInt     gkbdelay;    /* The delay window for the stopping criterion            */
  PetscReal    gkbnu;       /* Parameter for augmented Lagrangian H = A + nu*A01*A01' */
  PetscInt     gkbmaxit;    /* Maximum number of iterations for outer loop            */
  PetscBool    gkbmonitor;  /* Monitor for gkb iterations and the lower bound error   */
  PetscViewer  gkbviewer;   /* Viewer context for gkbmonitor                          */
  Vec          u, v, d, Hu; /* Work vectors for the GKB algorithm                     */
  PetscScalar *vecz;        /* Contains intermediate values, eg for lower bound       */

  PC_FieldSplitLink head;
  PetscBool         isrestrict;       /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */
  PetscBool         suboptionsset;    /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */
  PetscBool         dm_splits;        /* Whether to use DMCreateFieldDecomposition() whenever possible */
  PetscBool         diag_use_amat;    /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
  PetscBool         offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
  PetscBool         detect;           /* Whether to form 2-way split by finding zero diagonal entries */
  PetscBool         coordinates_set;  /* Whether PCSetCoordinates has been called */
} PC_FieldSplit;

/*
    Note:
    there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of
   inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the
   PC you could change this.
*/

/* This helper is so that setting a user-provided matrix is orthogonal to choosing to use it.  This way the
* application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */
static Mat FieldSplitSchurPre(PC_FieldSplit *jac)
{
  switch (jac->schurpre) {
  case PC_FIELDSPLIT_SCHUR_PRE_SELF:
    return jac->schur;
  case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
    return jac->schurp;
  case PC_FIELDSPLIT_SCHUR_PRE_A11:
    return jac->pmat[1];
  case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */
  case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */
  default:
    return jac->schur_user ? jac->schur_user : jac->pmat[1];
  }
}

#include <petscdraw.h>
static PetscErrorCode PCView_FieldSplit(PC pc, PetscViewer viewer)
{
  PC_FieldSplit    *jac = (PC_FieldSplit *)pc->data;
  PetscBool         isascii, isdraw;
  PetscInt          i, j;
  PC_FieldSplitLink ilink = jac->head;

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
  if (isascii) {
    if (jac->bs > 0) {
      PetscCall(PetscViewerASCIIPrintf(viewer, "  FieldSplit with %s composition: total splits = %" PetscInt_FMT ", blocksize = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits, jac->bs));
    } else {
      PetscCall(PetscViewerASCIIPrintf(viewer, "  FieldSplit with %s composition: total splits = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits));
    }
    if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, "  using Amat (not Pmat) as operator for blocks\n"));
    if (jac->diag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, "  using Amat (not Pmat) as operator for diagonal blocks\n"));
    if (jac->offdiag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, "  using Amat (not Pmat) as operator for off-diagonal blocks\n"));
    PetscCall(PetscViewerASCIIPrintf(viewer, "  Solver info for each split is in the following KSP objects:\n"));
    for (i = 0; i < jac->nsplits; i++) {
      if (ilink->fields) {
        PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Fields ", i));
        PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
        for (j = 0; j < ilink->nfields; j++) {
          if (j > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ","));
          PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j]));
        }
        PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
        PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
      } else {
        PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Defined by IS\n", i));
      }
      PetscCall(KSPView(ilink->ksp, viewer));
      ilink = ilink->next;
    }
  }

  if (isdraw) {
    PetscDraw draw;
    PetscReal x, y, w, wd;

    PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
    PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
    w  = 2 * PetscMin(1.0 - x, x);
    wd = w / (jac->nsplits + 1);
    x  = x - wd * (jac->nsplits - 1) / 2.0;
    for (i = 0; i < jac->nsplits; i++) {
      PetscCall(PetscDrawPushCurrentPoint(draw, x, y));
      PetscCall(KSPView(ilink->ksp, viewer));
      PetscCall(PetscDrawPopCurrentPoint(draw));
      x += wd;
      ilink = ilink->next;
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCView_FieldSplit_Schur(PC pc, PetscViewer viewer)
{
  PC_FieldSplit             *jac = (PC_FieldSplit *)pc->data;
  PetscBool                  isascii, isdraw;
  PetscInt                   i, j;
  PC_FieldSplitLink          ilink = jac->head;
  MatSchurComplementAinvType atype;

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
  if (isascii) {
    if (jac->bs > 0) {
      PetscCall(PetscViewerASCIIPrintf(viewer, "  FieldSplit with Schur preconditioner, blocksize = %" PetscInt_FMT ", factorization %s\n", jac->bs, PCFieldSplitSchurFactTypes[jac->schurfactorization]));
    } else {
      PetscCall(PetscViewerASCIIPrintf(viewer, "  FieldSplit with Schur preconditioner, factorization %s\n", PCFieldSplitSchurFactTypes[jac->schurfactorization]));
    }
    if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, "  using Amat (not Pmat) as operator for blocks\n"));
    switch (jac->schurpre) {
    case PC_FIELDSPLIT_SCHUR_PRE_SELF:
      PetscCall(PetscViewerASCIIPrintf(viewer, "  Preconditioner for the Schur complement formed from S itself\n"));
      break;
    case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
      if (jac->schur) {
        PetscCall(MatSchurComplementGetAinvType(jac->schur, &atype));
        PetscCall(PetscViewerASCIIPrintf(viewer, "  Preconditioner for the Schur complement formed from Sp, an assembled approximation to S, which uses A00's %sinverse\n", atype == MAT_SCHUR_COMPLEMENT_AINV_DIAG ? "diagonal's " : (atype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG ? "block diagonal's " : (atype == MAT_SCHUR_COMPLEMENT_AINV_FULL ? "full " : "lumped diagonal's "))));
      }
      break;
    case PC_FIELDSPLIT_SCHUR_PRE_A11:
      PetscCall(PetscViewerASCIIPrintf(viewer, "  Preconditioner for the Schur complement formed from A11\n"));
      break;
    case PC_FIELDSPLIT_SCHUR_PRE_FULL:
      PetscCall(PetscViewerASCIIPrintf(viewer, "  Preconditioner for the Schur complement formed from the exact Schur complement\n"));
      break;
    case PC_FIELDSPLIT_SCHUR_PRE_USER:
      if (jac->schur_user) {
        PetscCall(PetscViewerASCIIPrintf(viewer, "  Preconditioner for the Schur complement formed from user provided matrix\n"));
      } else {
        PetscCall(PetscViewerASCIIPrintf(viewer, "  Preconditioner for the Schur complement formed from A11\n"));
      }
      break;
    default:
      SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre);
    }
    PetscCall(PetscViewerASCIIPrintf(viewer, "  Split info:\n"));
    PetscCall(PetscViewerASCIIPushTab(viewer));
    for (i = 0; i < jac->nsplits; i++) {
      if (ilink->fields) {
        PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Fields ", i));
        PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
        for (j = 0; j < ilink->nfields; j++) {
          if (j > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ","));
          PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j]));
        }
        PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
        PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
      } else {
        PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Defined by IS\n", i));
      }
      ilink = ilink->next;
    }
    PetscCall(PetscViewerASCIIPrintf(viewer, "KSP solver for A00 block\n"));
    PetscCall(PetscViewerASCIIPushTab(viewer));
    if (jac->head) PetscCall(KSPView(jac->head->ksp, viewer));
    else PetscCall(PetscViewerASCIIPrintf(viewer, "  not yet available\n"));
    PetscCall(PetscViewerASCIIPopTab(viewer));
    if (jac->head && jac->kspupper != jac->head->ksp) {
      PetscCall(PetscViewerASCIIPrintf(viewer, "KSP solver for upper A00 in upper triangular factor\n"));
      PetscCall(PetscViewerASCIIPushTab(viewer));
      if (jac->kspupper) PetscCall(KSPView(jac->kspupper, viewer));
      else PetscCall(PetscViewerASCIIPrintf(viewer, "  not yet available\n"));
      PetscCall(PetscViewerASCIIPopTab(viewer));
    }
    PetscCall(PetscViewerASCIIPrintf(viewer, "KSP solver for S = A11 - A10 inv(A00) A01\n"));
    PetscCall(PetscViewerASCIIPushTab(viewer));
    if (jac->kspschur) {
      PetscCall(KSPView(jac->kspschur, viewer));
    } else {
      PetscCall(PetscViewerASCIIPrintf(viewer, "  not yet available\n"));
    }
    PetscCall(PetscViewerASCIIPopTab(viewer));
    PetscCall(PetscViewerASCIIPopTab(viewer));
  } else if (isdraw && jac->head) {
    PetscDraw draw;
    PetscReal x, y, w, wd, h;
    PetscInt  cnt = 2;
    char      str[32];

    PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
    PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
    if (jac->kspupper != jac->head->ksp) cnt++;
    w  = 2 * PetscMin(1.0 - x, x);
    wd = w / (cnt + 1);

    PetscCall(PetscSNPrintf(str, 32, "Schur fact. %s", PCFieldSplitSchurFactTypes[jac->schurfactorization]));
    PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h));
    y -= h;
    if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) {
      PetscCall(PetscSNPrintf(str, 32, "Prec. for Schur from %s", PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]));
    } else {
      PetscCall(PetscSNPrintf(str, 32, "Prec. for Schur from %s", PCFieldSplitSchurPreTypes[jac->schurpre]));
    }
    PetscCall(PetscDrawStringBoxed(draw, x + wd * (cnt - 1) / 2.0, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h));
    y -= h;
    x = x - wd * (cnt - 1) / 2.0;

    PetscCall(PetscDrawPushCurrentPoint(draw, x, y));
    PetscCall(KSPView(jac->head->ksp, viewer));
    PetscCall(PetscDrawPopCurrentPoint(draw));
    if (jac->kspupper != jac->head->ksp) {
      x += wd;
      PetscCall(PetscDrawPushCurrentPoint(draw, x, y));
      PetscCall(KSPView(jac->kspupper, viewer));
      PetscCall(PetscDrawPopCurrentPoint(draw));
    }
    x += wd;
    PetscCall(PetscDrawPushCurrentPoint(draw, x, y));
    PetscCall(KSPView(jac->kspschur, viewer));
    PetscCall(PetscDrawPopCurrentPoint(draw));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCView_FieldSplit_GKB(PC pc, PetscViewer viewer)
{
  PC_FieldSplit    *jac = (PC_FieldSplit *)pc->data;
  PetscBool         isascii, isdraw;
  PetscInt          i, j;
  PC_FieldSplitLink ilink = jac->head;

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
  if (isascii) {
    if (jac->bs > 0) {
      PetscCall(PetscViewerASCIIPrintf(viewer, "  FieldSplit with %s composition: total splits = %" PetscInt_FMT ", blocksize = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits, jac->bs));
    } else {
      PetscCall(PetscViewerASCIIPrintf(viewer, "  FieldSplit with %s composition: total splits = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits));
    }
    if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, "  using Amat (not Pmat) as operator for blocks\n"));
    if (jac->diag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, "  using Amat (not Pmat) as operator for diagonal blocks\n"));
    if (jac->offdiag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, "  using Amat (not Pmat) as operator for off-diagonal blocks\n"));

    PetscCall(PetscViewerASCIIPrintf(viewer, "  Stopping tolerance=%.1e, delay in error estimate=%" PetscInt_FMT ", maximum iterations=%" PetscInt_FMT "\n", (double)jac->gkbtol, jac->gkbdelay, jac->gkbmaxit));
    PetscCall(PetscViewerASCIIPrintf(viewer, "  Solver info for H = A00 + nu*A01*A01' matrix:\n"));
    PetscCall(PetscViewerASCIIPushTab(viewer));

    if (ilink->fields) {
      PetscCall(PetscViewerASCIIPrintf(viewer, "Split number 0 Fields "));
      PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
      for (j = 0; j < ilink->nfields; j++) {
        if (j > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ","));
        PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j]));
      }
      PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
      PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
    } else {
      PetscCall(PetscViewerASCIIPrintf(viewer, "Split number 0 Defined by IS\n"));
    }
    PetscCall(KSPView(ilink->ksp, viewer));

    PetscCall(PetscViewerASCIIPopTab(viewer));
  }

  if (isdraw) {
    PetscDraw draw;
    PetscReal x, y, w, wd;

    PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
    PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
    w  = 2 * PetscMin(1.0 - x, x);
    wd = w / (jac->nsplits + 1);
    x  = x - wd * (jac->nsplits - 1) / 2.0;
    for (i = 0; i < jac->nsplits; i++) {
      PetscCall(PetscDrawPushCurrentPoint(draw, x, y));
      PetscCall(KSPView(ilink->ksp, viewer));
      PetscCall(PetscDrawPopCurrentPoint(draw));
      x += wd;
      ilink = ilink->next;
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/* Precondition: jac->bs is set to a meaningful value or MATNEST */
static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
  PetscInt       bs, i, nfields, *ifields, nfields_col, *ifields_col;
  PetscBool      flg, flg_col, mnest;
  char           optionname[128], splitname[8], optionname_col[128];

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)pc->mat, MATNEST, &mnest));
  if (mnest) {
    PetscCall(MatNestGetSize(pc->pmat, &bs, NULL));
  } else {
    bs = jac->bs;
  }
  PetscCall(PetscMalloc2(bs, &ifields, bs, &ifields_col));
  for (i = 0, flg = PETSC_TRUE;; i++) {
    PetscCall(PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i));
    PetscCall(PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%" PetscInt_FMT "_fields", i));
    PetscCall(PetscSNPrintf(optionname_col, sizeof(optionname_col), "-pc_fieldsplit_%" PetscInt_FMT "_fields_col", i));
    nfields     = bs;
    nfields_col = bs;
    PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg));
    PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname_col, ifields_col, &nfields_col, &flg_col));
    if (!flg) break;
    else if (flg && !flg_col) {
      PetscCheck(nfields, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot list zero fields");
      PetscCall(PCFieldSplitSetFields(pc, splitname, nfields, ifields, ifields));
    } else {
      PetscCheck(nfields && nfields_col, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot list zero fields");
      PetscCheck(nfields == nfields_col, PETSC_COMM_SELF, PETSC_ERR_USER, "Number of row and column fields must match");
      PetscCall(PCFieldSplitSetFields(pc, splitname, nfields, ifields, ifields_col));
    }
  }
  if (i > 0) {
    /* Makes command-line setting of splits take precedence over setting them in code.
       Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would
       create new splits, which would probably not be what the user wanted. */
    jac->splitdefined = PETSC_TRUE;
  }
  PetscCall(PetscFree2(ifields, ifields_col));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
{
  PC_FieldSplit    *jac                = (PC_FieldSplit *)pc->data;
  PC_FieldSplitLink ilink              = jac->head;
  PetscBool         fieldsplit_default = PETSC_FALSE, coupling = PETSC_FALSE;
  PetscInt          i;

  PetscFunctionBegin;
  /*
   Kinda messy, but at least this now uses DMCreateFieldDecomposition().
   Should probably be rewritten.
   */
  if (!ilink) {
    PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_detect_coupling", &coupling, NULL));
    if (pc->dm && jac->dm_splits && !jac->detect && !coupling) {
      PetscInt  numFields, f, i, j;
      char    **fieldNames;
      IS       *fields;
      DM       *dms;
      DM        subdm[128];
      PetscBool flg;

      PetscCall(DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms));
      /* Allow the user to prescribe the splits */
      for (i = 0, flg = PETSC_TRUE;; i++) {
        PetscInt ifields[128];
        IS       compField;
        char     optionname[128], splitname[8];
        PetscInt nfields = numFields;

        PetscCall(PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%" PetscInt_FMT "_fields", i));
        PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg));
        if (!flg) break;
        PetscCheck(numFields <= 128, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot currently support %" PetscInt_FMT " > 128 fields", numFields);
        PetscCall(DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]));
        if (nfields == 1) {
          PetscCall(PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField));
        } else {
          PetscCall(PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i));
          PetscCall(PCFieldSplitSetIS(pc, splitname, compField));
        }
        PetscCall(ISDestroy(&compField));
        for (j = 0; j < nfields; ++j) {
          f = ifields[j];
          PetscCall(PetscFree(fieldNames[f]));
          PetscCall(ISDestroy(&fields[f]));
        }
      }
      if (i == 0) {
        for (f = 0; f < numFields; ++f) {
          PetscCall(PCFieldSplitSetIS(pc, fieldNames[f], fields[f]));
          PetscCall(PetscFree(fieldNames[f]));
          PetscCall(ISDestroy(&fields[f]));
        }
      } else {
        for (j = 0; j < numFields; j++) PetscCall(DMDestroy(dms + j));
        PetscCall(PetscFree(dms));
        PetscCall(PetscMalloc1(i, &dms));
        for (j = 0; j < i; ++j) dms[j] = subdm[j];
      }
      PetscCall(PetscFree(fieldNames));
      PetscCall(PetscFree(fields));
      if (dms) {
        PetscCall(PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n"));
        for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) {
          const char *prefix;
          PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ilink->ksp, &prefix));
          PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dms[i], prefix));
          PetscCall(KSPSetDM(ilink->ksp, dms[i]));
          PetscCall(KSPSetDMActive(ilink->ksp, PETSC_FALSE));
          PetscCall(PetscObjectIncrementTabLevel((PetscObject)dms[i], (PetscObject)ilink->ksp, 0));
          PetscCall(DMDestroy(&dms[i]));
        }
        PetscCall(PetscFree(dms));
      }
    } else {
      if (jac->bs <= 0) {
        if (pc->pmat) PetscCall(MatGetBlockSize(pc->pmat, &jac->bs));
        else jac->bs = 1;
      }

      if (jac->detect) {
        IS       zerodiags, rest;
        PetscInt nmin, nmax;

        PetscCall(MatGetOwnershipRange(pc->mat, &nmin, &nmax));
        if (jac->diag_use_amat) {
          PetscCall(MatFindZeroDiagonals(pc->mat, &zerodiags));
        } else {
          PetscCall(MatFindZeroDiagonals(pc->pmat, &zerodiags));
        }
        PetscCall(ISComplement(zerodiags, nmin, nmax, &rest));
        PetscCall(PCFieldSplitSetIS(pc, "0", rest));
        PetscCall(PCFieldSplitSetIS(pc, "1", zerodiags));
        PetscCall(ISDestroy(&zerodiags));
        PetscCall(ISDestroy(&rest));
      } else if (coupling) {
        IS       coupling, rest;
        PetscInt nmin, nmax;

        PetscCall(MatGetOwnershipRange(pc->mat, &nmin, &nmax));
        if (jac->offdiag_use_amat) {
          PetscCall(MatFindOffBlockDiagonalEntries(pc->mat, &coupling));
        } else {
          PetscCall(MatFindOffBlockDiagonalEntries(pc->pmat, &coupling));
        }
        PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc->mat), nmax - nmin, nmin, 1, &rest));
        PetscCall(ISSetIdentity(rest));
        PetscCall(PCFieldSplitSetIS(pc, "0", rest));
        PetscCall(PCFieldSplitSetIS(pc, "1", coupling));
        PetscCall(ISDestroy(&coupling));
        PetscCall(ISDestroy(&rest));
      } else {
        PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_default", &fieldsplit_default, NULL));
        if (!fieldsplit_default) {
          /* Allow user to set fields from command line,  if bs was known at the time of PCSetFromOptions_FieldSplit()
           then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */
          PetscCall(PCFieldSplitSetRuntimeSplits_Private(pc));
          if (jac->splitdefined) PetscCall(PetscInfo(pc, "Splits defined using the options database\n"));
        }
        if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) {
          Mat       M = pc->pmat;
          PetscBool isnest;
          PetscInt  nf;

          PetscCall(PetscInfo(pc, "Using default splitting of fields\n"));
          PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &isnest));
          if (!isnest) {
            M = pc->mat;
            PetscCall(PetscObjectTypeCompare((PetscObject)pc->mat, MATNEST, &isnest));
          }
          if (!isnest) nf = jac->bs;
          else PetscCall(MatNestGetSize(M, &nf, NULL));
          for (i = 0; i < nf; i++) {
            char splitname[8];

            PetscCall(PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i));
            PetscCall(PCFieldSplitSetFields(pc, splitname, 1, &i, &i));
          }
          jac->defaultsplit = PETSC_TRUE;
        }
      }
    }
  } else if (jac->nsplits == 1) {
    IS       is2;
    PetscInt nmin, nmax;

    PetscCheck(ilink->is, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Must provide at least two sets of fields to PCFieldSplit()");
    PetscCall(MatGetOwnershipRange(pc->mat, &nmin, &nmax));
    PetscCall(ISComplement(ilink->is, nmin, nmax, &is2));
    PetscCall(PCFieldSplitSetIS(pc, "1", is2));
    PetscCall(ISDestroy(&is2));
  }

  PetscCheck(jac->nsplits >= 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Unhandled case, must have at least two fields, not %" PetscInt_FMT, jac->nsplits);
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode MatGolubKahanComputeExplicitOperator(Mat A, Mat B, Mat C, Mat *H, PetscReal gkbnu)
{
  Mat       BT, T;
  PetscReal nrmT, nrmB;

  PetscFunctionBegin;
  PetscCall(MatHermitianTranspose(C, MAT_INITIAL_MATRIX, &T)); /* Test if augmented matrix is symmetric */
  PetscCall(MatAXPY(T, -1.0, B, DIFFERENT_NONZERO_PATTERN));
  PetscCall(MatNorm(T, NORM_1, &nrmT));
  PetscCall(MatNorm(B, NORM_1, &nrmB));
  PetscCheck(nrmB <= 0 || nrmT / nrmB < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix is not symmetric/Hermitian, GKB is not applicable.");

  /* Compute augmented Lagrangian matrix H = A00 + nu*A01*A01'. This corresponds to */
  /* setting N := 1/nu*I in [Ar13].                                                 */
  PetscCall(MatHermitianTranspose(B, MAT_INITIAL_MATRIX, &BT));
  PetscCall(MatMatMult(B, BT, MAT_INITIAL_MATRIX, PETSC_CURRENT, H)); /* H = A01*A01'          */
  PetscCall(MatAYPX(*H, gkbnu, A, DIFFERENT_NONZERO_PATTERN));        /* H = A00 + nu*A01*A01' */

  PetscCall(MatDestroy(&BT));
  PetscCall(MatDestroy(&T));
  PetscFunctionReturn(PETSC_SUCCESS);
}

PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions, const char pre[], const char name[], const char *option[], const char *value[], PetscBool *flg);

static PetscErrorCode PCSetUp_FieldSplit(PC pc)
{
  PC_FieldSplit    *jac = (PC_FieldSplit *)pc->data;
  PC_FieldSplitLink ilink;
  PetscInt          i, nsplit;
  PetscBool         matnest = PETSC_FALSE;

  PetscFunctionBegin;
  pc->failedreason = PC_NOERROR;
  PetscCall(PCFieldSplitSetDefaults(pc));
  nsplit = jac->nsplits;
  ilink  = jac->head;
  if (pc->pmat) PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &matnest));

  /* get the matrices for each split */
  if (!jac->issetup) {
    PetscInt rstart, rend, nslots, bs;

    jac->issetup = PETSC_TRUE;

    /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */
    if (jac->defaultsplit || !ilink->is) {
      if (jac->bs <= 0) jac->bs = nsplit;
    }

    /*  MatCreateSubMatrix() for [S]BAIJ matrices can only work if the indices include entire blocks of the matrix */
    PetscCall(MatGetBlockSize(pc->pmat, &bs));
    if (bs > 1 && (jac->bs <= bs || jac->bs % bs)) {
      PetscBool blk;

      PetscCall(PetscObjectTypeCompareAny((PetscObject)pc->pmat, &blk, MATBAIJ, MATSBAIJ, MATSEQBAIJ, MATSEQSBAIJ, MATMPIBAIJ, MATMPISBAIJ, NULL));
      PetscCheck(!blk, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Cannot use MATBAIJ with PCFIELDSPLIT and currently set matrix and PC blocksizes");
    }

    if (!matnest) { /* use the matrix blocksize and stride IS to determine the index sets that define the submatrices */
      bs = jac->bs;
      PetscCall(MatGetOwnershipRange(pc->pmat, &rstart, &rend));
      nslots = (rend - rstart) / bs;
      for (i = 0; i < nsplit; i++) {
        if (jac->defaultsplit) {
          PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + i, nsplit, &ilink->is));
          PetscCall(PetscObjectReference((PetscObject)ilink->is));
          ilink->is_col = ilink->is;
        } else if (!ilink->is) {
          PetscBool same_fields = PETSC_TRUE;

          for (PetscInt k = 0; k < ilink->nfields; k++) {
            if (ilink->fields[k] != ilink->fields_col[k]) same_fields = PETSC_FALSE;
          }

          if (ilink->nfields > 1) {
            PetscInt *ii, *jj, j, k, nfields = ilink->nfields, *fields = ilink->fields, *fields_col = ilink->fields_col;

            PetscCall(PetscMalloc1(ilink->nfields * nslots, &ii));
            if (!same_fields) PetscCall(PetscMalloc1(ilink->nfields * nslots, &jj));
            for (j = 0; j < nslots; j++) {
              for (k = 0; k < nfields; k++) {
                ii[nfields * j + k] = rstart + bs * j + fields[k];
                if (!same_fields) jj[nfields * j + k] = rstart + bs * j + fields_col[k];
              }
            }
            PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), nslots * nfields, ii, PETSC_OWN_POINTER, &ilink->is));
            if (!same_fields) PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), nslots * nfields, jj, PETSC_OWN_POINTER, &ilink->is_col));
            else {
              PetscCall(PetscObjectReference((PetscObject)ilink->is));
              ilink->is_col = ilink->is;
            }
            PetscCall(ISSetBlockSize(ilink->is, nfields));
            PetscCall(ISSetBlockSize(ilink->is_col, nfields));
          } else {
            PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + ilink->fields[0], bs, &ilink->is));
            if (!same_fields) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + ilink->fields_col[0], bs, &ilink->is_col));
            else {
              PetscCall(PetscObjectReference((PetscObject)ilink->is));
              ilink->is_col = ilink->is;
            }
          }
        }
        ilink = ilink->next;
      }
    } else { /* use the IS that define the MATNEST to determine the index sets that define the submatrices */
      IS      *rowis, *colis, *ises = NULL;
      PetscInt mis, nis;

      PetscCall(MatNestGetSize(pc->pmat, &mis, &nis));
      PetscCall(PetscMalloc2(mis, &rowis, nis, &colis));
      PetscCall(MatNestGetISs(pc->pmat, rowis, colis));
      if (!jac->defaultsplit) PetscCall(PetscMalloc1(mis, &ises));

      for (i = 0; i < nsplit; i++) {
        if (jac->defaultsplit) {
          PetscCall(ISDuplicate(rowis[i], &ilink->is));
          PetscCall(PetscObjectReference((PetscObject)ilink->is));
          ilink->is_col = ilink->is;
        } else if (!ilink->is) {
          if (ilink->nfields > 1) {
            for (PetscInt j = 0; j < ilink->nfields; j++) ises[j] = rowis[ilink->fields[j]];
            PetscCall(ISConcatenate(PetscObjectComm((PetscObject)pc), ilink->nfields, ises, &ilink->is));
          } else {
            PetscCall(ISDuplicate(rowis[ilink->fields[0]], &ilink->is));
          }
          PetscCall(PetscObjectReference((PetscObject)ilink->is));
          ilink->is_col = ilink->is;
        }
        ilink = ilink->next;
      }
      PetscCall(PetscFree2(rowis, colis));
      PetscCall(PetscFree(ises));
    }
  }

  ilink = jac->head;
  if (!jac->pmat) {
    Vec xtmp;

    PetscCall(MatCreateVecs(pc->pmat, &xtmp, NULL));
    PetscCall(PetscMalloc1(nsplit, &jac->pmat));
    PetscCall(PetscMalloc2(nsplit, &jac->x, nsplit, &jac->y));
    for (i = 0; i < nsplit; i++) {
      MatNullSpace sp;

      /* Check for matrix attached to IS */
      PetscCall(PetscObjectQuery((PetscObject)ilink->is, "pmat", (PetscObject *)&jac->pmat[i]));
      if (jac->pmat[i]) {
        PetscCall(PetscObjectReference((PetscObject)jac->pmat[i]));
        if (jac->type == PC_COMPOSITE_SCHUR) {
          jac->schur_user = jac->pmat[i];

          PetscCall(PetscObjectReference((PetscObject)jac->schur_user));
        }
      } else {
        const char *prefix;
        PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ilink->is_col, MAT_INITIAL_MATRIX, &jac->pmat[i]));
        PetscCall(MatGetOptionsPrefix(jac->pmat[i], &prefix));
        if (!prefix) {
          PetscCall(KSPGetOptionsPrefix(ilink->ksp, &prefix));
          PetscCall(MatSetOptionsPrefix(jac->pmat[i], prefix));
        }
        PetscCall(MatSetFromOptions(jac->pmat[i]));
        PetscCall(MatViewFromOptions(jac->pmat[i], NULL, "-mat_view"));
      }
      /* create work vectors for each split */
      PetscCall(MatCreateVecs(jac->pmat[i], &jac->x[i], &jac->y[i]));
      ilink->x = jac->x[i];
      ilink->y = jac->y[i];
      ilink->z = NULL;
      /* compute scatter contexts needed by multiplicative versions and non-default splits */
      PetscCall(VecScatterCreate(xtmp, ilink->is, jac->x[i], NULL, &ilink->sctx));
      PetscCall(PetscObjectQuery((PetscObject)ilink->is, "nearnullspace", (PetscObject *)&sp));
      if (sp) PetscCall(MatSetNearNullSpace(jac->pmat[i], sp));
      ilink = ilink->next;
    }
    PetscCall(VecDestroy(&xtmp));
  } else {
    MatReuse      scall;
    MatNullSpace *nullsp = NULL;

    if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
      PetscCall(MatGetNullSpaces(nsplit, jac->pmat, &nullsp));
      for (i = 0; i < nsplit; i++) PetscCall(MatDestroy(&jac->pmat[i]));
      scall = MAT_INITIAL_MATRIX;
    } else scall = MAT_REUSE_MATRIX;

    for (i = 0; i < nsplit; i++) {
      Mat pmat;

      /* Check for matrix attached to IS */
      PetscCall(PetscObjectQuery((PetscObject)ilink->is, "pmat", (PetscObject *)&pmat));
      if (!pmat) PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ilink->is_col, scall, &jac->pmat[i]));
      ilink = ilink->next;
    }
    if (nullsp) PetscCall(MatRestoreNullSpaces(nsplit, jac->pmat, &nullsp));
  }
  if (jac->diag_use_amat) {
    ilink = jac->head;
    if (!jac->mat) {
      PetscCall(PetscMalloc1(nsplit, &jac->mat));
      for (i = 0; i < nsplit; i++) {
        PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ilink->is_col, MAT_INITIAL_MATRIX, &jac->mat[i]));
        ilink = ilink->next;
      }
    } else {
      MatReuse      scall;
      MatNullSpace *nullsp = NULL;

      if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
        PetscCall(MatGetNullSpaces(nsplit, jac->mat, &nullsp));
        for (i = 0; i < nsplit; i++) PetscCall(MatDestroy(&jac->mat[i]));
        scall = MAT_INITIAL_MATRIX;
      } else scall = MAT_REUSE_MATRIX;

      for (i = 0; i < nsplit; i++) {
        PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ilink->is_col, scall, &jac->mat[i]));
        ilink = ilink->next;
      }
      if (nullsp) PetscCall(MatRestoreNullSpaces(nsplit, jac->mat, &nullsp));
    }
  } else {
    jac->mat = jac->pmat;
  }

  /* Check for null space attached to IS */
  ilink = jac->head;
  for (i = 0; i < nsplit; i++) {
    MatNullSpace sp;

    PetscCall(PetscObjectQuery((PetscObject)ilink->is, "nullspace", (PetscObject *)&sp));
    if (sp) PetscCall(MatSetNullSpace(jac->mat[i], sp));
    ilink = ilink->next;
  }

  if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR && jac->type != PC_COMPOSITE_GKB) {
    /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
    /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
    ilink = jac->head;
    if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) {
      /* special case need where Afield[0] is not needed and only certain columns of Afield[1] are needed since update is only on those rows of the solution */
      if (!jac->Afield) {
        PetscCall(PetscCalloc1(nsplit, &jac->Afield));
        if (jac->offdiag_use_amat) {
          PetscCall(MatCreateSubMatrix(pc->mat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->Afield[1]));
        } else {
          PetscCall(MatCreateSubMatrix(pc->pmat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->Afield[1]));
        }
      } else {
        MatReuse scall;

        if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
          PetscCall(MatDestroy(&jac->Afield[1]));
          scall = MAT_INITIAL_MATRIX;
        } else scall = MAT_REUSE_MATRIX;

        if (jac->offdiag_use_amat) {
          PetscCall(MatCreateSubMatrix(pc->mat, ilink->next->is, ilink->is, scall, &jac->Afield[1]));
        } else {
          PetscCall(MatCreateSubMatrix(pc->pmat, ilink->next->is, ilink->is, scall, &jac->Afield[1]));
        }
      }
    } else {
      if (!jac->Afield) {
        PetscCall(PetscMalloc1(nsplit, &jac->Afield));
        for (i = 0; i < nsplit; i++) {
          if (jac->offdiag_use_amat) {
            PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, NULL, MAT_INITIAL_MATRIX, &jac->Afield[i]));
          } else {
            PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, NULL, MAT_INITIAL_MATRIX, &jac->Afield[i]));
          }
          ilink = ilink->next;
        }
      } else {
        MatReuse scall;
        if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
          for (i = 0; i < nsplit; i++) PetscCall(MatDestroy(&jac->Afield[i]));
          scall = MAT_INITIAL_MATRIX;
        } else scall = MAT_REUSE_MATRIX;

        for (i = 0; i < nsplit; i++) {
          if (jac->offdiag_use_amat) {
            PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, NULL, scall, &jac->Afield[i]));
          } else {
            PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, NULL, scall, &jac->Afield[i]));
          }
          ilink = ilink->next;
        }
      }
    }
  }

  if (jac->type == PC_COMPOSITE_SCHUR) {
    PetscBool   isset, isspd = PETSC_FALSE, issym = PETSC_FALSE, flg;
    char        lscname[256];
    PetscObject LSC_L;

    PetscCheck(nsplit == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "To use Schur complement preconditioner you must have exactly 2 fields");

    /* If pc->mat is SPD, don't scale by -1 the Schur complement */
    PetscCall(MatIsSPDKnown(pc->pmat, &isset, &isspd));
    if (jac->schurscale == (PetscScalar)-1.0) jac->schurscale = (isset && isspd) ? 1.0 : -1.0;
    PetscCall(MatIsSymmetricKnown(pc->pmat, &isset, &issym));

    PetscCall(PetscObjectTypeCompareAny(jac->offdiag_use_amat ? (PetscObject)pc->mat : (PetscObject)pc->pmat, &flg, MATSEQSBAIJ, MATMPISBAIJ, ""));

    if (jac->schur) {
      KSP      kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper;
      MatReuse scall;

      if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
        scall = MAT_INITIAL_MATRIX;
        PetscCall(MatDestroy(&jac->B));
        PetscCall(MatDestroy(&jac->C));
      } else scall = MAT_REUSE_MATRIX;

      PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner));
      ilink = jac->head;
      PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->is, ilink->next->is, scall, &jac->B));
      if (!flg) PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->next->is, ilink->is, scall, &jac->C));
      else {
        PetscCall(MatIsHermitianKnown(jac->offdiag_use_amat ? pc->mat : pc->pmat, &isset, &flg));
        if (isset && flg) PetscCall(MatCreateHermitianTranspose(jac->B, &jac->C));
        else PetscCall(MatCreateTranspose(jac->B, &jac->C));
      }
      ilink = ilink->next;
      PetscCall(MatSchurComplementUpdateSubMatrices(jac->schur, jac->mat[0], jac->pmat[0], jac->B, jac->C, jac->mat[1]));
      if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
        PetscCall(MatDestroy(&jac->schurp));
        PetscCall(MatSchurComplementGetPmat(jac->schur, MAT_INITIAL_MATRIX, &jac->schurp));
      } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL && jac->kspupper != jac->head->ksp) {
        PetscCall(MatDestroy(&jac->schur_user));
        PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user));
      }
      if (kspA != kspInner) PetscCall(KSPSetOperators(kspA, jac->mat[0], jac->pmat[0]));
      if (kspUpper != kspA) PetscCall(KSPSetOperators(kspUpper, jac->mat[0], jac->pmat[0]));
      PetscCall(KSPSetOperators(jac->kspschur, jac->schur, FieldSplitSchurPre(jac)));
    } else {
      const char  *Dprefix;
      char         schurprefix[256], schurmatprefix[256];
      char         schurtestoption[256];
      MatNullSpace sp;
      KSP          kspt;

      /* extract the A01 and A10 matrices */
      ilink = jac->head;
      PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->is, ilink->next->is, MAT_INITIAL_MATRIX, &jac->B));
      if (!flg) PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->C));
      else {
        PetscCall(MatIsHermitianKnown(jac->offdiag_use_amat ? pc->mat : pc->pmat, &isset, &flg));
        if (isset && flg) PetscCall(MatCreateHermitianTranspose(jac->B, &jac->C));
        else PetscCall(MatCreateTranspose(jac->B, &jac->C));
      }
      ilink = ilink->next;
      /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
      PetscCall(MatCreate(((PetscObject)jac->mat[0])->comm, &jac->schur));
      PetscCall(MatSetType(jac->schur, MATSCHURCOMPLEMENT));
      PetscCall(MatSchurComplementSetSubMatrices(jac->schur, jac->mat[0], jac->pmat[0], jac->B, jac->C, jac->mat[1]));
      PetscCall(PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
      PetscCall(MatSetOptionsPrefix(jac->schur, schurmatprefix));
      PetscCall(MatSchurComplementGetKSP(jac->schur, &kspt));
      PetscCall(KSPSetOptionsPrefix(kspt, schurmatprefix));

      /* Note: this is not true in general */
      PetscCall(MatGetNullSpace(jac->mat[1], &sp));
      if (sp) PetscCall(MatSetNullSpace(jac->schur, sp));

      PetscCall(PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname));
      PetscCall(PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, schurtestoption, NULL, NULL, &flg));
      if (flg) {
        DM  dmInner;
        KSP kspInner;
        PC  pcInner;

        PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner));
        PetscCall(KSPReset(kspInner));
        PetscCall(KSPSetOperators(kspInner, jac->mat[0], jac->pmat[0]));
        PetscCall(PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
        /* Indent this deeper to emphasize the "inner" nature of this solver. */
        PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject)pc, 2));
        PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject)pc, 2));
        PetscCall(KSPSetOptionsPrefix(kspInner, schurprefix));

        /* Set DM for new solver */
        PetscCall(KSPGetDM(jac->head->ksp, &dmInner));
        PetscCall(KSPSetDM(kspInner, dmInner));
        PetscCall(KSPSetDMActive(kspInner, PETSC_FALSE));

        /* Defaults to PCKSP as preconditioner */
        PetscCall(KSPGetPC(kspInner, &pcInner));
        PetscCall(PCSetType(pcInner, PCKSP));
        PetscCall(PCKSPSetKSP(pcInner, jac->head->ksp));
      } else {
        /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
          * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
          * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
          * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
          * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
          * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
        PetscCall(KSPSetType(jac->head->ksp, KSPGMRES));
        PetscCall(MatSchurComplementSetKSP(jac->schur, jac->head->ksp));
      }
      PetscCall(KSPSetOperators(jac->head->ksp, jac->mat[0], jac->pmat[0]));
      PetscCall(KSPSetFromOptions(jac->head->ksp));
      PetscCall(MatSetFromOptions(jac->schur));

      PetscCall(PetscObjectTypeCompare((PetscObject)jac->schur, MATSCHURCOMPLEMENT, &flg));
      if (flg) { /* Need to do this otherwise PCSetUp_KSP will overwrite the amat of jac->head->ksp */
        KSP kspInner;
        PC  pcInner;

        PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner));
        PetscCall(KSPGetPC(kspInner, &pcInner));
        PetscCall(PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg));
        if (flg) {
          KSP ksp;

          PetscCall(PCKSPGetKSP(pcInner, &ksp));
          if (ksp == jac->head->ksp) PetscCall(PCSetUseAmat(pcInner, PETSC_TRUE));
        }
      }
      PetscCall(PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname));
      PetscCall(PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, schurtestoption, NULL, NULL, &flg));
      if (flg) {
        DM dmInner;

        PetscCall(PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
        PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper));
        PetscCall(KSPSetNestLevel(jac->kspupper, pc->kspnestlevel));
        PetscCall(KSPSetErrorIfNotConverged(jac->kspupper, pc->erroriffailure));
        PetscCall(KSPSetOptionsPrefix(jac->kspupper, schurprefix));
        PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject)pc, 1));
        PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject)pc, 1));
        PetscCall(KSPGetDM(jac->head->ksp, &dmInner));
        PetscCall(KSPSetDM(jac->kspupper, dmInner));
        PetscCall(KSPSetDMActive(jac->kspupper, PETSC_FALSE));
        PetscCall(KSPSetFromOptions(jac->kspupper));
        PetscCall(KSPSetOperators(jac->kspupper, jac->mat[0], jac->pmat[0]));
        PetscCall(VecDuplicate(jac->head->x, &jac->head->z));
      } else {
        jac->kspupper = jac->head->ksp;
        PetscCall(PetscObjectReference((PetscObject)jac->head->ksp));
      }

      if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) PetscCall(MatSchurComplementGetPmat(jac->schur, MAT_INITIAL_MATRIX, &jac->schurp));
      PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspschur));
      PetscCall(KSPSetNestLevel(jac->kspschur, pc->kspnestlevel));
      PetscCall(KSPSetErrorIfNotConverged(jac->kspschur, pc->erroriffailure));
      PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspschur, (PetscObject)pc, 1));
      if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
        PC pcschur;
        PetscCall(KSPGetPC(jac->kspschur, &pcschur));
        PetscCall(PCSetType(pcschur, PCNONE));
        /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
      } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
        if (jac->schurfactorization != PC_FIELDSPLIT_SCHUR_FACT_FULL || jac->kspupper != jac->head->ksp) PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user));
      }
      PetscCall(KSPSetOperators(jac->kspschur, jac->schur, FieldSplitSchurPre(jac)));
      PetscCall(KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix));
      PetscCall(KSPSetOptionsPrefix(jac->kspschur, Dprefix));
      /* propagate DM */
      {
        DM sdm;
        PetscCall(KSPGetDM(jac->head->next->ksp, &sdm));
        if (sdm) {
          PetscCall(KSPSetDM(jac->kspschur, sdm));
          PetscCall(KSPSetDMActive(jac->kspschur, PETSC_FALSE));
        }
      }
      /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
      /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
      PetscCall(KSPSetFromOptions(jac->kspschur));
    }
    PetscCall(MatAssemblyBegin(jac->schur, MAT_FINAL_ASSEMBLY));
    PetscCall(MatAssemblyEnd(jac->schur, MAT_FINAL_ASSEMBLY));
    if (issym) PetscCall(MatSetOption(jac->schur, MAT_SYMMETRIC, PETSC_TRUE));
    if (isspd) PetscCall(MatSetOption(jac->schur, MAT_SPD, PETSC_TRUE));

    /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
    PetscCall(PetscSNPrintf(lscname, sizeof(lscname), "%s_LSC_L", ilink->splitname));
    PetscCall(PetscObjectQuery((PetscObject)pc->mat, lscname, &LSC_L));
    if (!LSC_L) PetscCall(PetscObjectQuery((PetscObject)pc->pmat, lscname, &LSC_L));
    if (LSC_L) PetscCall(PetscObjectCompose((PetscObject)jac->schur, "LSC_L", LSC_L));
    PetscCall(PetscSNPrintf(lscname, sizeof(lscname), "%s_LSC_Lp", ilink->splitname));
    PetscCall(PetscObjectQuery((PetscObject)pc->pmat, lscname, &LSC_L));
    if (!LSC_L) PetscCall(PetscObjectQuery((PetscObject)pc->mat, lscname, &LSC_L));
    if (LSC_L) PetscCall(PetscObjectCompose((PetscObject)jac->schur, "LSC_Lp", LSC_L));
  } else if (jac->type == PC_COMPOSITE_GKB) {
    PetscCheck(nsplit == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "To use GKB preconditioner you must have exactly 2 fields");
    ilink = jac->head;
    PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->is, ilink->next->is, MAT_INITIAL_MATRIX, &jac->B));
    /* Create work vectors for GKB algorithm */
    PetscCall(VecDuplicate(ilink->x, &jac->u));
    PetscCall(VecDuplicate(ilink->x, &jac->Hu));
    PetscCall(VecDuplicate(ilink->x, &jac->w2));
    PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->C));
    ilink = ilink->next;
    /* Create work vectors for GKB algorithm */
    PetscCall(VecDuplicate(ilink->x, &jac->v));
    PetscCall(VecDuplicate(ilink->x, &jac->d));
    PetscCall(VecDuplicate(ilink->x, &jac->w1));
    PetscCall(MatGolubKahanComputeExplicitOperator(jac->mat[0], jac->B, jac->C, &jac->H, jac->gkbnu));
    PetscCall(PetscCalloc1(jac->gkbdelay, &jac->vecz));

    ilink = jac->head;
    PetscCall(KSPSetOperators(ilink->ksp, jac->H, jac->H));
    if (!jac->suboptionsset) PetscCall(KSPSetFromOptions(ilink->ksp));
    /* Create gkb_monitor context */
    if (jac->gkbmonitor) {
      PetscInt tablevel;
      PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &jac->gkbviewer));
      PetscCall(PetscViewerSetType(jac->gkbviewer, PETSCVIEWERASCII));
      PetscCall(PetscObjectGetTabLevel((PetscObject)ilink->ksp, &tablevel));
      PetscCall(PetscViewerASCIISetTab(jac->gkbviewer, tablevel));
      PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)ilink->ksp, 1));
    }
  } else {
    /* set up the individual splits' PCs */
    i     = 0;
    ilink = jac->head;
    while (ilink) {
      PetscCall(KSPSetOperators(ilink->ksp, jac->mat[i], jac->pmat[i]));
      /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
      if (!jac->suboptionsset) PetscCall(KSPSetFromOptions(ilink->ksp));
      i++;
      ilink = ilink->next;
    }
  }

  /* Set coordinates to the sub PC objects whenever these are set */
  if (jac->coordinates_set) {
    PC pc_coords;
    if (jac->type == PC_COMPOSITE_SCHUR) {
      // Head is first block.
      PetscCall(KSPGetPC(jac->head->ksp, &pc_coords));
      PetscCall(PCSetCoordinates(pc_coords, jac->head->dim, jac->head->ndofs, jac->head->coords));
      // Second one is Schur block, but its KSP object is in kspschur.
      PetscCall(KSPGetPC(jac->kspschur, &pc_coords));
      PetscCall(PCSetCoordinates(pc_coords, jac->head->next->dim, jac->head->next->ndofs, jac->head->next->coords));
    } else if (jac->type == PC_COMPOSITE_GKB) {
      PetscCall(PetscInfo(pc, "Warning: Setting coordinates does nothing for the GKB Fieldpslit preconditioner\n"));
    } else {
      ilink = jac->head;
      while (ilink) {
        PetscCall(KSPGetPC(ilink->ksp, &pc_coords));
        PetscCall(PCSetCoordinates(pc_coords, ilink->dim, ilink->ndofs, ilink->coords));
        ilink = ilink->next;
      }
    }
  }

  jac->suboptionsset = PETSC_TRUE;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCSetUpOnBlocks_FieldSplit_Schur(PC pc)
{
  PC_FieldSplit    *jac    = (PC_FieldSplit *)pc->data;
  PC_FieldSplitLink ilinkA = jac->head;
  KSP               kspA = ilinkA->ksp, kspUpper = jac->kspupper;

  PetscFunctionBegin;
  if (jac->schurfactorization == PC_FIELDSPLIT_SCHUR_FACT_FULL && kspUpper != kspA) {
    PetscCall(KSPSetUp(kspUpper));
    PetscCall(KSPSetUpOnBlocks(kspUpper));
  }
  PetscCall(KSPSetUp(kspA));
  PetscCall(KSPSetUpOnBlocks(kspA));
  if (jac->schurpre != PC_FIELDSPLIT_SCHUR_PRE_FULL) {
    PetscCall(KSPSetUp(jac->kspschur));
    PetscCall(KSPSetUpOnBlocks(jac->kspschur));
  } else if (kspUpper == kspA && jac->schurfactorization == PC_FIELDSPLIT_SCHUR_FACT_FULL) {
    Mat          A;
    PetscInt     m, M, N;
    VecType      vtype;
    PetscMemType mtype;
    PetscScalar *array;

    PetscCall(MatGetSize(jac->B, &M, &N));
    PetscCall(MatGetLocalSize(jac->B, &m, NULL));
    PetscCall(MatGetVecType(jac->B, &vtype));
    PetscCall(VecGetArrayAndMemType(ilinkA->x, &array, &mtype));
    PetscCall(VecRestoreArrayAndMemType(ilinkA->x, &array));
    PetscCall(PetscObjectQuery((PetscObject)jac->schur, "AinvB", (PetscObject *)&A));
    if (!A) {
      if (PetscMemTypeHost(mtype) || (!PetscDefined(HAVE_CUDA) && !PetscDefined(HAVE_HIP))) PetscCall(PetscMalloc1(m * (N + 1), &array));
#if PetscDefined(HAVE_CUDA)
      else if (PetscMemTypeCUDA(mtype)) PetscCallCUDA(cudaMalloc((void **)&array, sizeof(PetscScalar) * m * (N + 1)));
#endif
#if PetscDefined(HAVE_HIP)
      else if (PetscMemTypeHIP(mtype)) PetscCallHIP(hipMalloc((void **)&array, sizeof(PetscScalar) * m * (N + 1)));
#endif
      PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)jac->schur), vtype, m, PETSC_DECIDE, M, N + 1, -1, array, &A)); // number of columns of the Schur complement plus one
      PetscCall(PetscObjectCompose((PetscObject)jac->schur, "AinvB", (PetscObject)A));
      PetscCall(MatDestroy(&A));
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCSetUpOnBlocks_FieldSplit(PC pc)
{
  PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
  PC_FieldSplitLink ilink = jac->head;

  PetscFunctionBegin;
  while (ilink) {
    PetscCall(KSPSetUp(ilink->ksp));
    PetscCall(KSPSetUpOnBlocks(ilink->ksp));
    ilink = ilink->next;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCSetUpOnBlocks_FieldSplit_GKB(PC pc)
{
  PC_FieldSplit    *jac    = (PC_FieldSplit *)pc->data;
  PC_FieldSplitLink ilinkA = jac->head;
  KSP               ksp    = ilinkA->ksp;

  PetscFunctionBegin;
  PetscCall(KSPSetUp(ksp));
  PetscCall(KSPSetUpOnBlocks(ksp));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCApply_FieldSplit_Schur(PC pc, Vec x, Vec y)
{
  PC_FieldSplit    *jac    = (PC_FieldSplit *)pc->data;
  PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
  KSP               kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;
  Mat               AinvB = NULL;
  PetscInt          N, P;

  PetscFunctionBegin;
  switch (jac->schurfactorization) {
  case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
    /* [A00 0; 0 -S], positive definite, suitable for MINRES */
    PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
    PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y));
    PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    PetscCall(VecScale(ilinkD->y, jac->schurscale));
    PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    break;
  case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
    /* [A00 0; A10 S], suitable for left preconditioning */
    PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
    PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    PetscCall(MatMult(jac->C, ilinkA->y, ilinkD->x));
    PetscCall(VecScale(ilinkD->x, -1.));
    PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
    PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
    PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y));
    PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    break;
  case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
    /* [A00 A01; 0 S], suitable for right preconditioning */
    PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y));
    PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    PetscCall(MatMult(jac->B, ilinkD->y, ilinkA->x));
    PetscCall(VecScale(ilinkA->x, -1.));
    PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
    PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
    PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
    PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    break;
  case PC_FIELDSPLIT_SCHUR_FACT_FULL:
    /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1] */
    PetscCall(MatGetSize(jac->B, NULL, &P));
    N = P;
    PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(PetscLogEventBegin(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->y, NULL));
    if (kspUpper == kspA) {
      PetscCall(PetscObjectQuery((PetscObject)jac->schur, "AinvB", (PetscObject *)&AinvB));
      if (AinvB) {
        PetscCall(MatGetSize(AinvB, NULL, &N));
        if (N > P) { // first time PCApply_FieldSplit_Schur() is called
          PetscMemType mtype;
          Vec          c = NULL;
          PetscScalar *array;
          PetscInt     m, M;

          PetscCall(MatGetSize(jac->B, &M, NULL));
          PetscCall(MatGetLocalSize(jac->B, &m, NULL));
          PetscCall(MatDenseGetArrayAndMemType(AinvB, &array, &mtype));
          if (PetscMemTypeHost(mtype) || (!PetscDefined(HAVE_CUDA) && !PetscDefined(HAVE_HIP))) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)jac->schur), 1, m, M, array + m * P, &c));
#if PetscDefined(HAVE_CUDA)
          else if (PetscMemTypeCUDA(mtype)) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)jac->schur), 1, m, M, array + m * P, &c));
#endif
#if PetscDefined(HAVE_HIP)
          else if (PetscMemTypeHIP(mtype)) PetscCall(VecCreateMPIHIPWithArray(PetscObjectComm((PetscObject)jac->schur), 1, m, M, array + m * P, &c));
#endif
          PetscCall(MatDenseRestoreArrayAndMemType(AinvB, &array));
          PetscCall(VecCopy(ilinkA->x, c));
          PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user));
          PetscCall(KSPSetOperators(jac->kspschur, jac->schur, jac->schur_user));
          PetscCall(VecCopy(c, ilinkA->y)); // retrieve the solution as the last column of the composed Mat
          PetscCall(VecDestroy(&c));
        }
      }
    }
    if (N == P) PetscCall(KSPSolve(kspLower, ilinkA->x, ilinkA->y));
    PetscCall(KSPCheckSolve(kspLower, pc, ilinkA->y));
    PetscCall(PetscLogEventEnd(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->y, NULL));
    PetscCall(MatMult(jac->C, ilinkA->y, ilinkD->x));
    PetscCall(VecScale(ilinkD->x, -1.0));
    PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
    PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));

    PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y));
    PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));

    if (kspUpper == kspA) {
      if (!AinvB) {
        PetscCall(MatMult(jac->B, ilinkD->y, ilinkA->y));
        PetscCall(VecAXPY(ilinkA->x, -1.0, ilinkA->y));
        PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
        PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
        PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
        PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
      } else PetscCall(MatMultAdd(AinvB, ilinkD->y, ilinkA->y, ilinkA->y));
    } else {
      PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
      PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
      PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
      PetscCall(MatMult(jac->B, ilinkD->y, ilinkA->x));
      PetscCall(PetscLogEventBegin(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->z, NULL));
      PetscCall(KSPSolve(kspUpper, ilinkA->x, ilinkA->z));
      PetscCall(KSPCheckSolve(kspUpper, pc, ilinkA->z));
      PetscCall(PetscLogEventEnd(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->z, NULL));
      PetscCall(VecAXPY(ilinkA->y, -1.0, ilinkA->z));
    }
    PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*
  PCFieldSplitCreateWorkMats_Private - Allocate per-field dense work matrices for multi-RHS

  Input Parameters:
+ pc - the PC context
- X  - matrix to copy column-layout from

  Notes:
  If matrices already exist with correct column count, they are reused.
  If column count changed, old matrices are destroyed and new ones created.
*/
static PetscErrorCode PCFieldSplitCreateWorkMats_Private(PC pc, Mat X)
{
  PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
  PC_FieldSplitLink ilink = jac->head;
  PetscInt          mx, Mx, my, My, N;

  PetscFunctionBegin;
  while (ilink) {
    /* check if reallocation needed (previous allocation with wrong column count) */
    if (ilink->X) {
      PetscCall(MatGetSize(ilink->X, NULL, &N));
      if (N != X->cmap->N) {
        PetscCall(MatDestroy(&ilink->X));
        PetscCall(MatDestroy(&ilink->Y));
      }
    }
    /* create if needed */
    if (!ilink->X) {
      VecType xtype, ytype;

      PetscCall(VecGetType(ilink->x, &xtype));
      PetscCall(VecGetType(ilink->y, &ytype));
      PetscCall(VecGetLocalSize(ilink->x, &mx));
      PetscCall(VecGetSize(ilink->x, &Mx));
      PetscCall(VecGetLocalSize(ilink->y, &my));
      PetscCall(VecGetSize(ilink->y, &My));
      /* use default lda */
      PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)pc), xtype, mx, X->cmap->n, Mx, X->cmap->N, -1, NULL, &ilink->X));
      PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)pc), ytype, my, X->cmap->n, My, X->cmap->N, -1, NULL, &ilink->Y));
    }
    ilink = ilink->next;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCApplyTranspose_FieldSplit_Schur(PC pc, Vec x, Vec y)
{
  PC_FieldSplit    *jac    = (PC_FieldSplit *)pc->data;
  PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
  KSP               kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;

  PetscFunctionBegin;
  switch (jac->schurfactorization) {
  case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
    /* [A00 0; 0 -S], positive definite, suitable for MINRES */
    PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
    PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
    PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    PetscCall(VecScale(ilinkD->y, jac->schurscale));
    PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    break;
  case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
    PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
    PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    PetscCall(MatMultTranspose(jac->B, ilinkA->y, ilinkD->x));
    PetscCall(VecScale(ilinkD->x, -1.));
    PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
    PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
    PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
    PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    break;
  case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
    PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
    PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->x));
    PetscCall(VecScale(ilinkA->x, -1.));
    PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
    PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
    PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
    PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    break;
  case PC_FIELDSPLIT_SCHUR_FACT_FULL:
    PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(PetscLogEventBegin(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->y, NULL));
    PetscCall(KSPSolveTranspose(kspUpper, ilinkA->x, ilinkA->y));
    PetscCall(KSPCheckSolve(kspUpper, pc, ilinkA->y));
    PetscCall(PetscLogEventEnd(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->y, NULL));
    PetscCall(MatMultTranspose(jac->B, ilinkA->y, ilinkD->x));
    PetscCall(VecScale(ilinkD->x, -1.0));
    PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
    PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));

    PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
    PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));

    if (kspLower == kspA) {
      PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->y));
      PetscCall(VecAXPY(ilinkA->x, -1.0, ilinkA->y));
      PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
      PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
      PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
      PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    } else {
      PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
      PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
      PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
      PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->x));
      PetscCall(PetscLogEventBegin(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->z, NULL));
      PetscCall(KSPSolveTranspose(kspLower, ilinkA->x, ilinkA->z));
      PetscCall(KSPCheckSolve(kspLower, pc, ilinkA->z));
      PetscCall(PetscLogEventEnd(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->z, NULL));
      PetscCall(VecAXPY(ilinkA->y, -1.0, ilinkA->z));
    }
    PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

#define FieldSplitSplitSolveAdd(ilink, xx, yy) \
  ((PetscErrorCode)(VecScatterBegin(ilink->sctx, xx, ilink->x, INSERT_VALUES, SCATTER_FORWARD) || VecScatterEnd(ilink->sctx, xx, ilink->x, INSERT_VALUES, SCATTER_FORWARD) || PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL) || \
                    KSPSolve(ilink->ksp, ilink->x, ilink->y) || KSPCheckSolve(ilink->ksp, pc, ilink->y) || PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL) || VecScatterBegin(ilink->sctx, ilink->y, yy, ADD_VALUES, SCATTER_REVERSE) || \
                    VecScatterEnd(ilink->sctx, ilink->y, yy, ADD_VALUES, SCATTER_REVERSE)))

static PetscErrorCode PCApply_FieldSplit(PC pc, Vec x, Vec y)
{
  PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
  PC_FieldSplitLink ilink = jac->head;
  PetscInt          cnt, bs;

  PetscFunctionBegin;
  if (jac->type == PC_COMPOSITE_ADDITIVE) {
    PetscBool matnest;

    PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &matnest));
    if (jac->defaultsplit && !matnest) {
      PetscCall(VecGetBlockSize(x, &bs));
      PetscCheck(jac->bs <= 0 || bs == jac->bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Blocksize of x vector %" PetscInt_FMT " does not match fieldsplit blocksize %" PetscInt_FMT, bs, jac->bs);
      PetscCall(VecGetBlockSize(y, &bs));
      PetscCheck(jac->bs <= 0 || bs == jac->bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Blocksize of y vector %" PetscInt_FMT " does not match fieldsplit blocksize %" PetscInt_FMT, bs, jac->bs);
      PetscCall(VecStrideGatherAll(x, jac->x, INSERT_VALUES));
      while (ilink) {
        PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
        PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
        PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
        PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
        ilink = ilink->next;
      }
      PetscCall(VecStrideScatterAll(jac->y, y, INSERT_VALUES));
    } else {
      PetscCall(VecSet(y, 0.0));
      while (ilink) {
        PetscCall(FieldSplitSplitSolveAdd(ilink, x, y));
        ilink = ilink->next;
      }
    }
  } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
    PetscCall(VecSet(y, 0.0));
    /* solve on first block for first block variables */
    PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
    PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
    PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));

    /* compute the residual only onto second block variables using first block variables */
    PetscCall(MatMult(jac->Afield[1], ilink->y, ilink->next->x));
    ilink = ilink->next;
    PetscCall(VecScale(ilink->x, -1.0));
    PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
    PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));

    /* solve on second block variables */
    PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
    PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
    PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
    PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
  } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
    if (!jac->w1) {
      PetscCall(VecDuplicate(x, &jac->w1));
      PetscCall(VecDuplicate(x, &jac->w2));
    }
    PetscCall(VecSet(y, 0.0));
    PetscCall(FieldSplitSplitSolveAdd(ilink, x, y));
    cnt = 1;
    while (ilink->next) {
      ilink = ilink->next;
      /* compute the residual only over the part of the vector needed */
      PetscCall(MatMult(jac->Afield[cnt++], y, ilink->x));
      PetscCall(VecScale(ilink->x, -1.0));
      PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
      PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
      PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
      PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
      PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
      PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
      PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
      PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
    }
    if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
      cnt -= 2;
      while (ilink->previous) {
        ilink = ilink->previous;
        /* compute the residual only over the part of the vector needed */
        PetscCall(MatMult(jac->Afield[cnt--], y, ilink->x));
        PetscCall(VecScale(ilink->x, -1.0));
        PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
        PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
        PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
        PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
        PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
        PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
        PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
        PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
      }
    }
  } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Unsupported or unknown composition %d", (int)jac->type);
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCMatApply_FieldSplit(PC pc, Mat X, Mat Y)
{
  PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
  PC_FieldSplitLink ilink = jac->head;
  PetscInt          cnt;

  PetscFunctionBegin;
  /* create working matrices with the correct number of columns */
  PetscCall(PCFieldSplitCreateWorkMats_Private(pc, X));
  if (jac->type == PC_COMPOSITE_ADDITIVE) {
    PetscCall(MatZeroEntries(Y));
    while (ilink) {
      PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, INSERT_VALUES, SCATTER_FORWARD));
      PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
      PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y));
      PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
      PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE));
      ilink = ilink->next;
    }
  } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
    PetscCall(MatZeroEntries(Y));
    PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
    PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y));
    PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
    PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE));

    /* compute the residual only onto second block variables using first block variables */
    PetscCall(MatMatMult(jac->Afield[1], ilink->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilink->next->X));
    ilink = ilink->next;
    PetscCall(MatScale(ilink->X, -1.0));
    PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, ADD_VALUES, SCATTER_FORWARD));

    /* solve on second block variables */
    PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
    PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y));
    PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
    PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE));
  } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
    /* general multiplicative with any number of splits */
    PetscCall(MatZeroEntries(Y));
    /* first split */
    PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, INSERT_VALUES, SCATTER_FORWARD));
    PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
    PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y));
    PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
    PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE));
    cnt = 1;
    /* forward sweep */
    while (ilink->next) {
      ilink = ilink->next;
      /* compute the residual only over the part of the vector needed */
      PetscCall(MatMatMult(jac->Afield[cnt++], Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilink->X));
      PetscCall(MatScale(ilink->X, -1.0));
      PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, ADD_VALUES, SCATTER_FORWARD));
      PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
      PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y));
      PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
      PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE));
    }
    /* backward sweep for symmetric multiplicative */
    if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
      cnt -= 2;
      while (ilink->previous) {
        ilink = ilink->previous;
        /* compute the residual only over the part of the vector needed */
        PetscCall(MatMatMult(jac->Afield[cnt--], Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilink->X));
        PetscCall(MatScale(ilink->X, -1.0));
        PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, ADD_VALUES, SCATTER_FORWARD));
        PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
        PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y));
        PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
        PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE));
      }
    }
  } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PCMatApply() not implemented for this fieldsplit type");
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCApply_FieldSplit_GKB(PC pc, Vec x, Vec y)
{
  PC_FieldSplit    *jac    = (PC_FieldSplit *)pc->data;
  PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
  KSP               ksp = ilinkA->ksp;
  Vec               u, v, Hu, d, work1, work2;
  PetscScalar       alpha, z, nrmz2, *vecz;
  PetscReal         lowbnd, nu, beta;
  PetscInt          j, iterGKB;

  PetscFunctionBegin;
  PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
  PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
  PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
  PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));

  u     = jac->u;
  v     = jac->v;
  Hu    = jac->Hu;
  d     = jac->d;
  work1 = jac->w1;
  work2 = jac->w2;
  vecz  = jac->vecz;

  /* Change RHS to comply with matrix regularization H = A + nu*B*B' */
  /* Add q = q + nu*B*b */
  if (jac->gkbnu) {
    nu = jac->gkbnu;
    PetscCall(VecScale(ilinkD->x, jac->gkbnu));
    PetscCall(MatMultAdd(jac->B, ilinkD->x, ilinkA->x, ilinkA->x)); /* q = q + nu*B*b */
  } else {
    /* Situation when no augmented Lagrangian is used. Then we set inner  */
    /* matrix N = I in [Ar13], and thus nu = 1.                           */
    nu = 1;
  }

  /* Transform rhs from [q,tilde{b}] to [0,b] */
  PetscCall(PetscLogEventBegin(ilinkA->event, ksp, ilinkA->x, ilinkA->y, NULL));
  PetscCall(KSPSolve(ksp, ilinkA->x, ilinkA->y));
  PetscCall(KSPCheckSolve(ksp, pc, ilinkA->y));
  PetscCall(PetscLogEventEnd(ilinkA->event, ksp, ilinkA->x, ilinkA->y, NULL));
  PetscCall(MatMultHermitianTranspose(jac->B, ilinkA->y, work1));
  PetscCall(VecAXPBY(work1, 1.0 / nu, -1.0, ilinkD->x)); /* c = b - B'*x        */

  /* First step of algorithm */
  PetscCall(VecNorm(work1, NORM_2, &beta)); /* beta = sqrt(nu*c'*c)*/
  KSPCheckDot(ksp, beta);
  beta = PetscSqrtReal(nu) * beta;
  PetscCall(VecAXPBY(v, nu / beta, 0.0, work1)); /* v = nu/beta *c      */
  PetscCall(MatMult(jac->B, v, work2));          /* u = H^{-1}*B*v      */
  PetscCall(PetscLogEventBegin(ilinkA->event, ksp, work2, u, NULL));
  PetscCall(KSPSolve(ksp, work2, u));
  PetscCall(KSPCheckSolve(ksp, pc, u));
  PetscCall(PetscLogEventEnd(ilinkA->event, ksp, work2, u, NULL));
  PetscCall(MatMult(jac->H, u, Hu)); /* alpha = u'*H*u      */
  PetscCall(VecDot(Hu, u, &alpha));
  KSPCheckDot(ksp, alpha);
  PetscCheck(PetscRealPart(alpha) > 0.0, PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "GKB preconditioner diverged, H is not positive definite");
  alpha = PetscSqrtReal(PetscAbsScalar(alpha));
  PetscCall(VecScale(u, 1.0 / alpha));
  PetscCall(VecAXPBY(d, 1.0 / alpha, 0.0, v)); /* v = nu/beta *c      */

  z       = beta / alpha;
  vecz[1] = z;

  /* Computation of first iterate x(1) and p(1) */
  PetscCall(VecAXPY(ilinkA->y, z, u));
  PetscCall(VecCopy(d, ilinkD->y));
  PetscCall(VecScale(ilinkD->y, -z));

  iterGKB = 1;
  lowbnd  = 2 * jac->gkbtol;
  if (jac->gkbmonitor) PetscCall(PetscViewerASCIIPrintf(jac->gkbviewer, "%3" PetscInt_FMT " GKB Lower bound estimate %14.12e\n", iterGKB, (double)lowbnd));

  while (iterGKB < jac->gkbmaxit && lowbnd > jac->gkbtol) {
    iterGKB += 1;
    PetscCall(MatMultHermitianTranspose(jac->B, u, work1)); /* v <- nu*(B'*u-alpha/nu*v) */
    PetscCall(VecAXPBY(v, nu, -alpha, work1));
    PetscCall(VecNorm(v, NORM_2, &beta)); /* beta = sqrt(nu)*v'*v      */
    beta = beta / PetscSqrtReal(nu);
    PetscCall(VecScale(v, 1.0 / beta));
    PetscCall(MatMult(jac->B, v, work2)); /* u <- H^{-1}*(B*v-beta*H*u) */
    PetscCall(MatMult(jac->H, u, Hu));
    PetscCall(VecAXPY(work2, -beta, Hu));
    PetscCall(PetscLogEventBegin(ilinkA->event, ksp, work2, u, NULL));
    PetscCall(KSPSolve(ksp, work2, u));
    PetscCall(KSPCheckSolve(ksp, pc, u));
    PetscCall(PetscLogEventEnd(ilinkA->event, ksp, work2, u, NULL));
    PetscCall(MatMult(jac->H, u, Hu)); /* alpha = u'*H*u            */
    PetscCall(VecDot(Hu, u, &alpha));
    KSPCheckDot(ksp, alpha);
    PetscCheck(PetscRealPart(alpha) > 0.0, PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "GKB preconditioner diverged, H is not positive definite");
    alpha = PetscSqrtReal(PetscAbsScalar(alpha));
    PetscCall(VecScale(u, 1.0 / alpha));

    z       = -beta / alpha * z; /* z <- beta/alpha*z     */
    vecz[0] = z;

    /* Computation of new iterate x(i+1) and p(i+1) */
    PetscCall(VecAXPBY(d, 1.0 / alpha, -beta / alpha, v)); /* d = (v-beta*d)/alpha */
    PetscCall(VecAXPY(ilinkA->y, z, u));                   /* r = r + z*u          */
    PetscCall(VecAXPY(ilinkD->y, -z, d));                  /* p = p - z*d          */
    PetscCall(MatMult(jac->H, ilinkA->y, Hu));             /* ||u||_H = u'*H*u     */
    PetscCall(VecDot(Hu, ilinkA->y, &nrmz2));

    /* Compute Lower Bound estimate */
    if (iterGKB > jac->gkbdelay) {
      lowbnd = 0.0;
      for (j = 0; j < jac->gkbdelay; j++) lowbnd += PetscAbsScalar(vecz[j] * vecz[j]);
      lowbnd = PetscSqrtReal(lowbnd / PetscAbsScalar(nrmz2));
    }

    for (j = 0; j < jac->gkbdelay - 1; j++) vecz[jac->gkbdelay - j - 1] = vecz[jac->gkbdelay - j - 2];
    if (jac->gkbmonitor) PetscCall(PetscViewerASCIIPrintf(jac->gkbviewer, "%3" PetscInt_FMT " GKB Lower bound estimate %14.12e\n", iterGKB, (double)lowbnd));
  }

  PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
  PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
  PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
  PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
  PetscFunctionReturn(PETSC_SUCCESS);
}

#define FieldSplitSplitSolveAddTranspose(ilink, xx, yy) \
  ((PetscErrorCode)(VecScatterBegin(ilink->sctx, xx, ilink->y, INSERT_VALUES, SCATTER_FORWARD) || VecScatterEnd(ilink->sctx, xx, ilink->y, INSERT_VALUES, SCATTER_FORWARD) || PetscLogEventBegin(ilink->event, ilink->ksp, ilink->y, ilink->x, NULL) || \
                    KSPSolveTranspose(ilink->ksp, ilink->y, ilink->x) || KSPCheckSolve(ilink->ksp, pc, ilink->x) || PetscLogEventEnd(ilink->event, ilink->ksp, ilink->y, ilink->x, NULL) || VecScatterBegin(ilink->sctx, ilink->x, yy, ADD_VALUES, SCATTER_REVERSE) || \
                    VecScatterEnd(ilink->sctx, ilink->x, yy, ADD_VALUES, SCATTER_REVERSE)))

static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc, Vec x, Vec y)
{
  PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
  PC_FieldSplitLink ilink = jac->head;
  PetscInt          bs;

  PetscFunctionBegin;
  if (jac->type == PC_COMPOSITE_ADDITIVE) {
    PetscBool matnest;

    PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &matnest));
    if (jac->defaultsplit && !matnest) {
      PetscCall(VecGetBlockSize(x, &bs));
      PetscCheck(jac->bs <= 0 || bs == jac->bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Blocksize of x vector %" PetscInt_FMT " does not match fieldsplit blocksize %" PetscInt_FMT, bs, jac->bs);
      PetscCall(VecGetBlockSize(y, &bs));
      PetscCheck(jac->bs <= 0 || bs == jac->bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Blocksize of y vector %" PetscInt_FMT " does not match fieldsplit blocksize %" PetscInt_FMT, bs, jac->bs);
      PetscCall(VecStrideGatherAll(x, jac->x, INSERT_VALUES));
      while (ilink) {
        PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
        PetscCall(KSPSolveTranspose(ilink->ksp, ilink->x, ilink->y));
        PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
        PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
        ilink = ilink->next;
      }
      PetscCall(VecStrideScatterAll(jac->y, y, INSERT_VALUES));
    } else {
      PetscCall(VecSet(y, 0.0));
      while (ilink) {
        PetscCall(FieldSplitSplitSolveAddTranspose(ilink, x, y));
        ilink = ilink->next;
      }
    }
  } else {
    if (!jac->w1) {
      PetscCall(VecDuplicate(x, &jac->w1));
      PetscCall(VecDuplicate(x, &jac->w2));
    }
    PetscCall(VecSet(y, 0.0));
    if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
      PetscCall(FieldSplitSplitSolveAddTranspose(ilink, x, y));
      while (ilink->next) {
        ilink = ilink->next;
        PetscCall(MatMultTranspose(pc->mat, y, jac->w1));
        PetscCall(VecWAXPY(jac->w2, -1.0, jac->w1, x));
        PetscCall(FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y));
      }
      while (ilink->previous) {
        ilink = ilink->previous;
        PetscCall(MatMultTranspose(pc->mat, y, jac->w1));
        PetscCall(VecWAXPY(jac->w2, -1.0, jac->w1, x));
        PetscCall(FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y));
      }
    } else {
      while (ilink->next) { /* get to last entry in linked list */
        ilink = ilink->next;
      }
      PetscCall(FieldSplitSplitSolveAddTranspose(ilink, x, y));
      while (ilink->previous) {
        ilink = ilink->previous;
        PetscCall(MatMultTranspose(pc->mat, y, jac->w1));
        PetscCall(VecWAXPY(jac->w2, -1.0, jac->w1, x));
        PetscCall(FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y));
      }
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCReset_FieldSplit(PC pc)
{
  PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
  PC_FieldSplitLink ilink = jac->head, next;

  PetscFunctionBegin;
  while (ilink) {
    PetscCall(KSPDestroy(&ilink->ksp));
    PetscCall(VecDestroy(&ilink->x));
    PetscCall(VecDestroy(&ilink->y));
    PetscCall(VecDestroy(&ilink->z));
    PetscCall(MatDestroy(&ilink->X));
    PetscCall(MatDestroy(&ilink->Y));
    PetscCall(VecScatterDestroy(&ilink->sctx));
    PetscCall(ISDestroy(&ilink->is));
    PetscCall(ISDestroy(&ilink->is_col));
    PetscCall(PetscFree(ilink->splitname));
    PetscCall(PetscFree(ilink->fields));
    PetscCall(PetscFree(ilink->fields_col));
    next = ilink->next;
    PetscCall(PetscFree(ilink));
    ilink = next;
  }
  jac->head = NULL;
  PetscCall(PetscFree2(jac->x, jac->y));
  if (jac->mat && jac->mat != jac->pmat) {
    PetscCall(MatDestroyMatrices(jac->nsplits, &jac->mat));
  } else if (jac->mat) {
    jac->mat = NULL;
  }
  if (jac->pmat) PetscCall(MatDestroyMatrices(jac->nsplits, &jac->pmat));
  if (jac->Afield) PetscCall(MatDestroyMatrices(jac->nsplits, &jac->Afield));
  jac->nsplits = 0;
  PetscCall(VecDestroy(&jac->w1));
  PetscCall(VecDestroy(&jac->w2));
  if (jac->schur) PetscCall(PetscObjectCompose((PetscObject)jac->schur, "AinvB", NULL));
  PetscCall(MatDestroy(&jac->schur));
  PetscCall(MatDestroy(&jac->schurp));
  PetscCall(MatDestroy(&jac->schur_user));
  PetscCall(KSPDestroy(&jac->kspschur));
  PetscCall(KSPDestroy(&jac->kspupper));
  PetscCall(MatDestroy(&jac->B));
  PetscCall(MatDestroy(&jac->C));
  PetscCall(MatDestroy(&jac->H));
  PetscCall(VecDestroy(&jac->u));
  PetscCall(VecDestroy(&jac->v));
  PetscCall(VecDestroy(&jac->Hu));
  PetscCall(VecDestroy(&jac->d));
  PetscCall(PetscFree(jac->vecz));
  PetscCall(PetscViewerDestroy(&jac->gkbviewer));
  jac->isrestrict = PETSC_FALSE;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCDestroy_FieldSplit(PC pc)
{
  PetscFunctionBegin;
  PetscCall(PCReset_FieldSplit(pc));
  PetscCall(PetscFree(pc->data));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetFields_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetType_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetBlockSize_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitRestrictIS_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSchurGetSubKSP_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", NULL));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc, PetscOptionItems PetscOptionsObject)
{
  PetscInt        bs;
  PetscBool       flg;
  PC_FieldSplit  *jac = (PC_FieldSplit *)pc->data;
  PCCompositeType ctype;

  PetscFunctionBegin;
  PetscOptionsHeadBegin(PetscOptionsObject, "FieldSplit options");
  PetscCall(PetscOptionsBool("-pc_fieldsplit_dm_splits", "Whether to use DMCreateFieldDecomposition() for splits", "PCFieldSplitSetDMSplits", jac->dm_splits, &jac->dm_splits, NULL));
  PetscCall(PetscOptionsInt("-pc_fieldsplit_block_size", "Blocksize that defines number of fields", "PCFieldSplitSetBlockSize", jac->bs, &bs, &flg));
  if (flg) PetscCall(PCFieldSplitSetBlockSize(pc, bs));
  jac->diag_use_amat = pc->useAmat;
  PetscCall(PetscOptionsBool("-pc_fieldsplit_diag_use_amat", "Use Amat (not Pmat) to extract diagonal fieldsplit blocks", "PCFieldSplitSetDiagUseAmat", jac->diag_use_amat, &jac->diag_use_amat, NULL));
  jac->offdiag_use_amat = pc->useAmat;
  PetscCall(PetscOptionsBool("-pc_fieldsplit_off_diag_use_amat", "Use Amat (not Pmat) to extract off-diagonal fieldsplit blocks", "PCFieldSplitSetOffDiagUseAmat", jac->offdiag_use_amat, &jac->offdiag_use_amat, NULL));
  PetscCall(PetscOptionsBool("-pc_fieldsplit_detect_saddle_point", "Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint", jac->detect, &jac->detect, NULL));
  PetscCall(PCFieldSplitSetDetectSaddlePoint(pc, jac->detect)); /* Sets split type and Schur PC type */
  PetscCall(PetscOptionsEnum("-pc_fieldsplit_type", "Type of composition", "PCFieldSplitSetType", PCCompositeTypes, (PetscEnum)jac->type, (PetscEnum *)&ctype, &flg));
  if (flg) PetscCall(PCFieldSplitSetType(pc, ctype));
  /* Only setup fields once */
  if (jac->bs > 0 && jac->nsplits == 0) {
    /* only allow user to set fields from command line.
       otherwise user can set them in PCFieldSplitSetDefaults() */
    PetscCall(PCFieldSplitSetRuntimeSplits_Private(pc));
    if (jac->splitdefined) PetscCall(PetscInfo(pc, "Splits defined using the options database\n"));
  }
  if (jac->type == PC_COMPOSITE_SCHUR) {
    PetscCall(PetscOptionsGetEnum(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_schur_factorization_type", PCFieldSplitSchurFactTypes, (PetscEnum *)&jac->schurfactorization, &flg));
    if (flg) PetscCall(PetscInfo(pc, "Deprecated use of -pc_fieldsplit_schur_factorization_type\n"));
    PetscCall(PetscOptionsEnum("-pc_fieldsplit_schur_fact_type", "Which off-diagonal parts of the block factorization to use", "PCFieldSplitSetSchurFactType", PCFieldSplitSchurFactTypes, (PetscEnum)jac->schurfactorization, (PetscEnum *)&jac->schurfactorization, NULL));
    PetscCall(PetscOptionsEnum("-pc_fieldsplit_schur_precondition", "How to build preconditioner for Schur complement", "PCFieldSplitSetSchurPre", PCFieldSplitSchurPreTypes, (PetscEnum)jac->schurpre, (PetscEnum *)&jac->schurpre, NULL));
    PetscCall(PetscOptionsScalar("-pc_fieldsplit_schur_scale", "Scale Schur complement", "PCFieldSplitSetSchurScale", jac->schurscale, &jac->schurscale, NULL));
  } else if (jac->type == PC_COMPOSITE_GKB) {
    PetscCall(PetscOptionsReal("-pc_fieldsplit_gkb_tol", "The tolerance for the lower bound stopping criterion", "PCFieldSplitSetGKBTol", jac->gkbtol, &jac->gkbtol, NULL));
    PetscCall(PetscOptionsInt("-pc_fieldsplit_gkb_delay", "The delay value for lower bound criterion", "PCFieldSplitSetGKBDelay", jac->gkbdelay, &jac->gkbdelay, NULL));
    PetscCall(PetscOptionsBoundedReal("-pc_fieldsplit_gkb_nu", "Parameter in augmented Lagrangian approach", "PCFieldSplitSetGKBNu", jac->gkbnu, &jac->gkbnu, NULL, 0.0));
    PetscCall(PetscOptionsInt("-pc_fieldsplit_gkb_maxit", "Maximum allowed number of iterations", "PCFieldSplitSetGKBMaxit", jac->gkbmaxit, &jac->gkbmaxit, NULL));
    PetscCall(PetscOptionsBool("-pc_fieldsplit_gkb_monitor", "Prints number of GKB iterations and error", "PCFieldSplitGKB", jac->gkbmonitor, &jac->gkbmonitor, NULL));
  }
  /*
    In the initial call to this routine the sub-solver data structures do not exist so we cannot call KSPSetFromOptions() on them yet.
    But after the initial setup of ALL the layers of sub-solvers is completed we do want to call KSPSetFromOptions() on the sub-solvers every time it
    is called on the outer solver in case changes were made in the options database

    But even after PCSetUp_FieldSplit() is called all the options inside the inner levels of sub-solvers may still not have been set thus we only call the KSPSetFromOptions()
    if we know that the entire stack of sub-solvers below this have been complete instantiated, we check this by seeing if any solver iterations are complete.
    Without this extra check test p2p1fetidp_olof_full and others fail with incorrect matrix types.

    There could be a negative side effect of calling the KSPSetFromOptions() below.

    If one captured the PetscObjectState of the options database one could skip these calls if the database has not changed from the previous call
  */
  if (jac->issetup) {
    PC_FieldSplitLink ilink = jac->head;
    if (jac->type == PC_COMPOSITE_SCHUR) {
      if (jac->kspupper && jac->kspupper->totalits > 0) PetscCall(KSPSetFromOptions(jac->kspupper));
      if (jac->kspschur && jac->kspschur->totalits > 0) PetscCall(KSPSetFromOptions(jac->kspschur));
    }
    while (ilink) {
      if (ilink->ksp->totalits > 0) PetscCall(KSPSetFromOptions(ilink->ksp));
      ilink = ilink->next;
    }
  }
  PetscOptionsHeadEnd();
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc, const char splitname[], PetscInt n, const PetscInt *fields, const PetscInt *fields_col)
{
  PC_FieldSplit    *jac = (PC_FieldSplit *)pc->data;
  PC_FieldSplitLink ilink, next = jac->head;
  char              prefix[128];
  PetscInt          i;
  PetscLogEvent     nse;

  PetscFunctionBegin;
  if (jac->splitdefined) {
    PetscCall(PetscInfo(pc, "Ignoring new split \"%s\" because the splits have already been defined\n", splitname));
    PetscFunctionReturn(PETSC_SUCCESS);
  }
  for (i = 0; i < n; i++) PetscCheck(fields[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative field %" PetscInt_FMT " requested", fields[i]);
  PetscCall(PetscNew(&ilink));
  if (splitname) {
    PetscCall(PetscStrallocpy(splitname, &ilink->splitname));
  } else {
    PetscCall(PetscMalloc1(3, &ilink->splitname));
    PetscCall(PetscSNPrintf(ilink->splitname, 2, "%" PetscInt_FMT, jac->nsplits));
  }
  PetscCall(PetscMPIIntCast(jac->nsplits, &nse));
  ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + nse : KSP_Solve_FS_0 + 4; /* Splits greater than 4 logged in 4th split */
  PetscCall(PetscMalloc1(n, &ilink->fields));
  PetscCall(PetscArraycpy(ilink->fields, fields, n));
  PetscCall(PetscMalloc1(n, &ilink->fields_col));
  PetscCall(PetscArraycpy(ilink->fields_col, fields_col, n));

  ilink->nfields = n;
  ilink->next    = NULL;
  PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &ilink->ksp));
  PetscCall(KSPSetNestLevel(ilink->ksp, pc->kspnestlevel));
  PetscCall(KSPSetErrorIfNotConverged(ilink->ksp, pc->erroriffailure));
  PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)pc, 1));
  PetscCall(KSPSetType(ilink->ksp, KSPPREONLY));

  PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
  PetscCall(KSPSetOptionsPrefix(ilink->ksp, prefix));

  if (!next) {
    jac->head       = ilink;
    ilink->previous = NULL;
  } else {
    while (next->next) next = next->next;
    next->next      = ilink;
    ilink->previous = next;
  }
  jac->nsplits++;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc, PetscInt *n, KSP **subksp)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

  PetscFunctionBegin;
  *subksp = NULL;
  if (n) *n = 0;
  if (jac->type == PC_COMPOSITE_SCHUR) {
    PetscInt nn;

    PetscCheck(jac->schur, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()");
    PetscCheck(jac->nsplits == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Unexpected number of splits %" PetscInt_FMT " != 2", jac->nsplits);
    nn = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0);
    PetscCall(PetscMalloc1(nn, subksp));
    (*subksp)[0] = jac->head->ksp;
    (*subksp)[1] = jac->kspschur;
    if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper;
    if (n) *n = nn;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc, PetscInt *n, KSP **subksp)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

  PetscFunctionBegin;
  PetscCheck(jac->schur, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitGetSubKSP()");
  PetscCall(PetscMalloc1(jac->nsplits, subksp));
  PetscCall(MatSchurComplementGetKSP(jac->schur, *subksp));

  (*subksp)[1] = jac->kspschur;
  if (n) *n = jac->nsplits;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc, PetscInt *n, KSP **subksp)
{
  PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
  PetscInt          cnt   = 0;
  PC_FieldSplitLink ilink = jac->head;

  PetscFunctionBegin;
  PetscCall(PetscMalloc1(jac->nsplits, subksp));
  while (ilink) {
    (*subksp)[cnt++] = ilink->ksp;
    ilink            = ilink->next;
  }
  PetscCheck(cnt == jac->nsplits, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupt PCFIELDSPLIT object: number of splits in linked list %" PetscInt_FMT " does not match number in object %" PetscInt_FMT, cnt, jac->nsplits);
  if (n) *n = jac->nsplits;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitRestrictIS - Restricts the fieldsplit `IS`s to be within a given `IS`.

  Input Parameters:
+ pc  - the preconditioner context
- isy - the index set that defines the indices to which the fieldsplit is to be restricted

  Level: advanced

  Developer Notes:
  It seems the resulting `IS`s will not cover the entire space, so
  how can they define a convergent preconditioner? Needs explaining.

.seealso: [](sec_block_matrices), `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
@*/
PetscErrorCode PCFieldSplitRestrictIS(PC pc, IS isy)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscValidHeaderSpecific(isy, IS_CLASSID, 2);
  PetscTryMethod(pc, "PCFieldSplitRestrictIS_C", (PC, IS), (pc, isy));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
{
  PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
  PC_FieldSplitLink ilink = jac->head, next;
  PetscInt          localsize, size, sizez, i;
  const PetscInt   *ind, *indz;
  PetscInt         *indc, *indcz;
  PetscBool         flg;

  PetscFunctionBegin;
  PetscCall(ISGetLocalSize(isy, &localsize));
  PetscCallMPI(MPI_Scan(&localsize, &size, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)isy)));
  size -= localsize;
  while (ilink) {
    IS isrl, isr;
    PC subpc;
    PetscCall(ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl));
    PetscCall(ISGetLocalSize(isrl, &localsize));
    PetscCall(PetscMalloc1(localsize, &indc));
    PetscCall(ISGetIndices(isrl, &ind));
    PetscCall(PetscArraycpy(indc, ind, localsize));
    PetscCall(ISRestoreIndices(isrl, &ind));
    PetscCall(ISDestroy(&isrl));
    for (i = 0; i < localsize; i++) *(indc + i) += size;
    PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)isy), localsize, indc, PETSC_OWN_POINTER, &isr));
    PetscCall(PetscObjectReference((PetscObject)isr));
    PetscCall(ISDestroy(&ilink->is));
    ilink->is = isr;
    PetscCall(PetscObjectReference((PetscObject)isr));
    PetscCall(ISDestroy(&ilink->is_col));
    ilink->is_col = isr;
    PetscCall(ISDestroy(&isr));
    PetscCall(KSPGetPC(ilink->ksp, &subpc));
    PetscCall(PetscObjectTypeCompare((PetscObject)subpc, PCFIELDSPLIT, &flg));
    if (flg) {
      IS       iszl, isz;
      MPI_Comm comm;
      PetscCall(ISGetLocalSize(ilink->is, &localsize));
      comm = PetscObjectComm((PetscObject)ilink->is);
      PetscCall(ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl));
      PetscCallMPI(MPI_Scan(&localsize, &sizez, 1, MPIU_INT, MPI_SUM, comm));
      sizez -= localsize;
      PetscCall(ISGetLocalSize(iszl, &localsize));
      PetscCall(PetscMalloc1(localsize, &indcz));
      PetscCall(ISGetIndices(iszl, &indz));
      PetscCall(PetscArraycpy(indcz, indz, localsize));
      PetscCall(ISRestoreIndices(iszl, &indz));
      PetscCall(ISDestroy(&iszl));
      for (i = 0; i < localsize; i++) *(indcz + i) += sizez;
      PetscCall(ISCreateGeneral(comm, localsize, indcz, PETSC_OWN_POINTER, &isz));
      PetscCall(PCFieldSplitRestrictIS(subpc, isz));
      PetscCall(ISDestroy(&isz));
    }
    next  = ilink->next;
    ilink = next;
  }
  jac->isrestrict = PETSC_TRUE;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc, const char splitname[], IS is)
{
  PC_FieldSplit    *jac = (PC_FieldSplit *)pc->data;
  PC_FieldSplitLink ilink, next = jac->head;
  char              prefix[128];
  PetscLogEvent     nse;

  PetscFunctionBegin;
  if (jac->splitdefined) {
    PetscCall(PetscInfo(pc, "Ignoring new split \"%s\" because the splits have already been defined\n", splitname));
    PetscFunctionReturn(PETSC_SUCCESS);
  }
  PetscCall(PetscNew(&ilink));
  if (splitname) {
    PetscCall(PetscStrallocpy(splitname, &ilink->splitname));
  } else {
    PetscCall(PetscMalloc1(8, &ilink->splitname));
    PetscCall(PetscSNPrintf(ilink->splitname, 7, "%" PetscInt_FMT, jac->nsplits));
  }
  PetscCall(PetscMPIIntCast(jac->nsplits, &nse));
  ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + nse : KSP_Solve_FS_0 + 4; /* Splits greater than 4 logged in 4th split */
  PetscCall(PetscObjectReference((PetscObject)is));
  PetscCall(ISDestroy(&ilink->is));
  ilink->is = is;
  PetscCall(PetscObjectReference((PetscObject)is));
  PetscCall(ISDestroy(&ilink->is_col));
  ilink->is_col = is;
  ilink->next   = NULL;
  PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &ilink->ksp));
  PetscCall(KSPSetNestLevel(ilink->ksp, pc->kspnestlevel));
  PetscCall(KSPSetErrorIfNotConverged(ilink->ksp, pc->erroriffailure));
  PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)pc, 1));
  PetscCall(KSPSetType(ilink->ksp, KSPPREONLY));

  PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
  PetscCall(KSPSetOptionsPrefix(ilink->ksp, prefix));

  if (!next) {
    jac->head       = ilink;
    ilink->previous = NULL;
  } else {
    while (next->next) next = next->next;
    next->next      = ilink;
    ilink->previous = next;
  }
  jac->nsplits++;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitSetFields - Sets the fields that define one particular split in `PCFIELDSPLIT`

  Logically Collective

  Input Parameters:
+ pc         - the preconditioner context
. splitname  - name of this split, if `NULL` the number of the split is used
. n          - the number of fields in this split
. fields     - the fields in this split
- fields_col - generally the same as `fields`, if it does not match `fields` then the submatrix that is solved for this set of fields comes from an off-diagonal block
               of the matrix and `fields_col` provides the column indices for that block

  Options Database Key:
. -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the `%d`'th split

  Level: intermediate

  Notes:
  Use `PCFieldSplitSetIS()` to set a  general set of indices as a split.

  If the matrix used to construct the preconditioner is `MATNEST` then field i refers to the `is_row[i]` `IS` passed to `MatCreateNest()`.

  If the matrix used to construct the preconditioner is not `MATNEST` then
  `PCFieldSplitSetFields()` is for defining fields as strided blocks (based on the block size provided to the matrix with `MatSetBlockSize()` or
  to the `PC` with `PCFieldSplitSetBlockSize()`). For example, if the block
  size is three then one can define a split as 0, or 1 or 2 or 0,1 or 0,2 or 1,2 which mean
  0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x23x56x8.. x12x45x78x....
  where the numbered entries indicate what is in the split.

  This function is called once per split (it creates a new split each time).  Solve options
  for this split will be available under the prefix `-fieldsplit_SPLITNAME_`.

  `PCFieldSplitSetIS()` does not support having a `fields_col` different from `fields`

  Developer Notes:
  This routine does not actually create the `IS` representing the split, that is delayed
  until `PCSetUp_FieldSplit()`, because information about the vector/matrix layouts may not be
  available when this routine is called.

.seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`, `PCFieldSplitSetIS()`, `PCFieldSplitRestrictIS()`,
          `MatSetBlockSize()`, `MatCreateNest()`
@*/
PetscErrorCode PCFieldSplitSetFields(PC pc, const char splitname[], PetscInt n, const PetscInt fields[], const PetscInt fields_col[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscAssertPointer(splitname, 2);
  PetscCheck(n >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Provided number of fields %" PetscInt_FMT " in split \"%s\" not positive", n, splitname);
  PetscAssertPointer(fields, 4);
  PetscTryMethod(pc, "PCFieldSplitSetFields_C", (PC, const char[], PetscInt, const PetscInt *, const PetscInt *), (pc, splitname, n, fields, fields_col));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) to build
  the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators.

  Logically Collective

  Input Parameters:
+ pc  - the preconditioner object
- flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from

  Options Database Key:
. -pc_fieldsplit_diag_use_amat - use the Amat to provide the diagonal blocks

  Level: intermediate

.seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitGetDiagUseAmat()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFIELDSPLIT`
@*/
PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc, PetscBool flg)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
  PetscBool      isfs;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
  PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT);
  jac->diag_use_amat = flg;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) to build
  the sub-matrices associated with each split.  Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators.

  Logically Collective

  Input Parameter:
. pc - the preconditioner object

  Output Parameter:
. flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from

  Level: intermediate

.seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitSetDiagUseAmat()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFIELDSPLIT`
@*/
PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc, PetscBool *flg)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
  PetscBool      isfs;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscAssertPointer(flg, 2);
  PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
  PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT);
  *flg = jac->diag_use_amat;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) to build
  the sub-matrices associated with each split.  Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators.

  Logically Collective

  Input Parameters:
+ pc  - the preconditioner object
- flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from

  Options Database Key:
. -pc_fieldsplit_off_diag_use_amat <bool> - use the Amat to extract the off-diagonal blocks

  Level: intermediate

.seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFieldSplitSetDiagUseAmat()`, `PCFIELDSPLIT`
@*/
PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc, PetscBool flg)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
  PetscBool      isfs;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
  PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT);
  jac->offdiag_use_amat = flg;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) to build
  the sub-matrices associated with each split.  Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators.

  Logically Collective

  Input Parameter:
. pc - the preconditioner object

  Output Parameter:
. flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from

  Level: intermediate

.seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFieldSplitGetDiagUseAmat()`, `PCFIELDSPLIT`
@*/
PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc, PetscBool *flg)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
  PetscBool      isfs;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscAssertPointer(flg, 2);
  PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
  PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT);
  *flg = jac->offdiag_use_amat;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitSetIS - Sets the exact elements for a split in a `PCFIELDSPLIT`

  Logically Collective

  Input Parameters:
+ pc        - the preconditioner context
. splitname - name of this split, if `NULL` the number of the split is used
- is        - the index set that defines the elements in this split

  Level: intermediate

  Notes:
  Use `PCFieldSplitSetFields()`, for splits defined by strided `IS` based on the matrix block size or the `is_rows[]` passed into `MATNEST`

  This function is called once per split (it creates a new split each time).  Solve options
  for this split will be available under the prefix -fieldsplit_SPLITNAME_.

.seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`, `PCFieldSplitSetFields()`
@*/
PetscErrorCode PCFieldSplitSetIS(PC pc, const char splitname[], IS is)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  if (splitname) PetscAssertPointer(splitname, 2);
  PetscValidHeaderSpecific(is, IS_CLASSID, 3);
  PetscTryMethod(pc, "PCFieldSplitSetIS_C", (PC, const char[], IS), (pc, splitname, is));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitGetIS - Retrieves the elements for a split as an `IS`

  Logically Collective

  Input Parameters:
+ pc        - the preconditioner context
- splitname - name of this split

  Output Parameter:
. is - the index set that defines the elements in this split, or `NULL` if the split is not found

  Level: intermediate

.seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetIS()`, `PCFieldSplitGetISByIndex()`
@*/
PetscErrorCode PCFieldSplitGetIS(PC pc, const char splitname[], IS *is)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscAssertPointer(splitname, 2);
  PetscAssertPointer(is, 3);
  {
    PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
    PC_FieldSplitLink ilink = jac->head;
    PetscBool         found;

    *is = NULL;
    while (ilink) {
      PetscCall(PetscStrcmp(ilink->splitname, splitname, &found));
      if (found) {
        *is = ilink->is;
        break;
      }
      ilink = ilink->next;
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitGetISByIndex - Retrieves the elements for a given split as an `IS`

  Logically Collective

  Input Parameters:
+ pc    - the preconditioner context
- index - index of this split

  Output Parameter:
. is - the index set that defines the elements in this split

  Level: intermediate

.seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitGetIS()`, `PCFieldSplitSetIS()`,

@*/
PetscErrorCode PCFieldSplitGetISByIndex(PC pc, PetscInt index, IS *is)
{
  PetscFunctionBegin;
  PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative field %" PetscInt_FMT " requested", index);
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscAssertPointer(is, 3);
  {
    PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
    PC_FieldSplitLink ilink = jac->head;
    PetscInt          i     = 0;
    PetscCheck(index < jac->nsplits, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT " requested but only %" PetscInt_FMT " exist", index, jac->nsplits);

    while (i < index) {
      ilink = ilink->next;
      ++i;
    }
    PetscCall(PCFieldSplitGetIS(pc, ilink->splitname, is));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the
  fieldsplit preconditioner when calling `PCFieldSplitSetFields()`. If not set the matrix block size is used.

  Logically Collective

  Input Parameters:
+ pc - the preconditioner context
- bs - the block size

  Level: intermediate

  Note:
  If the matrix is a `MATNEST` then the `is_rows[]` passed to `MatCreateNest()` determines the fields.

.seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
@*/
PetscErrorCode PCFieldSplitSetBlockSize(PC pc, PetscInt bs)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscValidLogicalCollectiveInt(pc, bs, 2);
  PetscTryMethod(pc, "PCFieldSplitSetBlockSize_C", (PC, PetscInt), (pc, bs));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PCFieldSplitGetSubKSP - Gets the `KSP` contexts for all splits

  Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameters:
+ n      - the number of splits
- subksp - the array of `KSP` contexts

  Level: advanced

  Notes:
  After `PCFieldSplitGetSubKSP()` the array of `KSP`s is to be freed by the user with `PetscFree()`
  (not the `KSP`, just the array that contains them).

  You must call `PCSetUp()` before calling `PCFieldSplitGetSubKSP()`.

  If the fieldsplit is of type `PC_COMPOSITE_SCHUR`, it returns the `KSP` object used inside the
  Schur complement and the `KSP` object used to iterate over the Schur complement.
  To access all the `KSP` objects used in `PC_COMPOSITE_SCHUR`, use `PCFieldSplitSchurGetSubKSP()`.

  If the fieldsplit is of type `PC_COMPOSITE_GKB`, it returns the `KSP` object used to solve the
  inner linear system defined by the matrix H in each loop.

  Fortran Note:
  Call `PCFieldSplitRestoreSubKSP()` when the array of `KSP` is no longer needed

  Developer Notes:
  There should be a `PCFieldSplitRestoreSubKSP()` instead of requiring the user to call `PetscFree()`

.seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`, `PCFieldSplitSchurGetSubKSP()`
@*/
PetscErrorCode PCFieldSplitGetSubKSP(PC pc, PetscInt *n, KSP *subksp[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  if (n) PetscAssertPointer(n, 2);
  PetscUseMethod(pc, "PCFieldSplitGetSubKSP_C", (PC, PetscInt *, KSP **), (pc, n, subksp));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PCFieldSplitSchurGetSubKSP - Gets the `KSP` contexts used inside the Schur complement based `PCFIELDSPLIT`

  Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameters:
+ n      - the number of splits
- subksp - the array of `KSP` contexts

  Level: advanced

  Notes:
  After `PCFieldSplitSchurGetSubKSP()` the array of `KSP`s is to be freed by the user with `PetscFree()`
  (not the `KSP` just the array that contains them).

  You must call `PCSetUp()` before calling `PCFieldSplitSchurGetSubKSP()`.

  If the fieldsplit type is of type `PC_COMPOSITE_SCHUR`, it returns (in order)
+  1  - the `KSP` used for the (1,1) block
.  2  - the `KSP` used for the Schur complement (not the one used for the interior Schur solver)
-  3  - the `KSP` used for the (1,1) block in the upper triangular factor (if different from that of the (1,1) block).

  It returns a null array if the fieldsplit is not of type `PC_COMPOSITE_SCHUR`; in this case, you should use `PCFieldSplitGetSubKSP()`.

  Fortran Note:
  Call `PCFieldSplitSchurRestoreSubKSP()` when the array of `KSP` is no longer needed

  Developer Notes:
  There should be a `PCFieldSplitRestoreSubKSP()` instead of requiring the user to call `PetscFree()`

  Should the functionality of `PCFieldSplitSchurGetSubKSP()` and `PCFieldSplitGetSubKSP()` be merged?

.seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`, `PCFieldSplitGetSubKSP()`
@*/
PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc, PetscInt *n, KSP *subksp[])
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  if (n) PetscAssertPointer(n, 2);
  PetscUseMethod(pc, "PCFieldSplitSchurGetSubKSP_C", (PC, PetscInt *, KSP **), (pc, n, subksp));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitSetSchurPre -  Indicates from what operator the preconditioner is constructed for the Schur complement.
  The default is the A11 matrix.

  Collective

  Input Parameters:
+ pc    - the preconditioner context
. ptype - which matrix to use for preconditioning the Schur complement: `PC_FIELDSPLIT_SCHUR_PRE_A11` (default),
              `PC_FIELDSPLIT_SCHUR_PRE_SELF`, `PC_FIELDSPLIT_SCHUR_PRE_USER`,
              `PC_FIELDSPLIT_SCHUR_PRE_SELFP`, and `PC_FIELDSPLIT_SCHUR_PRE_FULL`
- pre   - matrix to use for preconditioning, or `NULL`

  Options Database Keys:
+ -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is `a11`. See notes for meaning of various arguments
- -fieldsplit_1_pc_type <pctype>                               - the preconditioner algorithm that is used to construct the preconditioner from the operator

  Level: intermediate

  Notes:
  If ptype is
+     a11 - the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner
  matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix
.     self - the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix:
  The only preconditioners that currently work with this symbolic representation matrix object are `PCLSC` and `PCHPDDM`
.     user - the preconditioner for the Schur complement is generated from the user provided matrix (pre argument
  to this function).
.     selfp - the preconditioning for the Schur complement is generated from an explicitly-assembled approximation $ Sp = A11 - A10 inv(diag(A00)) A01 $
  This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
  lumped before extracting the diagonal using the additional option `-fieldsplit_1_mat_schur_complement_ainv_type lump`
-     full - the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation
  computed internally by `PCFIELDSPLIT` (this is expensive)
  useful mostly as a test that the Schur complement approach can work for your problem

  When solving a saddle point problem, where the A11 block is identically zero, using `a11` as the ptype only makes sense
  with the additional option `-fieldsplit_1_pc_type none`. Usually for saddle point problems one would use a `ptype` of `self` and
  `-fieldsplit_1_pc_type lsc` which uses the least squares commutator to compute a preconditioner for the Schur complement.

  Developer Note:
  The name of this function and the option `-pc_fieldsplit_schur_precondition` are inconsistent; precondition should be used everywhere.

.seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`,
          `MatSchurComplementSetAinvType()`, `PCLSC`, `PCFieldSplitSetSchurFactType()`
@*/
PetscErrorCode PCFieldSplitSetSchurPre(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscTryMethod(pc, "PCFieldSplitSetSchurPre_C", (PC, PCFieldSplitSchurPreType, Mat), (pc, ptype, pre));
  PetscFunctionReturn(PETSC_SUCCESS);
}

PetscErrorCode PCFieldSplitSchurPrecondition(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
{
  return PCFieldSplitSetSchurPre(pc, ptype, pre);
} /* Deprecated name */

/*@
  PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be
  preconditioned.  See `PCFieldSplitSetSchurPre()` for details.

  Logically Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameters:
+ ptype - which matrix to use for preconditioning the Schur complement: `PC_FIELDSPLIT_SCHUR_PRE_A11`, `PC_FIELDSPLIT_SCHUR_PRE_SELF`, `PC_FIELDSPLIT_SCHUR_PRE_USER`
- pre   - matrix to use for preconditioning (with `PC_FIELDSPLIT_SCHUR_PRE_USER`), or `NULL`

  Level: intermediate

.seealso: [](sec_block_matrices), `PC`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCLSC`
@*/
PetscErrorCode PCFieldSplitGetSchurPre(PC pc, PCFieldSplitSchurPreType *ptype, Mat *pre)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscUseMethod(pc, "PCFieldSplitGetSchurPre_C", (PC, PCFieldSplitSchurPreType *, Mat *), (pc, ptype, pre));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitSchurGetS -  extract the `MATSCHURCOMPLEMENT` object used by this `PCFIELDSPLIT` in case it needs to be configured separately

  Not Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameter:
. S - the Schur complement matrix

  Level: advanced

  Note:
  This matrix should not be destroyed using `MatDestroy()`; rather, use `PCFieldSplitSchurRestoreS()`.

.seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MATSCHURCOMPLEMENT`, `PCFieldSplitSchurRestoreS()`,
          `MatCreateSchurComplement()`, `MatSchurComplementGetKSP()`, `MatSchurComplementComputeExplicitOperator()`, `MatGetSchurComplement()`
@*/
PetscErrorCode PCFieldSplitSchurGetS(PC pc, Mat *S)
{
  const char    *t;
  PetscBool      isfs;
  PC_FieldSplit *jac;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscCall(PetscObjectGetType((PetscObject)pc, &t));
  PetscCall(PetscStrcmp(t, PCFIELDSPLIT, &isfs));
  PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PC of type PCFIELDSPLIT, got %s instead", t);
  jac = (PC_FieldSplit *)pc->data;
  PetscCheck(jac->type == PC_COMPOSITE_SCHUR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PCFIELDSPLIT of type SCHUR, got %d instead", jac->type);
  if (S) *S = jac->schur;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitSchurRestoreS -  returns the `MATSCHURCOMPLEMENT` matrix used by this `PC`

  Not Collective

  Input Parameters:
+ pc - the preconditioner context
- S  - the Schur complement matrix

  Level: advanced

.seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MatSchurComplement`, `PCFieldSplitSchurGetS()`
@*/
PetscErrorCode PCFieldSplitSchurRestoreS(PC pc, Mat *S)
{
  const char    *t;
  PetscBool      isfs;
  PC_FieldSplit *jac;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscCall(PetscObjectGetType((PetscObject)pc, &t));
  PetscCall(PetscStrcmp(t, PCFIELDSPLIT, &isfs));
  PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PC of type PCFIELDSPLIT, got %s instead", t);
  jac = (PC_FieldSplit *)pc->data;
  PetscCheck(jac->type == PC_COMPOSITE_SCHUR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PCFIELDSPLIT of type SCHUR, got %d instead", jac->type);
  PetscCheck(S && (*S == jac->schur), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatSchurComplement restored is not the same as gotten");
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

  PetscFunctionBegin;
  jac->schurpre = ptype;
  if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
    PetscCall(MatDestroy(&jac->schur_user));
    jac->schur_user = pre;
    PetscCall(PetscObjectReference((PetscObject)jac->schur_user));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc, PCFieldSplitSchurPreType *ptype, Mat *pre)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

  PetscFunctionBegin;
  if (ptype) *ptype = jac->schurpre;
  if (pre) *pre = jac->schur_user;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitSetSchurFactType -  sets which blocks of the approximate block factorization to retain in the preconditioner {cite}`murphy2000note` and {cite}`ipsen2001note`

  Collective

  Input Parameters:
+ pc    - the preconditioner context
- ftype - which blocks of factorization to retain, `PC_FIELDSPLIT_SCHUR_FACT_FULL` is default

  Options Database Key:
. -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - default is `full`

  Level: intermediate

  Notes:
  The `full` factorization is

  ```{math}
  \left(\begin{array}{cc} A & B \\
  C & E \\
  \end{array}\right) =
  \left(\begin{array}{cc} I & 0 \\
  C A^{-1} & I \\
  \end{array}\right)
  \left(\begin{array}{cc} A & 0 \\
  0 & S \\
  \end{array}\right)
  \left(\begin{array}{cc} I & A^{-1}B \\
  0 & I \\
  \end{array}\right) = L D U,
  ```

  where $ S = E - C A^{-1} B $. In practice, the full factorization is applied via block triangular solves with the grouping $L(DU)$. `upper` uses $DU$, `lower` uses $LD$,
  and `diag` is the diagonal part with the sign of $S$ flipped (because this makes the preconditioner positive definite for many formulations,
  thus allowing the use of `KSPMINRES)`. Sign flipping of $S$ can be turned off with `PCFieldSplitSetSchurScale()`.

  If $A$ and $S$ are solved exactly
+  1 - `full` factorization is a direct solver.
.  2 - The preconditioned operator with `lower` or `upper` has all eigenvalues equal to 1 and minimal polynomial of degree 2, so `KSPGMRES` converges in 2 iterations.
-  3 - With `diag`, the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so `KSPGMRES` converges in at most 4 iterations.

  If the iteration count is very low, consider using `KSPFGMRES` or `KSPGCR` which can use one less preconditioner
  application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice.

  For symmetric problems in which $A$ is positive definite and $S$ is negative definite, `diag` can be used with `KSPMINRES`.

  A flexible method like `KSPFGMRES` or `KSPGCR`, [](sec_flexibleksp), must be used if the fieldsplit preconditioner is nonlinear (e.g., a few iterations of a Krylov method is used to solve with $A$ or $S$).

.seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurScale()`,
          [](sec_flexibleksp), `PCFieldSplitSetSchurPre()`
@*/
PetscErrorCode PCFieldSplitSetSchurFactType(PC pc, PCFieldSplitSchurFactType ftype)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscTryMethod(pc, "PCFieldSplitSetSchurFactType_C", (PC, PCFieldSplitSchurFactType), (pc, ftype));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc, PCFieldSplitSchurFactType ftype)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

  PetscFunctionBegin;
  jac->schurfactorization = ftype;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitSetSchurScale -  Controls the sign flip of S for `PC_FIELDSPLIT_SCHUR_FACT_DIAG`.

  Collective

  Input Parameters:
+ pc    - the preconditioner context
- scale - scaling factor for the Schur complement

  Options Database Key:
. -pc_fieldsplit_schur_scale <scale> - default is -1.0

  Level: intermediate

.seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurFactType`, `PCFieldSplitSetSchurFactType()`
@*/
PetscErrorCode PCFieldSplitSetSchurScale(PC pc, PetscScalar scale)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscValidLogicalCollectiveScalar(pc, scale, 2);
  PetscTryMethod(pc, "PCFieldSplitSetSchurScale_C", (PC, PetscScalar), (pc, scale));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc, PetscScalar scale)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

  PetscFunctionBegin;
  jac->schurscale = scale;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement

  Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameters:
+ A00 - the (0,0) block
. A01 - the (0,1) block
. A10 - the (1,0) block
- A11 - the (1,1) block

  Level: advanced

  Note:
  Use `NULL` for any unneeded output arguments

.seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `MatSchurComplementGetSubMatrices()`, `MatSchurComplementSetSubMatrices()`
@*/
PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc, Mat *A00, Mat *A01, Mat *A10, Mat *A11)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscCheck(jac->type == PC_COMPOSITE_SCHUR, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
  if (A00) *A00 = jac->pmat[0];
  if (A01) *A01 = jac->B;
  if (A10) *A10 = jac->C;
  if (A11) *A11 = jac->pmat[1];
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitSetGKBTol -  Sets the solver tolerance for the generalized Golub-Kahan bidiagonalization preconditioner {cite}`arioli2013` in `PCFIELDSPLIT`

  Collective

  Input Parameters:
+ pc        - the preconditioner context
- tolerance - the solver tolerance

  Options Database Key:
. -pc_fieldsplit_gkb_tol <tolerance> - default is 1e-5

  Level: intermediate

  Note:
  The generalized GKB algorithm {cite}`arioli2013` uses a lower bound estimate of the error in energy norm as stopping criterion.
  It stops once the lower bound estimate undershoots the required solver tolerance. Although the actual error might be bigger than
  this estimate, the stopping criterion is satisfactory in practical cases.

.seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBMaxit()`
@*/
PetscErrorCode PCFieldSplitSetGKBTol(PC pc, PetscReal tolerance)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscValidLogicalCollectiveReal(pc, tolerance, 2);
  PetscTryMethod(pc, "PCFieldSplitSetGKBTol_C", (PC, PetscReal), (pc, tolerance));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCFieldSplitSetGKBTol_FieldSplit(PC pc, PetscReal tolerance)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

  PetscFunctionBegin;
  jac->gkbtol = tolerance;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitSetGKBMaxit -  Sets the maximum number of iterations for the generalized Golub-Kahan bidiagonalization preconditioner {cite}`arioli2013` in `PCFIELDSPLIT`

  Collective

  Input Parameters:
+ pc    - the preconditioner context
- maxit - the maximum number of iterations

  Options Database Key:
. -pc_fieldsplit_gkb_maxit <maxit> - default is 100

  Level: intermediate

.seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBNu()`
@*/
PetscErrorCode PCFieldSplitSetGKBMaxit(PC pc, PetscInt maxit)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscValidLogicalCollectiveInt(pc, maxit, 2);
  PetscTryMethod(pc, "PCFieldSplitSetGKBMaxit_C", (PC, PetscInt), (pc, maxit));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCFieldSplitSetGKBMaxit_FieldSplit(PC pc, PetscInt maxit)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

  PetscFunctionBegin;
  jac->gkbmaxit = maxit;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitSetGKBDelay -  Sets the delay in the lower bound error estimate in the generalized Golub-Kahan bidiagonalization {cite}`arioli2013` in `PCFIELDSPLIT`
  preconditioner.

  Collective

  Input Parameters:
+ pc    - the preconditioner context
- delay - the delay window in the lower bound estimate

  Options Database Key:
. -pc_fieldsplit_gkb_delay <delay> - default is 5

  Level: intermediate

  Notes:
  The algorithm uses a lower bound estimate of the error in energy norm as stopping criterion. The lower bound of the error $ ||u-u^k||_H $
  is expressed as a truncated sum. The error at iteration k can only be measured at iteration (k + `delay`), and thus the algorithm needs
  at least (`delay` + 1) iterations to stop.

  For more details on the generalized Golub-Kahan bidiagonalization method and its lower bound stopping criterion, please refer to {cite}`arioli2013`

.seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()`
@*/
PetscErrorCode PCFieldSplitSetGKBDelay(PC pc, PetscInt delay)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscValidLogicalCollectiveInt(pc, delay, 2);
  PetscTryMethod(pc, "PCFieldSplitSetGKBDelay_C", (PC, PetscInt), (pc, delay));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCFieldSplitSetGKBDelay_FieldSplit(PC pc, PetscInt delay)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

  PetscFunctionBegin;
  jac->gkbdelay = delay;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitSetGKBNu -  Sets the scalar value nu >= 0 in the transformation H = A00 + nu*A01*A01' of the (1,1) block in the
  Golub-Kahan bidiagonalization preconditioner {cite}`arioli2013` in `PCFIELDSPLIT`

  Collective

  Input Parameters:
+ pc - the preconditioner context
- nu - the shift parameter

  Options Database Key:
. -pc_fieldsplit_gkb_nu <nu> - default is 1

  Level: intermediate

  Notes:
  This shift is in general done to obtain better convergence properties for the outer loop of the algorithm. This is often achieved by choosing `nu` sufficiently large. However,
  if `nu` is chosen too large, the matrix H might be badly conditioned and the solution of the linear system $Hx = b$ in the inner loop becomes difficult. It is therefore
  necessary to find a good balance in between the convergence of the inner and outer loop.

  For `nu` = 0, no shift is done. In this case A00 has to be positive definite. The matrix N in {cite}`arioli2013` is then chosen as identity.

.seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()`
@*/
PetscErrorCode PCFieldSplitSetGKBNu(PC pc, PetscReal nu)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscValidLogicalCollectiveReal(pc, nu, 2);
  PetscTryMethod(pc, "PCFieldSplitSetGKBNu_C", (PC, PetscReal), (pc, nu));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCFieldSplitSetGKBNu_FieldSplit(PC pc, PetscReal nu)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

  PetscFunctionBegin;
  jac->gkbnu = nu;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc, PCCompositeType type)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

  PetscFunctionBegin;
  jac->type = type;
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", NULL));

  if (type == PC_COMPOSITE_SCHUR) {
    pc->ops->apply          = PCApply_FieldSplit_Schur;
    pc->ops->applytranspose = PCApplyTranspose_FieldSplit_Schur;
    pc->ops->matapply       = NULL;
    pc->ops->view           = PCView_FieldSplit_Schur;
    pc->ops->setuponblocks  = PCSetUpOnBlocks_FieldSplit_Schur;

    PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit_Schur));
    PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", PCFieldSplitSetSchurPre_FieldSplit));
    PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", PCFieldSplitGetSchurPre_FieldSplit));
    PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", PCFieldSplitSetSchurFactType_FieldSplit));
    PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", PCFieldSplitSetSchurScale_FieldSplit));
  } else if (type == PC_COMPOSITE_GKB) {
    pc->ops->apply          = PCApply_FieldSplit_GKB;
    pc->ops->applytranspose = NULL;
    pc->ops->matapply       = NULL;
    pc->ops->view           = PCView_FieldSplit_GKB;
    pc->ops->setuponblocks  = PCSetUpOnBlocks_FieldSplit_GKB;

    PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit));
    PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", PCFieldSplitSetGKBTol_FieldSplit));
    PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", PCFieldSplitSetGKBMaxit_FieldSplit));
    PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", PCFieldSplitSetGKBNu_FieldSplit));
    PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", PCFieldSplitSetGKBDelay_FieldSplit));
  } else {
    pc->ops->apply          = PCApply_FieldSplit;
    pc->ops->applytranspose = PCApplyTranspose_FieldSplit;
    pc->ops->matapply       = PCMatApply_FieldSplit;
    pc->ops->view           = PCView_FieldSplit;
    pc->ops->setuponblocks  = PCSetUpOnBlocks_FieldSplit;

    PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc, PetscInt bs)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

  PetscFunctionBegin;
  PetscCheck(bs >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Blocksize must be positive, you gave %" PetscInt_FMT, bs);
  PetscCheck(jac->bs <= 0 || jac->bs == bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot change fieldsplit blocksize from %" PetscInt_FMT " to %" PetscInt_FMT " after it has been set", jac->bs, bs);
  jac->bs = bs;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCSetCoordinates_FieldSplit(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
{
  PC_FieldSplit    *jac           = (PC_FieldSplit *)pc->data;
  PC_FieldSplitLink ilink_current = jac->head;
  IS                is_owned;

  PetscFunctionBegin;
  jac->coordinates_set = PETSC_TRUE; // Internal flag
  PetscCall(MatGetOwnershipIS(pc->mat, &is_owned, NULL));

  while (ilink_current) {
    // For each IS, embed it to get local coords indces
    IS              is_coords;
    PetscInt        ndofs_block;
    const PetscInt *block_dofs_enumeration; // Numbering of the dofs relevant to the current block

    // Setting drop to true for safety. It should make no difference.
    PetscCall(ISEmbed(ilink_current->is, is_owned, PETSC_TRUE, &is_coords));
    PetscCall(ISGetLocalSize(is_coords, &ndofs_block));
    PetscCall(ISGetIndices(is_coords, &block_dofs_enumeration));

    // Allocate coordinates vector and set it directly
    PetscCall(PetscMalloc1(ndofs_block * dim, &ilink_current->coords));
    for (PetscInt dof = 0; dof < ndofs_block; ++dof) {
      for (PetscInt d = 0; d < dim; ++d) (ilink_current->coords)[dim * dof + d] = coords[dim * block_dofs_enumeration[dof] + d];
    }
    ilink_current->dim   = dim;
    ilink_current->ndofs = ndofs_block;
    PetscCall(ISRestoreIndices(is_coords, &block_dofs_enumeration));
    PetscCall(ISDestroy(&is_coords));
    ilink_current = ilink_current->next;
  }
  PetscCall(ISDestroy(&is_owned));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitSetType - Sets the type, `PCCompositeType`, of a `PCFIELDSPLIT`

  Collective

  Input Parameters:
+ pc   - the preconditioner context
- type - `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE` (default), `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`,
         `PC_COMPOSITE_GKB`

  Options Database Key:
. -pc_fieldsplit_type <one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type

  Level: intermediate

.seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCCompositeType`, `PCCompositeGetType()`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`,
          `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`, `PCFieldSplitSetSchurFactType()`
@*/
PetscErrorCode PCFieldSplitSetType(PC pc, PCCompositeType type)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscTryMethod(pc, "PCFieldSplitSetType_C", (PC, PCCompositeType), (pc, type));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitGetType - Gets the type, `PCCompositeType`, of a `PCFIELDSPLIT`

  Not collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameter:
. type - `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE` (default), `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`

  Level: intermediate

.seealso: [](sec_block_matrices), `PC`, `PCCompositeSetType()`, `PCFIELDSPLIT`, `PCCompositeType`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`,
          `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`
@*/
PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscAssertPointer(type, 2);
  *type = jac->type;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitSetDMSplits - Flags whether `DMCreateFieldDecomposition()` should be used to define the splits in a `PCFIELDSPLIT`, whenever possible.

  Logically Collective

  Input Parameters:
+ pc  - the preconditioner context
- flg - boolean indicating whether to use field splits defined by the `DM`

  Options Database Key:
. -pc_fieldsplit_dm_splits <bool> - use the field splits defined by the `DM`

  Level: intermediate

  Developer Note:
  The name should be `PCFieldSplitSetUseDMSplits()`, similar change to options database

.seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitGetDMSplits()`, `DMCreateFieldDecomposition()`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
@*/
PetscErrorCode PCFieldSplitSetDMSplits(PC pc, PetscBool flg)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
  PetscBool      isfs;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscValidLogicalCollectiveBool(pc, flg, 2);
  PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
  if (isfs) jac->dm_splits = flg;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitGetDMSplits - Returns flag indicating whether `DMCreateFieldDecomposition()` should be used to define the splits in a `PCFIELDSPLIT`, whenever possible.

  Logically Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameter:
. flg - boolean indicating whether to use field splits defined by the `DM`

  Level: intermediate

  Developer Note:
  The name should be `PCFieldSplitGetUseDMSplits()`

.seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetDMSplits()`, `DMCreateFieldDecomposition()`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
@*/
PetscErrorCode PCFieldSplitGetDMSplits(PC pc, PetscBool *flg)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
  PetscBool      isfs;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscAssertPointer(flg, 2);
  PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
  if (isfs) {
    if (flg) *flg = jac->dm_splits;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitGetDetectSaddlePoint - Returns flag indicating whether `PCFIELDSPLIT` will attempt to automatically determine fields based on zero diagonal entries.

  Logically Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameter:
. flg - boolean indicating whether to detect fields or not

  Level: intermediate

.seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetDetectSaddlePoint()`
@*/
PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc, PetscBool *flg)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

  PetscFunctionBegin;
  *flg = jac->detect;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCFieldSplitSetDetectSaddlePoint - Sets flag indicating whether `PCFIELDSPLIT` will attempt to automatically determine fields based on zero diagonal entries.

  Logically Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameter:
. flg - boolean indicating whether to detect fields or not

  Options Database Key:
. -pc_fieldsplit_detect_saddle_point <bool> - detect and use the saddle point

  Level: intermediate

  Note:
  Also sets the split type to `PC_COMPOSITE_SCHUR` (see `PCFieldSplitSetType()`) and the Schur preconditioner type to `PC_FIELDSPLIT_SCHUR_PRE_SELF` (see `PCFieldSplitSetSchurPre()`).

.seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitGetDetectSaddlePoint()`, `PCFieldSplitSetType()`, `PCFieldSplitSetSchurPre()`, `PC_FIELDSPLIT_SCHUR_PRE_SELF`
@*/
PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc, PetscBool flg)
{
  PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

  PetscFunctionBegin;
  jac->detect = flg;
  if (jac->detect) {
    PetscCall(PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR));
    PetscCall(PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_SELF, NULL));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*MC
  PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
  collections of variables (that may overlap) called fields or splits. Each field often represents a different continuum variable
  represented on a grid, such as velocity, pressure, or temperature.
  In the literature these are sometimes called block preconditioners; but should not be confused with `PCBJACOBI`.
  See [the users manual section on "Solving Block Matrices"](sec_block_matrices) for more details.

  Options Database Keys:
+   -pc_fieldsplit_%d_fields <a,b,..>                                                - indicates the fields to be used in the `%d`'th split
.   -pc_fieldsplit_default                                                           - automatically add any fields to additional splits that have not
                                                                                       been supplied explicitly by `-pc_fieldsplit_%d_fields`
.   -pc_fieldsplit_block_size <bs>                                                   - size of block that defines fields (i.e. there are bs fields)
                                                                                       when the matrix is not of `MatType` `MATNEST`
.   -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur,gkb> - type of relaxation or factorization splitting
.   -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full>                     - default is `a11`; see `PCFieldSplitSetSchurPre()`
.   -pc_fieldsplit_schur_fact_type <diag,lower,upper,full>                           - set factorization type when using `-pc_fieldsplit_type schur`;
                                                                                       see `PCFieldSplitSetSchurFactType()`
.   -pc_fieldsplit_dm_splits <true,false> (default is true)                          - Whether to use `DMCreateFieldDecomposition()` for splits
-   -pc_fieldsplit_detect_saddle_point                                               - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver

  Options prefixes for inner solvers when using the Schur complement preconditioner are `-fieldsplit_0_` and `-fieldsplit_1_` .
  The options prefix for the inner solver when using the Golub-Kahan biadiagonalization preconditioner is `-fieldsplit_0_`
  For all other solvers they are `-fieldsplit_%d_` for the `%d`'th field; use `-fieldsplit_` for all fields.

  To set options on the solvers for all blocks, prepend `-fieldsplit_` to all the `PC`
  options database keys. For example, `-fieldsplit_pc_type ilu` `-fieldsplit_pc_factor_levels 1`.

  To set the options on the solvers separate for each block call `PCFieldSplitGetSubKSP()`
  and set the options directly on the resulting `KSP` object

  Level: intermediate

  Notes:
  Use `PCFieldSplitSetFields()` to set splits defined by "strided" entries or with a `MATNEST` and `PCFieldSplitSetIS()`
  to define a split by an arbitrary collection of entries.

  If no splits are set, the default is used. If a `DM` is associated with the `PC` and it supports
  `DMCreateFieldDecomposition()`, then that is used for the default. Otherwise if the matrix is not `MATNEST`, the splits are defined by entries strided by bs,
  beginning at 0 then 1, etc to bs-1. The block size can be set with `PCFieldSplitSetBlockSize()`,
  if this is not called the block size defaults to the blocksize of the second matrix passed
  to `KSPSetOperators()`/`PCSetOperators()`.

  For the Schur complement preconditioner if
  ```{math}
    J = \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{10} & A_{11} \end{array}\right]
  ```

  the preconditioner using `full` factorization is logically
  ```{math}
    \left[\begin{array}{cc} I & -\text{ksp}(A_{00}) A_{01} \\ 0 & I \end{array}\right] \left[\begin{array}{cc} \text{ksp}(A_{00}) & 0 \\ 0 & \text{ksp}(S) \end{array}\right] \left[\begin{array}{cc} I & 0 \\ -A_{10} \text{ksp}(A_{00}) & I \end{array}\right]
      ```
  where the action of $\text{ksp}(A_{00})$ is applied using the `KSP` solver with prefix `-fieldsplit_0_`.  $S$ is the Schur complement
  ```{math}
     S = A_{11} - A_{10} \text{ksp}(A_{00}) A_{01}
  ```
  which is usually dense and not stored explicitly.  The action of $\text{ksp}(S)$ is computed using the `KSP` solver with prefix `-fieldsplit_splitname_` (where `splitname`
  was given in providing the SECOND split or 1 if not given). Accordingly, if using `PCFieldSplitGetSubKSP()`, the array of sub-`KSP` contexts will hold two `KSP`s: at its
  0th index, the `KSP` associated with `-fieldsplit_0_`, and at its 1st index, the `KSP` corresponding to `-fieldsplit_1_`.
  By default, $A_{11}$ is used to construct a preconditioner for $S$, use `PCFieldSplitSetSchurPre()` for all the possible ways to construct the preconditioner for $S$.

  The factorization type is set using `-pc_fieldsplit_schur_fact_type <diag, lower, upper, full>`. `full` is shown above,
  `diag` gives
  ```{math}
    \left[\begin{array}{cc} \text{ksp}(A_{00}) & 0 \\  0 & -\text{ksp}(S) \end{array}\right]
  ```
  Note that, slightly counter intuitively, there is a negative in front of the $\text{ksp}(S)$  so that the preconditioner is positive definite. For SPD matrices $J$, the sign flip
  can be turned off with `PCFieldSplitSetSchurScale()` or by command line `-pc_fieldsplit_schur_scale 1.0`. The `lower` factorization is the inverse of
  ```{math}
    \left[\begin{array}{cc} A_{00} & 0 \\  A_{10} & S \end{array}\right]
  ```
  where the inverses of $A_{00}$ and $S$ are applied using `KSP`s. The upper factorization is the inverse of
  ```{math}
    \left[\begin{array}{cc} A_{00} & A_{01} \\  0 & S \end{array}\right]
  ```
  where again the inverses of $A_{00}$ and $S$ are applied using `KSP`s.

  If only one set of indices (one `IS`) is provided with `PCFieldSplitSetIS()` then the complement of that `IS`
  is used automatically for a second submatrix.

  The fieldsplit preconditioner cannot currently be used with the `MATBAIJ` or `MATSBAIJ` data formats if the blocksize is larger than 1.
  Generally it should be used with the `MATAIJ` or `MATNEST` `MatType`

  The forms of these preconditioners are closely related, if not identical, to forms derived as "Distributive Iterations", see,
  for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling {cite}`wesseling2009`.
  One can also use `PCFIELDSPLIT` inside a smoother resulting in "Distributive Smoothers".

  See "A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations" {cite}`elman2008tcp`.

  The Constrained Pressure Preconditioner (CPR) can be implemented using `PCCOMPOSITE` with `PCGALERKIN`. CPR first solves an $R A P$ subsystem, updates the
  residual on all variables (`PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)`), and then applies a simple ILU like preconditioner on all the variables.

  The generalized Golub-Kahan bidiagonalization preconditioner (GKB) can be applied to symmetric $2 \times 2$ block matrices of the shape
  ```{math}
    \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{01}' & 0 \end{array}\right]
  ```
  with $A_{00}$ positive semi-definite. The implementation follows {cite}`arioli2013`. Therein, we choose $N := 1/\nu * I$ and the $(1,1)$-block of the matrix is modified to $H = _{A00} + \nu*A_{01}*A_{01}'$.
  A linear system $Hx = b$ has to be solved in each iteration of the GKB algorithm. This solver is chosen with the option prefix `-fieldsplit_0_`.

  Some `PCFIELDSPLIT` variants are called physics-based preconditioners, since the preconditioner takes into account the underlying physics of the
  problem. But this nomenclature is not well-defined.

  Developer Note:
  The Schur complement functionality of `PCFIELDSPLIT` should likely be factored into its own `PC` thus simplifying the implementation of the preconditioners and their
  user API.

.seealso: [](sec_block_matrices), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCLSC`,
          `PCFieldSplitGetSubKSP()`, `PCFieldSplitSchurGetSubKSP()`, `PCFieldSplitSetFields()`,
          `PCFieldSplitSetType()`, `PCFieldSplitSetIS()`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitSetSchurFactType()`,
          `MatSchurComplementSetAinvType()`, `PCFieldSplitSetSchurScale()`, `PCFieldSplitSetDetectSaddlePoint()`
M*/

PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
{
  PC_FieldSplit *jac;

  PetscFunctionBegin;
  PetscCall(PetscNew(&jac));

  jac->bs                 = -1;
  jac->type               = PC_COMPOSITE_MULTIPLICATIVE;
  jac->schurpre           = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
  jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
  jac->schurscale         = -1.0;
  jac->dm_splits          = PETSC_TRUE;
  jac->gkbtol             = 1e-5;
  jac->gkbdelay           = 5;
  jac->gkbnu              = 1;
  jac->gkbmaxit           = 100;

  pc->data = (void *)jac;

  pc->ops->setup           = PCSetUp_FieldSplit;
  pc->ops->reset           = PCReset_FieldSplit;
  pc->ops->destroy         = PCDestroy_FieldSplit;
  pc->ops->setfromoptions  = PCSetFromOptions_FieldSplit;
  pc->ops->applyrichardson = NULL;

  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSchurGetSubKSP_C", PCFieldSplitSchurGetSubKSP_FieldSplit));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetFields_C", PCFieldSplitSetFields_FieldSplit));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", PCFieldSplitSetIS_FieldSplit));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetType_C", PCFieldSplitSetType_FieldSplit));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetBlockSize_C", PCFieldSplitSetBlockSize_FieldSplit));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitRestrictIS_C", PCFieldSplitRestrictIS_FieldSplit));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_FieldSplit));

  /* Initialize function pointers */
  PetscCall(PCFieldSplitSetType(pc, jac->type));
  PetscFunctionReturn(PETSC_SUCCESS);
}