#ifdef PETSC_RCS_HEADER
static char vcid[] = "$Id: schur.c,v 1.15 2000/10/04 02:16:21 knepley Exp $";
#endif
/*
  Defines the Schur preconditioner due to Ahmed Sameh for AIJ matrices
*/
#include "src/sles/pc/pcimpl.h"   /*I "pc.h" I*/
#include "schur.h"

extern int GridResetConstrainedMultiply_Private(Grid, GMat);

/*
  Math:

  There are several possible options for this type of preconditioner. I call it
  a Schur complement preconditioner because we split the problem into blocks

        /  A   B \
        \ B^T  0 /

  and generally must provide a mechanism for solving (preconditioning) the
  Schur complement B^T A^{-1} B at each itertion of the overall preconditioner.
  We could include Uzawa type preconditioners, but since they are so simple, I
  have already abstracted them using GMatCreateUzawa() and related functions.
    Here I am concerned with preconditioners for the entire saddle point problem
  above, the simplest being a block diagonal (Jacobi)

        / \hat A      0   \
        \    0    -\hat S /

  or the slightly more sophisticated block triangular

        / \hat A      B   \
        \    0    -\hat S /

  culminating in the full block LU factorization

        / \hat A  0 \ / I   \hat A^{-1} B \ 
        \   B^T   I / \ 0      -\hat S    /

  These all depend on solving the Schur complement S, so we need an effective
  preconditioner for it. Elman, Wathen and Silvester propose the "BFBT"
  preconditioner which essentially projects the problem onto the pressure space
  and then applies the convection diffusion operator. So we have

        (B^T A^{-1} B)^{-1} = (B^T B)^{-1} B^T A B (B^T B)^{-1}
                            = A_p (B^T B)^{-1}

  where in the last step they assume a commutation relation of the form

        A B  = B A_p

  so that the momentum equations on the pressure mesh are equivalent in some
  sense to those on the velocity mesh. In addition, they recommend several cycles
  of multigrid for \hat A^{-1}.
    Sameh instead proposes to precondition S with the pressure lapalcian L_p and
  perform fast solves with \hat A using the balanced method. I do not understand
  why Wathen, et.al. do not do the full block LU instead of the triangular thing,
  but we can test this. We can afford to do an iterative method for S since
  applications of \hat A^{-1} are cheap, but this may be compared with the sparse
  factorizations of Saad.
    Lastly, it must be noted that you can solve a preconditioned system for \hat A
  or \hat S, or merely apply the preconditioners. The cost of solving must be
  weighed against the extra outer iterations incured by the weaker approximation.

  Programming:

    We allow several modes for this preconditioner since parts of the matrix may
  not be available, e.g. if the matrix is given as implicitly P^T A P. The flag
  explicitConstraints signals that a constrained matrix has been explicitly formed.
  Otherwise, we use GridGetConstraintMatrix() to retrieve the projector P onto the
  constraint space.
    Also, we may have new variables introduced by constraints. These will be present
  already in explicit matrices, but are not yet included in the implicit formulation.
  Thus, we have functions which add this contribution to a given vector, with functions
  like GVecEvaluateJacobianConstrained(). This function operates in the projected space
  so that we get an action like

        P^T A P + G

  where G incorporates the contribution of the new fields.
*/

/*--------------------------------------------- Public Functions ----------------------------------------------------*/
#undef  __FUNCT__
#define __FUNCT__ "PCSchurSetGradientOperator"
/*@C PCSchurSetGradientOperator
	This function sets the operator and its transpose, which define the
  constraint manifold in the saddle-point problem. For instance, in
  the Stokes problem, this operator is the gradient, and its transpose
  is the divergence. It must be called prior to PCSetUp().

  Collective on PC

  Input Parameters:
+ pc     - The preconditioner context
. gradOp - The gradient operator
- divOp  - The divergence operator

  Level: intermediate

.keywords schur, set, gradient
.seealso PCSchurGetIterationNumber()
@*/
int PCSchurSetGradientOperator(PC pc, int gradOp, int divOp)
{
  PC_Schur *s;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_COOKIE);
  s = (PC_Schur *) pc->data;
  s->gradOp = gradOp;
  s->divOp  = divOp;
  PetscFunctionReturn(0);
}
#undef  __FUNCT__
#define __FUNCT__ "PCSchurGetIterationNumber"
/*@C
  PCSchurGetIterationNumber - Get the iteration numbers since the solve was started.

  Not collective

  Input Parameter:
. pc       - The preconditioner context

  Output
+ its      - The total number of iterations used to solve A
- schurIts - The number of iterations used to solve the Schur complement

  Level: intermediate

.keywords schur, iteration
.seealso PCSchurSetGradientOperator()
@*/
int PCSchurGetIterationNumber(PC pc, int *its, int *schurIts)
{
  PC_Schur *s;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_COOKIE);
  PetscValidPointer(its);
  PetscValidPointer(schurIts);
  s = (PC_Schur *) pc->data;
  *its      = s->iter;
  *schurIts = s->schurIter;
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCSchurInnerMonitor"
/*@C
  PCSchurInnerMonitor - Print the residual norm of the solver for the inner system A.

  Collective on KSP

  Input Parameters:
. ksp   - The KSP context
. n     - The iteration number
. rnorm - The 2-norm of residual
. ctx   - The Schur context

  Level: intermediate

.keywords: KSP, schur, monitor, norm
.seealso: KSPSetMonitor(), PCSchurMonitor()
@*/
int PCSchurInnerMonitor(KSP ksp, int n, PetscReal rnorm, void *ctx)
{
  MPI_Comm comm;
  int      ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ksp, KSP_COOKIE);
  ierr = PetscObjectGetComm((PetscObject) ksp, &comm);                                                   CHKERRQ(ierr);
  PetscPrintf(comm, "  inner iter = %d, residual norm %g \n", n, rnorm);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCSchurMonitor"
/*@C
  PCSchurMonitor - Print the residual norm of the solver for the Schur complement system S.

  Collective on KSP

  Input Parameters:
. ksp   - The KSP context
. n     - The iteration number
. rnorm - The 2-norm of residual
. ctx   - The Schur context

  Level: intermediate

.keywords: KSP, schur, monitor, norm
.seealso: KSPSetMonitor(), PCSchurInnerMonitor()
@*/
int PCSchurMonitor(KSP ksp, int n, PetscReal rnorm, void *ctx)
{
  MPI_Comm comm;
  int      ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ksp, KSP_COOKIE);
  ierr = PetscObjectGetComm((PetscObject) ksp, &comm);                                                   CHKERRQ(ierr);
  PetscPrintf(comm, "  schur iter = %d, residual norm %g \n", n, rnorm);
  PetscFunctionReturn(0);
}


#undef __FUNCT__  
#define __FUNCT__ "PCSchurSolveMonitor"
/*@C
  PCSchurSolveMonitor - Print the number of iterates for each inner solve in the Schur complement system.

  Collective on KSP

  Input Parameters:
. ksp   - The KSP context
. n     - The iteration number
. rnorm - The 2-norm of residual
. ctx   - The Schur context

  Level: intermediate

.keywords: KSP, schur, monitor, norm
.seealso: KSPSetMonitor(), PCSchurMonitor(), PCSchurInnerMonitor()
@*/
int PCSchurSolveMonitor(KSP ksp, int n, PetscReal rnorm, void *ctx)
{
  PC_Schur *s = (PC_Schur *) ctx;
  KSP       innerKSP;
  int       its;
  int       ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ksp, KSP_COOKIE);
  ierr = SLESGetKSP(s->sles, &innerKSP);                                                                 CHKERRQ(ierr);
  ierr = KSPGetIterationNumber(innerKSP, &its);                                                          CHKERRQ(ierr);
  if (n > 0)
    PetscLogInfo(ksp, "PCSchur: %d iterations in A solve %d\n", its, n);
  PetscFunctionReturn(0);
}

/*--------------------------------------------- Private Functions ---------------------------------------------------*/
#undef __FUNCT__  
#define __FUNCT__ "PCDestroyStructures_Schur"
static int PCDestroyStructures_Schur(PC pc)
{
  PC_Schur *s = (PC_Schur *) pc->data;
  int       ierr;

  PetscFunctionBegin;
  /* Destroy operator structures */
  ierr = PetscFree(s->momOps);                                                                           CHKERRQ(ierr);
  ierr = PetscFree(s->momOpIsALE);                                                                       CHKERRQ(ierr);
  ierr = PetscFree(s->momOpAlphas);                                                                      CHKERRQ(ierr);

  /* Destroy variable orderings */
  ierr = VarOrderingDestroy(s->sOrder);                                                                  CHKERRQ(ierr);
  ierr = VarOrderingDestroy(s->tOrder);                                                                  CHKERRQ(ierr);
  ierr = LocalVarOrderingDestroy(s->sLocOrder);                                                          CHKERRQ(ierr);
  ierr = LocalVarOrderingDestroy(s->tLocOrder);                                                          CHKERRQ(ierr);

  /* Destroy matrices and reorderings */
  ierr = GMatDestroyUzawa(s->S);                                                                         CHKERRQ(ierr);
  ierr = GMatDestroy(s->A);                                                                              CHKERRQ(ierr);
  if (s->sparseA != PETSC_NULL) {
    ierr = GMatDestroy(s->sparseA);                                                                      CHKERRQ(ierr);
  }
  ierr = GMatDestroy(s->B);                                                                              CHKERRQ(ierr);
  if (s->lap != PETSC_NULL) {
    ierr = GMatDestroy(s->lap);                                                                          CHKERRQ(ierr);
  }
  if (s->rowPerm != PETSC_NULL) {
    ierr = ISDestroy(s->rowPerm);                                                                        CHKERRQ(ierr);
    ierr = ISDestroy(s->colPerm);                                                                        CHKERRQ(ierr);
  }
  /* Destroy field structures */
  ierr = VecDestroy(s->u);                                                                               CHKERRQ(ierr);
  ierr = VecDestroy(s->uNew);                                                                            CHKERRQ(ierr);
  ierr = VecDestroy(s->uNew2);                                                                           CHKERRQ(ierr);
  ierr = VecDestroy(s->p);                                                                               CHKERRQ(ierr);
  ierr = VecDestroy(s->pNew);                                                                            CHKERRQ(ierr);
  ierr = VecScatterDestroy(s->uScatter);                                                                 CHKERRQ(ierr);
  ierr = VecScatterDestroy(s->pScatter);                                                                 CHKERRQ(ierr);
  if ((s->isConstrained == PETSC_TRUE) && (s->explicitConstraints == PETSC_FALSE)) {
    ierr = VecDestroy(s->projX);                                                                         CHKERRQ(ierr);
    ierr = VecDestroy(s->projY);                                                                         CHKERRQ(ierr);
    ierr = VecDestroy(s->projZ);                                                                         CHKERRQ(ierr);
  }

  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "PCValidQ_Schur"
int PCValidQ_Schur(PC pc)
{
  PC_Schur *s = (PC_Schur *) pc->data;

  PetscFunctionBegin;
  if (pc->setupcalled < 2)
    PetscFunctionReturn(1);

  PetscValidPointer(s);
  /* Check dimensions */

  /* Check fields */
  PetscValidHeaderSpecific(s->u,         VEC_COOKIE);
  PetscValidHeaderSpecific(s->uNew,      VEC_COOKIE);
  PetscValidHeaderSpecific(s->uNew2,     VEC_COOKIE);
  PetscValidHeaderSpecific(s->p,         VEC_COOKIE);
  PetscValidHeaderSpecific(s->pNew,      VEC_COOKIE);
  if ((s->isConstrained == PETSC_TRUE) && (s->explicitConstraints == PETSC_FALSE)) {
    PetscValidHeaderSpecific(s->projX,   VEC_COOKIE);
    PetscValidHeaderSpecific(s->projY,   VEC_COOKIE);
    PetscValidHeaderSpecific(s->projZ,   VEC_COOKIE);
  }
  PetscValidHeaderSpecific(s->uScatter,  VEC_SCATTER_COOKIE);
  PetscValidHeaderSpecific(s->pScatter,  VEC_SCATTER_COOKIE);

  /* Check inner solvers */
  PetscValidHeaderSpecific(s->sles,      SLES_COOKIE);
  PetscValidHeaderSpecific(s->schurSles, SLES_COOKIE);

  /* Check matrices */
  PetscValidHeaderSpecific(s->A,         MAT_COOKIE);
  PetscValidHeaderSpecific(s->B,         MAT_COOKIE);
  PetscValidHeaderSpecific(s->S,         MAT_COOKIE);
  if (s->sparseA != PETSC_NULL) {
    PetscValidHeaderSpecific(s->sparseA, MAT_COOKIE);
  }

  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCDebug_Schur"
static int PCDebug_Schur(PC pc)
{
  PC_Schur    *s = (PC_Schur *) pc->data;
  Vec          x, y;
  Mat          P;
  Vec          xx;
  Vec          yy;
  PetscScalar *xArray;
  PetscScalar  zero = 0.0, one = 1.0, minusOne = -1.0;
  PetscReal    uNorm, pNorm;
  int          row, size;
  int          ierr;

  PetscFunctionBegin;
  ierr = VecDuplicate(pc->vec, &x);                                                                       CHKERRQ(ierr);
  ierr = VecDuplicate(pc->vec, &y);                                                                       CHKERRQ(ierr);
  if (s->explicitConstraints == PETSC_FALSE) {
    ierr = GridGetConstraintMatrix(s->grid, &P);                                                          CHKERRQ(ierr);
    xx   = s->projX;
    yy   = s->projY;
  } else {
    xx   = x;
    yy   = y;
  }

  /* Check the matrix-vector product */
  ierr = VecGetSize(x, &size);                                                                            CHKERRQ(ierr);
  ierr = VecGetArray(x, &xArray);                                                                         CHKERRQ(ierr);
  for(row = 0; row < size; row++) {
    /* Start with e_row */
    ierr = VecSet(&zero, x);                                                                              CHKERRQ(ierr);
    xArray[row] = 1.0;
    /* \hat A x */
    ierr = MatMult(pc->mat, x, y);                                                                        CHKERRQ(ierr);
    if (s->explicitConstraints == PETSC_FALSE) {
      /* Project vectors */
      ierr = MatMult(P, x, xx);                                                                           CHKERRQ(ierr);
      ierr = MatMult(P, y, yy);                                                                           CHKERRQ(ierr);
    }
    /* New \hat A x */
    ierr = VecScatterBegin(xx, s->u, INSERT_VALUES, SCATTER_FORWARD, s->uScatter);                        CHKERRQ(ierr);
    ierr = VecScatterEnd(xx, s->u, INSERT_VALUES, SCATTER_FORWARD, s->uScatter);                          CHKERRQ(ierr);
    ierr = VecScatterBegin(xx, s->p, INSERT_VALUES, SCATTER_FORWARD, s->pScatter);                        CHKERRQ(ierr);
    ierr = VecScatterEnd(xx, s->p, INSERT_VALUES, SCATTER_FORWARD, s->pScatter);                          CHKERRQ(ierr);
    if (s->colPerm != PETSC_NULL) {
      ierr = VecPermute(s->u, s->colPerm, PETSC_FALSE);                                                   CHKERRQ(ierr);
    }
    /* uNew2 = \hat A u + B p, pNew = B u */
    ierr = MatMult(s->A, s->u, s->uNew);                                                                  CHKERRQ(ierr);
    ierr = MatMult(s->B, s->p, s->uNew2);                                                                 CHKERRQ(ierr);
    ierr = VecAXPY(&one, s->uNew2, s->uNew);                                                              CHKERRQ(ierr);
    ierr = MatMultTranspose(s->B, s->u, s->pNew);                                                         CHKERRQ(ierr);
    if (s->rowPerm != PETSC_NULL) {
      ierr = VecPermute(s->uNew, s->rowPerm, PETSC_TRUE);                                                 CHKERRQ(ierr);
    }
    if (s->explicitConstraints == PETSC_FALSE) {
      /* Add in contribution from new variables */
      /* ierr = GVecEvaluateJacobianConstrained(xx, xx);                                                    CHKERRQ(ierr); */
      ierr = VecScatterBegin(xx, s->uNew, ADD_VALUES, SCATTER_FORWARD, s->uScatter);                      CHKERRQ(ierr);
      ierr = VecScatterEnd(xx, s->uNew, ADD_VALUES, SCATTER_FORWARD, s->uScatter);                        CHKERRQ(ierr);
      ierr = VecScatterBegin(xx, s->pNew, ADD_VALUES, SCATTER_FORWARD, s->pScatter);                      CHKERRQ(ierr);
      ierr = VecScatterEnd(xx, s->pNew, ADD_VALUES, SCATTER_FORWARD, s->pScatter);                        CHKERRQ(ierr);
    }
    /* Compare field vectors */
    ierr = VecScatterBegin(yy, s->u, INSERT_VALUES, SCATTER_FORWARD, s->uScatter);                        CHKERRQ(ierr);
    ierr = VecScatterEnd(yy, s->u, INSERT_VALUES, SCATTER_FORWARD, s->uScatter);                          CHKERRQ(ierr);
    ierr = VecScatterBegin(yy, s->p, INSERT_VALUES, SCATTER_FORWARD, s->pScatter);                        CHKERRQ(ierr);
    ierr = VecScatterEnd(yy, s->p, INSERT_VALUES, SCATTER_FORWARD, s->pScatter);                          CHKERRQ(ierr);
    ierr = VecAXPY(&minusOne, s->u, s->uNew);                                                             CHKERRQ(ierr);
    ierr = VecAXPY(&minusOne, s->p, s->pNew);                                                             CHKERRQ(ierr);
    ierr = VecNorm(s->uNew, NORM_2, &uNorm);                                                              CHKERRQ(ierr);
    ierr = VecNorm(s->pNew, NORM_2, &pNorm);                                                              CHKERRQ(ierr);
    if ((uNorm > 1.0e-8) || (pNorm > 1.0e-8)) {
      SETERRQ(PETSC_ERR_PLIB, "Invalid field matrix");
    }
  }

  ierr = VecRestoreArray(x, &xArray);                                                                     CHKERRQ(ierr);
  ierr = VecDestroy(x);                                                                                   CHKERRQ(ierr);
  ierr = VecDestroy(y);                                                                                   CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCSetFromOptions_Inner"
int PCSetFromOptions_Inner(PC pc, SLES sles, SLES schurSles)
{
  PC_Schur  *s = (PC_Schur *) pc->data;
  KSP        ksp;
  char      *prefix;
  PetscTruth opt;
  int        ierr;

  PetscFunctionBegin;
  ierr = SLESGetOptionsPrefix(sles, &prefix);                                                             CHKERRQ(ierr);
  ierr = SLESGetKSP(sles, &ksp);                                                                          CHKERRQ(ierr);
  ierr = PetscOptionsHasName(prefix, "-ksp_monitor", &opt);                                               CHKERRQ(ierr);
  if (opt == PETSC_TRUE) {
    ierr = KSPClearMonitor(ksp);                                                                          CHKERRQ(ierr);
    ierr = KSPSetMonitor(ksp, PCSchurInnerMonitor, s, PETSC_NULL);                                        CHKERRQ(ierr);
  }
  ierr = SLESGetOptionsPrefix(schurSles, &prefix);                                                        CHKERRQ(ierr);
  ierr = SLESGetKSP(schurSles, &ksp);                                                                     CHKERRQ(ierr);
  ierr = PetscOptionsHasName(prefix, "-ksp_monitor", &opt);                                               CHKERRQ(ierr);
  if (opt == PETSC_TRUE) {
    ierr = KSPClearMonitor(ksp);                                                                          CHKERRQ(ierr);
    ierr = KSPSetMonitor(ksp, PCSchurMonitor, s, PETSC_NULL);                                             CHKERRQ(ierr);
  }
  ierr = PetscOptionsHasName(prefix, "-ksp_solve_monitor", &opt);                                         CHKERRQ(ierr);
  if (opt == PETSC_TRUE) {
    ierr = KSPClearMonitor(ksp);                                                                          CHKERRQ(ierr);
    ierr = KSPSetMonitor(ksp, PCSchurSolveMonitor, s, PETSC_NULL);                                        CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCSetFromOptions_Schur"
static int PCSetFromOptions_Schur(PC pc)
{
  PC_Schur  *s = (PC_Schur *) pc->data;
  PetscTruth opt;
  int        ierr;

  PetscFunctionBegin;
  /* Enable explicit constraints */
  ierr = PetscOptionsHasName(pc->prefix, "-pc_schur_explicit", &opt);                                     CHKERRQ(ierr);
  if (opt == PETSC_TRUE)
    s->explicitConstraints = PETSC_TRUE;

  /* Enable laplacian preconditioning of the Schur complement */
  ierr = PetscOptionsHasName(pc->prefix, "-pc_schur_lap", &opt);                                          CHKERRQ(ierr);
  if (opt == PETSC_TRUE)
    s->useLaplacian = PETSC_TRUE;

  /* Enable Mathematica support */
  ierr = PetscOptionsHasName(pc->prefix, "-pc_schur_math", &opt);                                         CHKERRQ(ierr);
  if (opt == PETSC_TRUE) {
    s->useMath = PETSC_TRUE;
    ierr = PetscViewerMathematicaOpen(pc->comm, PETSC_DECIDE, PETSC_NULL, PETSC_NULL, &s->mathViewer);    CHKERRQ(ierr);
  }

  /* Setup inner solvers */
  ierr = PCSetFromOptions_Inner(pc, s->sles, s->schurSles);                                               CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCSchurGetMomentumOperators"
static int PCSchurGetMomentumOperators(PC pc)
{
  PC_Schur   *s = (PC_Schur *) pc->data;
  int         numOps;         /* Grid op query: The number of grid matrix operators used for A */
  int         sField, tField; /* Grid op query: The shape and test function fields */
  PetscScalar alpha;          /* Grid op query: The scalar multiplier */
  PetscTruth  isALE;          /* Grid op query: The flag for ALE operators */
  int         op;
  int         ierr;

  PetscFunctionBegin;
  ierr = GridGetNumMatOperators(s->grid, &s->numMomOps);                                                  CHKERRQ(ierr);
  s->numMomOps  -= 2;
  ierr = PetscMalloc(s->numMomOps * sizeof(int),        &s->momOps);                                      CHKERRQ(ierr);
  ierr = PetscMalloc(s->numMomOps * sizeof(PetscTruth), &s->momOpIsALE);                                  CHKERRQ(ierr);
  ierr = PetscMalloc(s->numMomOps * sizeof(double),     &s->momOpAlphas);                                 CHKERRQ(ierr);
  /* Find gradient fields */
  ierr   = GridGetMatOperatorStart(s->grid, &op, &sField, &tField, &alpha, &isALE);                       CHKERRQ(ierr);
  numOps = 0;
  while(op >= 0) {
    if (op == s->gradOp) {
      s->sField      = sField;
      s->tField      = tField;
      s->gradOpAlpha = alpha;
    } else if (op != s->divOp) {
      s->momOps[numOps]      = op;
      s->momOpIsALE[numOps]  = isALE;
      s->momOpAlphas[numOps] = alpha;
      numOps++;
    }
    ierr = GridGetMatOperatorNext(s->grid, &op, &sField, &tField, &alpha, &isALE);                        CHKERRQ(ierr);
  }
  if ((s->sField < 0) || (s->tField < 0)) {
    SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Gradient operator not found in grid.");
  }
  if (numOps != s->numMomOps) {
    SETERRQ(PETSC_ERR_PLIB, "Inconsistent operator setup in grid");
  }

  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCSchurReorderMatrices"
static int PCSchurReorderMatrices(PC pc)
{
  PC_Schur  *s = (PC_Schur *) pc->data;
  Mat        newA, newB;     /* The reordering matrices */
  char       orderType[256]; /* The type of reordering */
  IS         tempIS;         /* The identity ordering for B */
  int        bw;             /* The bandwidth limitation on the reordered matrix */
  PetscReal  frac;           /* The fractional bandwidth limitation on the reordered matrix */
  double     tol;            /* The drop tolerance */
  int        m, n;
  PetscTruth opt, issparse;
  int        ierr;

  PetscFunctionBegin;
  /* Reduce the bandwidth of A */
  ierr = PetscOptionsHasName(pc->prefix, "-pc_schur_reorder", &opt);                                      CHKERRQ(ierr);
  if (opt == PETSC_TRUE) {
    ierr = MatGetOrdering(s->A, MATORDERING_RCM, &s->rowPerm, &s->colPerm);                               CHKERRQ(ierr);
    ierr = PetscOptionsGetString(pc->prefix, "-pc_schur_reorder", orderType, 255, &opt);                  CHKERRQ(ierr);
    if (opt == PETSC_TRUE) {
      ierr = PetscStrcasecmp(orderType, "sparse", &issparse);                                             CHKERRQ(ierr);
      if (issparse) {
        bw   = PETSC_DECIDE;
        ierr = PetscOptionsGetInt(pc->prefix, "-pc_schur_sparsify_bw", &bw, &opt);                        CHKERRQ(ierr);
        frac = 0.0;
        ierr = PetscOptionsGetReal(pc->prefix, "-pc_schur_sparsify_frac", &frac, &opt);                   CHKERRQ(ierr);
        tol  = 0.0;
        ierr = PetscOptionsGetReal(pc->prefix, "-pc_schur_sparsify_tol", &frac, &opt);                    CHKERRQ(ierr);
        ierr = MatPermuteSparsify(s->A, bw, frac, tol, s->rowPerm, s->colPerm, &s->sparseA);              CHKERRQ(ierr);
      }
      ierr = MatPermute(s->A, s->rowPerm, s->colPerm, &newA);                                             CHKERRQ(ierr);
    }
    ierr = MatDestroy(s->A);                                                                              CHKERRQ(ierr);
    ierr = MatGetSize(s->B, &m, &n);                                                                      CHKERRQ(ierr);
    ierr = ISCreateStride(pc->comm, n, 0, 1, &tempIS);                                                    CHKERRQ(ierr);
    ierr = ISSetPermutation(tempIS);                                                                      CHKERRQ(ierr);
    ierr = MatPermute(s->B, s->rowPerm, tempIS, &newB);                                                   CHKERRQ(ierr);
    ierr = ISDestroy(tempIS);                                                                             CHKERRQ(ierr);
    ierr = MatDestroy(s->B);                                                                              CHKERRQ(ierr);
    s->A = newA;
    s->B = newB;
    ierr = PetscObjectCompose((PetscObject) s->A, "Grid", (PetscObject) s->grid);                         CHKERRQ(ierr);
    ierr = PetscObjectCompose((PetscObject) s->B, "Grid", (PetscObject) s->grid);                         CHKERRQ(ierr);
  }

  /* View the new matrices */
  ierr = PetscOptionsHasName(PETSC_NULL, "-pc_schur_view", &opt);                                         CHKERRQ(ierr);
  if (opt == PETSC_TRUE) {
    PetscViewer viewer = PETSC_VIEWER_DRAW_(pc->comm);
    PetscDraw   draw;

    ierr = PetscViewerDrawGetDraw(viewer, 0, &draw);                                                      CHKERRQ(ierr);
    ierr = PetscDrawSetPause(draw, -1);                                                                   CHKERRQ(ierr);
    ierr = MatView(s->A, viewer);                                                                         CHKERRQ(ierr);
    if (s->sparseA != PETSC_NULL) {
      ierr = MatView(s->sparseA, viewer);                                                                 CHKERRQ(ierr);
    }
    ierr = MatView(s->B, viewer);                                                                         CHKERRQ(ierr);
  }

  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCSchurCreateMatrices"
static int PCSchurCreateMatrices(PC pc)
{
  PC_Schur     *s    = (PC_Schur *) pc->data;
  Grid          grid = s->grid;
  VarOrdering   tempOrder;
  int           op;
  int           ierr;

  PetscFunctionBegin;
  ierr = GridIsConstrained(grid, &s->isConstrained);                                                     CHKERRQ(ierr);
  /* Create test function variable orderings */
  ierr = VarOrderingCreateGeneral(grid, 1, &s->tField, &s->tOrder);                                      CHKERRQ(ierr);
  ierr = LocalVarOrderingCreate(grid, 1, &s->tField, &s->tLocOrder);                                     CHKERRQ(ierr);

  /* Create matrices */
  if ((s->isConstrained == PETSC_FALSE) || (s->explicitConstraints == PETSC_FALSE)) {
    /* Create shape function variable orderings */
    ierr = VarOrderingCreateGeneral(grid, 1, &s->sField, &s->sOrder);                                    CHKERRQ(ierr);
    ierr = LocalVarOrderingCreate(grid, 1, &s->sField, &s->sLocOrder);                                   CHKERRQ(ierr);
    ierr = GMatCreateRectangular(grid, s->tOrder, s->tOrder, &s->A);                                     CHKERRQ(ierr);
    ierr = GMatCreateRectangular(grid, s->sOrder, s->tOrder, &s->B);                                     CHKERRQ(ierr);
    if (s->useLaplacian == PETSC_TRUE) {
      ierr = GMatCreateRectangular(grid, s->sOrder, s->sOrder, &s->lap);                                 CHKERRQ(ierr);
    }
  } else {
    /* Constrain test function variable ordering */
    ierr = VarOrderingConstrain(grid, s->tOrder, &tempOrder);                                            CHKERRQ(ierr);
    ierr = VarOrderingDestroy(s->tOrder);                                                                CHKERRQ(ierr);
    s->tOrder = tempOrder;
    /* Create shape function variable orderings */
    ierr = VarOrderingCreateGeneral(grid, 1, &s->sField, &s->sOrder);                                    CHKERRQ(ierr);
    ierr = LocalVarOrderingCreate(grid, 1, &s->sField, &s->sLocOrder);                                   CHKERRQ(ierr);
    ierr = GMatCreateRectangular(grid, s->tOrder, s->tOrder, &s->A);                                     CHKERRQ(ierr);
    ierr = GMatCreateRectangular(grid, s->sOrder, s->tOrder, &s->B);                                     CHKERRQ(ierr);
    if (s->useLaplacian == PETSC_TRUE) {
      ierr = GMatCreateRectangular(grid, s->sOrder, s->sOrder, &s->lap);                                 CHKERRQ(ierr);
    }
  }
  ierr = MatSetOption(s->A, MAT_NEW_NONZERO_ALLOCATION_ERR);                                             CHKERRQ(ierr);
  ierr = MatSetOption(s->B, MAT_NEW_NONZERO_ALLOCATION_ERR);                                             CHKERRQ(ierr);
  if (s->useLaplacian == PETSC_TRUE) {
    ierr = MatSetOption(s->lap, MAT_NEW_NONZERO_ALLOCATION_ERR);                                         CHKERRQ(ierr);
  }

  /* Construct the operators */
  if ((s->isConstrained == PETSC_FALSE) || (s->explicitConstraints == PETSC_FALSE)) {
    for(op = 0; op < s->numMomOps; op++) {
      if (s->momOpIsALE[op] == PETSC_FALSE) {
        ierr = GMatEvaluateOperatorGalerkin(s->A, PETSC_NULL, 1, &s->tField, &s->tField, s->momOps[op], s->momOpAlphas[op],
                                            MAT_FINAL_ASSEMBLY, s);
        CHKERRQ(ierr);
      } else {
        ierr = GMatEvaluateALEOperatorGalerkin(s->A, 1, &s->tField, s->tOrder, s->tLocOrder, &s->tField, s->tOrder,
                                               s->tLocOrder, s->momOps[op], s->momOpAlphas[op],
                                               MAT_FINAL_ASSEMBLY, s);
        CHKERRQ(ierr);
      }
    }
    ierr = GMatEvaluateOperatorGalerkin(s->B, PETSC_NULL, 1, &s->sField, &s->tField, s->gradOp, s->gradOpAlpha, MAT_FINAL_ASSEMBLY, s);
    CHKERRQ(ierr);
    if (s->useLaplacian == PETSC_TRUE) {
      ierr = GMatEvaluateOperatorGalerkin(s->lap, PETSC_NULL, 1, &s->sField, &s->sField, LAPLACIAN, 1.0, MAT_FINAL_ASSEMBLY, s);
      CHKERRQ(ierr);
    }
  } else {
    for(op = 0; op < s->numMomOps; op++) {
      if (s->momOpIsALE[op] == PETSC_FALSE) {
        ierr = GMatEvaluateOperatorGalerkin(s->A, PETSC_NULL, 1, &s->tField, &s->tField, s->momOps[op], s->momOpAlphas[op],
                                            MAT_FINAL_ASSEMBLY, s);
        CHKERRQ(ierr);
      } else {
        ierr = GMatEvaluateALEConstrainedOperatorGalerkin(s->A, 1, &s->tField, s->tOrder, s->tLocOrder, &s->tField, s->tOrder,
                                                          s->tLocOrder, s->momOps[op], s->momOpAlphas[op], MAT_FINAL_ASSEMBLY, s);
        CHKERRQ(ierr);
      }
    }
    ierr = GMatEvaluateNewFields(s->A, 1, &s->tField, s->tOrder, s->tLocOrder, &s->tField, s->tOrder, s->tLocOrder, 1.0,
                                 MAT_FINAL_ASSEMBLY);
    CHKERRQ(ierr);
    ierr = GMatEvaluateOperatorGalerkin(s->B, PETSC_NULL, 1, &s->sField, &s->tField, s->gradOp, s->gradOpAlpha, MAT_FINAL_ASSEMBLY, s);
    CHKERRQ(ierr);
    if (s->useLaplacian == PETSC_TRUE) {
      ierr = GMatEvaluateOperatorGalerkin(s->lap, PETSC_NULL, 1, &s->sField, &s->sField, LAPLACIAN, 1.0, MAT_FINAL_ASSEMBLY, s);
      CHKERRQ(ierr);
    }
    /* Reset matrix multiplication */
    ierr = GridResetConstrainedMultiply_Private(grid, s->A);                                              CHKERRQ(ierr);
    ierr = GridResetConstrainedMultiply_Private(grid, s->B);                                              CHKERRQ(ierr);
    if (s->useLaplacian == PETSC_TRUE) {
      ierr = GridResetConstrainedMultiply_Private(grid, s->lap);                                          CHKERRQ(ierr);
    }
  }

  /* Reduce the bandwidth of A */
  ierr = PCSchurReorderMatrices(pc);                                                                      CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCSchurCreateInnerSLES"
static int PCSchurCreateInnerSLES(PC pc)
{
  PC_Schur *s = (PC_Schur *) pc->data;
  int       ierr;

  PetscFunctionBegin;
  /* Create the momentum solver */
  ierr = SLESCreate(pc->comm, &s->sles);                                                                 CHKERRQ(ierr);
  ierr = SLESAppendOptionsPrefix(s->sles, "pc_schur_inner_");                                            CHKERRQ(ierr);

  /* Create the Schur complement solver */
  ierr = SLESCreate(pc->comm, &s->schurSles);                                                            CHKERRQ(ierr);
  ierr = SLESAppendOptionsPrefix(s->schurSles, "pc_schur_");                                             CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCSchurSetupInnerSLES"
static int PCSchurSetupInnerSLES(PC pc)
{
  PC_Schur *s = (PC_Schur *) pc->data;
  KSP       innerKsp;
  PC        innerPc;
  int       ierr;

  PetscFunctionBegin;
  /* Create the momentum solver */
  ierr = SLESCreate(pc->comm, &s->sles);                                                                  CHKERRQ(ierr);
  if (s->sparseA != PETSC_NULL) {
    ierr = SLESSetOperators(s->sles, s->A, s->sparseA, DIFFERENT_NONZERO_PATTERN);                        CHKERRQ(ierr);
  } else {
    ierr = SLESSetOperators(s->sles, s->A, s->A,       SAME_NONZERO_PATTERN);                             CHKERRQ(ierr);
  }
  ierr = SLESAppendOptionsPrefix(s->sles, "pc_schur_inner_");                                             CHKERRQ(ierr);
  ierr = SLESGetKSP(s->sles, &innerKsp);                                                                  CHKERRQ(ierr);
  ierr = SLESGetPC(s->sles, &innerPc);                                                                    CHKERRQ(ierr);
  ierr = KSPSetType(innerKsp, KSPPREONLY);                                                                CHKERRQ(ierr);
  ierr = PCSetType(innerPc, PCLU);                                                                        CHKERRQ(ierr);
  ierr = SLESSetFromOptions(s->sles);                                                                     CHKERRQ(ierr);

  /* Create the Schur complement solver */
  ierr = SLESCreate(pc->comm, &s->schurSles);                                                             CHKERRQ(ierr);
  ierr = GMatCreateUzawa(s->A, s->B, s->sles, &s->S);                                                     CHKERRQ(ierr);
  ierr = SLESSetOperators(s->schurSles, s->S, s->S, SAME_NONZERO_PATTERN);                                CHKERRQ(ierr);
  ierr = SLESAppendOptionsPrefix(s->schurSles, "pc_schur_");                                              CHKERRQ(ierr);
  ierr = SLESGetKSP(s->schurSles, &innerKsp);                                                             CHKERRQ(ierr);
  ierr = SLESGetPC(s->schurSles, &innerPc);                                                               CHKERRQ(ierr);
  ierr = KSPSetType(innerKsp, KSPGMRES);                                                                  CHKERRQ(ierr);
  if (s->useLaplacian == PETSC_TRUE) {
    ierr = PCSetOperators(innerPc, s->S, s->lap, DIFFERENT_NONZERO_PATTERN);                              CHKERRQ(ierr);
    ierr = PCSetType(innerPc, PCLU);                                                                      CHKERRQ(ierr);
  } else {
    ierr = PCSetType(innerPc, PCNONE);                                                                    CHKERRQ(ierr);
  }
  ierr = SLESSetFromOptions(s->schurSles);                                                                CHKERRQ(ierr);

  /* Setup monitors */
  ierr = PCSetFromOptions_Inner(pc, s->sles, s->schurSles);                                               CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCCreateStructures_Schur"
static int PCCreateStructures_Schur(PC pc)
{
  PC_Schur *s    = (PC_Schur *) pc->data;
  Grid      grid = s->grid;
  int       ierr;

  PetscFunctionBegin;
  /* Create matrices */
  ierr = PCSchurCreateMatrices(pc);                                                                      CHKERRQ(ierr);

  /* Create the work vectors and field scatters */
  ierr = GVecCreateRectangular(grid, s->tOrder, &s->u);                                                  CHKERRQ(ierr);
  ierr = GVecDuplicate(s->u, &s->uNew);                                                                  CHKERRQ(ierr);
  ierr = GVecDuplicate(s->u, &s->uNew2);                                                                 CHKERRQ(ierr);
  ierr = GVecCreateRectangular(grid, s->sOrder, &s->p);                                                  CHKERRQ(ierr);
  ierr = GVecDuplicate(s->p, &s->pNew);                                                                  CHKERRQ(ierr);
  if ((s->isConstrained == PETSC_TRUE) && (s->explicitConstraints == PETSC_FALSE)) {
    ierr = GVecCreate(grid, &s->projX);                                                                  CHKERRQ(ierr);
    ierr = GVecDuplicate(s->projX, &s->projY);                                                           CHKERRQ(ierr);
    ierr = GVecCreateConstrained(grid, &s->projZ);                                                       CHKERRQ(ierr);
    ierr = GVecScatterCreate(s->projX, s->u, &s->uScatter);                                              CHKERRQ(ierr);
    ierr = GVecScatterCreate(s->projX, s->p, &s->pScatter);                                              CHKERRQ(ierr);
  } else {
    ierr = GVecScatterCreate(pc->vec, s->u, &s->uScatter);                                               CHKERRQ(ierr);
    ierr = GVecScatterCreate(pc->vec, s->p, &s->pScatter);                                               CHKERRQ(ierr);
  }

  /* Setup Inner Solvers */
  ierr = PCSchurSetupInnerSLES(pc);                                                                      CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCSetUp_Schur"
static int PCSetUp_Schur(PC pc)
{
  PC_Schur  *s = (PC_Schur *) pc->data;
  PetscTruth opt;
  int        ierr;

  PetscFunctionBegin;
  /* This lets the factorization become stale, that is uses the same preconditioner for several Newton steps */
  if (pc->setupcalled > 0) {
    PetscFunctionReturn(0);
  }
  pc->setupcalled = 2;

  /* Get the grid */
  ierr = GVecGetGrid(pc->vec, &s->grid);                                                                  CHKERRQ(ierr);

  /* Divide up the problem */
  ierr = PCSchurGetMomentumOperators(pc);                                                                 CHKERRQ(ierr);

  /* Create matrices, vector storage, scatters, and inner solvers */
  ierr = PCCreateStructures_Schur(pc);                                                                    CHKERRQ(ierr);

#ifdef PETSC_USE_BOPT_g
  ierr = PetscTrValid(__LINE__, __FUNCT__, __FILE__, __SDIR__);                                            CHKERRQ(ierr);
#endif
  ierr = PetscOptionsHasName(PETSC_NULL, "-pc_schur_debug", &opt);                                        CHKERRQ(ierr);
  if (opt == PETSC_TRUE) {
    ierr = PCDebug_Schur(pc);                                                                             CHKERRQ(ierr);
  }

  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCPreSolve_Schur"
static int PCPreSolve_Schur(PC pc, KSP ksp, Vec x,  Vec rhs)
{
  PC_Schur *s = (PC_Schur *) pc->data;

  PetscFunctionBegin;
  s->iter      = 0;
  s->schurIter = 0;
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCPostSolve_Schur"
static int PCPostSolve_Schur(PC pc, KSP ksp, Vec x, Vec rhs)
{
  PetscFunctionBegin;
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCApply_Schur"
static int PCApply_Schur(PC pc, Vec x, Vec y)
{
  PC_Schur   *s  = (PC_Schur *) pc->data;
  Vec         xx = x;
  Vec         yy = y;
  Vec         z  = s->projZ;
  Vec         tempSol;
  Mat         P, invP;
  PetscScalar minusOne = -1.0;
  int         its, schurIts;
  PetscTruth  opt;
  int         ierr;

  PetscFunctionBegin;
  ierr = PetscOptionsHasName(PETSC_NULL, "-pc_schur_test", &opt);                                         CHKERRQ(ierr);
  if (opt == PETSC_TRUE) {
    /* Set the input vector to Ax */
    ierr = VecDuplicate(x, &tempSol);                                                                     CHKERRQ(ierr);
    ierr = VecCopy(x, tempSol);                                                                           CHKERRQ(ierr);
    ierr = MatMult(pc->mat, x, y);                                                                        CHKERRQ(ierr);
    ierr = VecCopy(y, x);                                                                                 CHKERRQ(ierr);
  }

  /* Project: Apply (P^+)^T = P (P^T P)^{-1} */
  if ((s->isConstrained == PETSC_TRUE) && (s->explicitConstraints == PETSC_FALSE)) {
    xx   = s->projX;
    yy   = s->projY;
    ierr = GridGetConstraintMatrix(s->grid, &P);                                                          CHKERRQ(ierr);
    ierr = GridGetConstraintInverse(s->grid, &invP);                                                      CHKERRQ(ierr);
    ierr = MatMult(invP, x, z);                                                                           CHKERRQ(ierr);
    ierr = MatMult(P, z, xx);                                                                             CHKERRQ(ierr);
  }

  /*            /    A^{-1}    0      0   \ / u \   /     A^{-1} u     \   / w_u \
     L^{-1} x = | -B^T A^{-1}  I      0   | | p | = | p - B^T A^{-1} u | = | w_p |
                \      0       0   G^{-1} / \ U /   \     G^{-1} U     /   \ w_U /
  */
  ierr = VecScatterBegin(xx, s->u, INSERT_VALUES, SCATTER_FORWARD, s->uScatter);                          CHKERRQ(ierr);
  ierr = VecScatterEnd(xx, s->u, INSERT_VALUES, SCATTER_FORWARD, s->uScatter);                            CHKERRQ(ierr);
  /* Permute u */
  if (s->colPerm != PETSC_NULL) {
    ierr = VecPermute(s->u, s->colPerm, PETSC_FALSE);                                                     CHKERRQ(ierr);
  }
  /* Velocity solve: Now w_u = s->uNew */
  ierr = SLESSolve(s->sles, s->u, s->uNew, &its);                                                         CHKERRQ(ierr);
  s->iter += its;
  PetscLogInfo(pc, "PCSchur: %d iterations in A solve %d\n", its, 0);
  /* Pressure solve: Now w_p = s->p */
  ierr = MatMultTranspose(s->B, s->uNew, s->p);                                                           CHKERRQ(ierr);
  ierr = VecScale(&minusOne, s->p);                                                                       CHKERRQ(ierr);
  ierr = VecScatterBegin(xx, s->p, ADD_VALUES, SCATTER_FORWARD, s->pScatter);                             CHKERRQ(ierr);
  ierr = VecScatterEnd(xx, s->p, ADD_VALUES, SCATTER_FORWARD, s->pScatter);                               CHKERRQ(ierr);

  /*            /   I    A^{-1} B (B^T A^{-1} B)^{-1}   0   \ / w_u \   / w_u + A^{-1} B (B^T A^{-1} B)^{-1} w_p \
     U^{-1} w = |   0       -(B^T A^{-1} B)^{-1}        0   | | w_p | = |         -(B^T A^{-1} B)^{-1} w_p       |
                \   0                 0                 I   / \ w_U /   \                      w_U               /
  */
  /* Schur solve: Now s->pNew = -S^{-1} w_p */
  ierr = SLESSolve(s->schurSles, s->p, s->pNew, &schurIts);                                               CHKERRQ(ierr);
  ierr = VecScale(&minusOne, s->pNew);                                                                    CHKERRQ(ierr);
  s->schurIter += schurIts;
  PetscLogInfo(pc, "PCSchur: %d iterations in B^T A^{-1} B solve\n", schurIts);
  /* Velocity solve: Now s->uNew2 = -A^{-1} B S^{-1} w_p */
  ierr = MatMult(s->B, s->pNew, s->u);                                                                    CHKERRQ(ierr);
  ierr = SLESSolve(s->sles, s->u, s->uNew2, &its);                                                        CHKERRQ(ierr);
  s->iter += its;
  PetscLogInfo(pc, "PCSchur: %d iterations in A solve %d\n", its, schurIts+1);
  /* Now s->uNew = w_u - A^{-1} B S^{-1} w_p */
  ierr = VecAXPY(&minusOne, s->uNew2, s->uNew);                                                           CHKERRQ(ierr);
  /* Permute uNew */
  if (s->rowPerm != PETSC_NULL) {
    ierr = VecPermute(s->uNew, s->rowPerm, PETSC_TRUE);                                                   CHKERRQ(ierr);
  }
  /* Put u and p in y */
  ierr = VecScatterBegin(s->uNew, yy, INSERT_VALUES, SCATTER_REVERSE, s->uScatter);                       CHKERRQ(ierr);
  ierr = VecScatterEnd(s->uNew, yy, INSERT_VALUES, SCATTER_REVERSE, s->uScatter);                         CHKERRQ(ierr);
  ierr = VecScatterBegin(s->pNew, yy, INSERT_VALUES, SCATTER_REVERSE, s->pScatter);                       CHKERRQ(ierr);
  ierr = VecScatterEnd(s->pNew, yy, INSERT_VALUES, SCATTER_REVERSE, s->pScatter);                         CHKERRQ(ierr);

  /* Project back: Apply P^+ = (P^T P)^{-1} P^T */
  if ((s->isConstrained == PETSC_TRUE) && (s->explicitConstraints == PETSC_FALSE)) {
    /* Project solution back onto constrained unknowns */
    ierr = MatMultTranspose(P, yy, z);                                                                    CHKERRQ(ierr);
    /* Particle solve: Add G^{-1} U = w_U to z */
    ierr = GVecSolveJacobianConstrained(z, x);                                                            CHKERRQ(ierr);
    /* Projection Solve: Now y = (P^T P)^{-1} z */
    ierr = MatMult(invP, z, y);                                                                           CHKERRQ(ierr);
  }

  ierr = PetscOptionsHasName(PETSC_NULL, "-pc_schur_test", &opt);                                         CHKERRQ(ierr);
  if (opt == PETSC_TRUE) {
    PetscReal solNorm;

    /* Check that we get A^{-1} A x = x */
    ierr = VecAXPY(&minusOne, y, tempSol);                                                                CHKERRQ(ierr);
    ierr = VecNorm(tempSol, NORM_2, &solNorm);                                                            CHKERRQ(ierr);
    ierr = VecDestroy(tempSol);                                                                           CHKERRQ(ierr);
    PetscPrintf(pc->comm, "Error in test solution: %g\n", solNorm);
  }
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCApplyTrans_Schur"
static int PCApplyTrans_Schur(PC pc, Vec x, Vec y)
{
  PetscFunctionBegin;
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCApplySymmetricLeft_Schur"
static int PCApplySymmetricLeft_Schur(PC pc, Vec x, Vec y)
{
  PetscFunctionBegin;
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCApplySymmetricRight_Schur"
static int PCApplySymmetricRight_Schur(PC pc, Vec x, Vec y)
{
  PetscFunctionBegin;
  PetscFunctionReturn(0);
}

#undef __FUNCT__  
#define __FUNCT__ "PCDestroy_Schur"
static int PCDestroy_Schur(PC pc)
{
  PC_Schur *s = (PC_Schur *) pc->data;
  int       ierr;

  PetscFunctionBegin;
  /* Destroy inner solvers */
  ierr = SLESDestroy(s->sles);                                                                            CHKERRQ(ierr);
  ierr = SLESDestroy(s->schurSles);                                                                       CHKERRQ(ierr);
  if (pc->setupcalled) {
    ierr = PCDestroyStructures_Schur(pc);                                                                 CHKERRQ(ierr);
  }
  if (s->mathViewer != PETSC_NULL) {
    ierr = PetscViewerDestroy(s->mathViewer);                                                             CHKERRQ(ierr);
  }
  PetscFree(s);
  PetscFunctionReturn(0);
}

EXTERN_C_BEGIN
#undef __FUNCT__  
#define __FUNCT__ "PCCreate_Schur"
int PCCreate_Schur(PC pc)
{
  PC_Schur *s;
  int       ierr;

  PetscFunctionBegin;
  ierr = PetscNew(PC_Schur, &s);                                                                          CHKERRQ(ierr);
  PetscLogObjectMemory(pc, sizeof(PC_Schur));
  pc->data                = (void *) s;

  pc->ops->setup               = PCSetUp_Schur;
  pc->ops->apply               = PCApply_Schur;
  pc->ops->applyrichardson     = PETSC_NULL;
  pc->ops->applyBA             = PETSC_NULL;
  pc->ops->applytranspose      = PCApplyTrans_Schur;
  pc->ops->applyBAtranspose    = PETSC_NULL;
  pc->ops->setfromoptions      = PCSetFromOptions_Schur;
  pc->ops->presolve            = PCPreSolve_Schur;
  pc->ops->postsolve           = PCPostSolve_Schur;
  pc->ops->getfactoredmatrix   = PETSC_NULL;
  pc->ops->applysymmetricleft  = PCApplySymmetricLeft_Schur;
  pc->ops->applysymmetricright = PCApplySymmetricRight_Schur;
  pc->ops->setuponblocks       = PETSC_NULL;
  pc->ops->destroy             = PCDestroy_Schur;
  pc->ops->view                = PETSC_NULL;
  pc->modifysubmatrices        = PETSC_NULL;

  s->grid                 = PETSC_NULL;
  s->isConstrained        = PETSC_FALSE;
  s->explicitConstraints  = PETSC_FALSE;
  s->useMath              = PETSC_FALSE;
  s->mathViewer           = PETSC_NULL;

  s->u                    = PETSC_NULL;
  s->uNew                 = PETSC_NULL;
  s->uNew2                = PETSC_NULL;
  s->p                    = PETSC_NULL;
  s->pNew                 = PETSC_NULL;
  s->uScatter             = PETSC_NULL;
  s->pScatter             = PETSC_NULL;
  s->projX                = PETSC_NULL;
  s->projY                = PETSC_NULL;
  s->projZ                = PETSC_NULL;

  s->sles                 = PETSC_NULL;
  s->schurSles            = PETSC_NULL;
  s->iter                 = 0;
  s->schurIter            = 0;

  s->numMomOps            = -1;
  s->gradOp               = GRADIENT;
  s->divOp                = DIVERGENCE;
  s->gradOpAlpha          = 0.0;
  s->momOps               = PETSC_NULL;
  s->momOpIsALE           = PETSC_NULL;
  s->momOpAlphas          = PETSC_NULL;
  s->sField               = -1;
  s->tField               = -1;
  s->sOrder               = PETSC_NULL;
  s->sLocOrder            = PETSC_NULL;
  s->tOrder               = PETSC_NULL;
  s->tLocOrder            = PETSC_NULL;
  s->A                    = PETSC_NULL;
  s->sparseA              = PETSC_NULL;
  s->B                    = PETSC_NULL;
  s->S                    = PETSC_NULL;
  s->lap                  = PETSC_NULL;
  s->rowPerm              = PETSC_NULL;
  s->colPerm              = PETSC_NULL;

  /* Create the inner solvers */
  ierr = PCSchurCreateInnerSLES(pc);                                                                      CHKERRQ(ierr);

  PetscFunctionReturn(0);
}
EXTERN_C_END