diff --git a/firedrake/assemble.py b/firedrake/assemble.py
index 4a7f5d76b6..9cf853279c 100644
--- a/firedrake/assemble.py
+++ b/firedrake/assemble.py
@@ -693,8 +693,6 @@ def _make_parloops(expr, tensor, bcs, diagonal, fc_params, assembly_rank):
     domains = expr.ufl_domains()
 
     if isinstance(expr, slate.TensorBase):
-        if diagonal:
-            raise NotImplementedError("Diagonal + slate not supported")
         kernels = slac.compile_expression(expr, compiler_parameters=form_compiler_parameters)
     else:
         kernels = tsfc_interface.compile_form(expr, "form", parameters=form_compiler_parameters, diagonal=diagonal)
diff --git a/firedrake/slate/slac/compiler.py b/firedrake/slate/slac/compiler.py
index 383a1d9ee3..45d6afd790 100644
--- a/firedrake/slate/slac/compiler.py
+++ b/firedrake/slate/slac/compiler.py
@@ -167,7 +167,7 @@ def generate_loopy_kernel(slate_expr, compiler_parameters=None):
 
     # Create a loopy builder for the Slate expression,
     # e.g. contains the loopy kernels coming from TSFC
-    gem_expr, var2terminal = slate_to_gem(slate_expr)
+    gem_expr, var2terminal = slate_to_gem(slate_expr, compiler_parameters["slate_compiler"])
 
     scalar_type = compiler_parameters["form_compiler"]["scalar_type"]
     slate_loopy, output_arg = gem_to_loopy(gem_expr, var2terminal, scalar_type)
diff --git a/firedrake/slate/slac/optimise.py b/firedrake/slate/slac/optimise.py
index 7f15b3eed1..4c6ce853b1 100644
--- a/firedrake/slate/slac/optimise.py
+++ b/firedrake/slate/slac/optimise.py
@@ -27,9 +27,12 @@ def optimise(expression, parameters):
 
     Returns: An optimised Slate expression
     """
-    # 1) Block optimisation
+    # 0) Block optimisation
     expression = push_block(expression)
 
+    # 1) DiagonalTensor optimisation
+    expression = push_diag(expression)
+
     # 2) Multiplication optimisation
     if expression.rank < 2:
         expression = push_mul(expression, parameters)
@@ -70,6 +73,8 @@ def _push_block_transpose(expr, self, indices):
 
 @_push_block.register(Add)
 @_push_block.register(Negative)
+@_push_block.register(DiagonalTensor)
+@_push_block.register(Reciprocal)
 def _push_block_distributive(expr, self, indices):
     """Distributes Blocks for these nodes"""
     return type(expr)(*map(self, expr.children, repeat(indices))) if indices else expr
@@ -111,6 +116,66 @@ def _push_block_block(expr, self, indices):
     return block
 
 
+def push_diag(expression):
+    """Executes a Slate compiler optimisation pass.
+    The optimisation is achieved by pushing DiagonalTensor from the outside to the inside of an expression.
+
+    :arg expression: A (potentially unoptimised) Slate expression.
+
+    Returns: An optimised Slate expression, where DiagonalTensors are sitting
+    on terminal tensors whereever possible.
+    """
+    mapper = MemoizerArg(_push_diag)
+    return mapper(expression, False)
+
+
+@singledispatch
+def _push_diag(expr, self, diag):
+    raise AssertionError("Cannot handle terminal type: %s" % type(expr))
+
+
+@_push_diag.register(Transpose)
+@_push_diag.register(Add)
+@_push_diag.register(Negative)
+def _push_diag_distributive(expr, self, diag):
+    """Distributes the DiagonalTensors into these nodes"""
+    return type(expr)(*map(self, expr.children, repeat(diag)))
+
+
+@_push_diag.register(Factorization)
+@_push_diag.register(Inverse)
+@_push_diag.register(Solve)
+@_push_diag.register(Mul)
+@_push_diag.register(Tensor)
+def _push_diag_stop(expr, self, diag):
+    """Diagonal Tensors cannot be pushed further into this set of nodes."""
+    expr = type(expr)(*map(self, expr.children, repeat(False))) if not expr.terminal else expr
+    return DiagonalTensor(expr) if diag else expr
+
+
+@_push_diag.register(Block)
+def _push_diag_block(expr, self, diag):
+    """Diagonal Tensors cannot be pushed further into this set of nodes."""
+    expr = type(expr)(*map(self, expr.children, repeat(False)), expr._indices) if not expr.terminal else expr
+    return DiagonalTensor(expr) if diag else expr
+
+
+@_push_diag.register(AssembledVector)
+@_push_diag.register(Reciprocal)
+def _push_diag_vectors(expr, self, diag):
+    """DiagonalTensors should not be pushed onto rank-1 tensors."""
+    if diag:
+        raise AssertionError("It is not legal to define DiagonalTensors on rank-1 tensors.")
+    else:
+        return expr
+
+
+@_push_diag.register(DiagonalTensor)
+def _push_diag_diag(expr, self, diag):
+    """DiagonalTensors are either pushed down or ignored when wrapped into another DiagonalTensor."""
+    return self(*expr.children, not diag)
+
+
 def push_mul(tensor, options):
     """Executes a Slate compiler optimisation pass.
     The optimisation is achieved by pushing coefficients from
@@ -179,6 +244,8 @@ def _drop_double_transpose_transpose(expr, self):
 @_drop_double_transpose.register(Mul)
 @_drop_double_transpose.register(Solve)
 @_drop_double_transpose.register(Inverse)
+@_drop_double_transpose.register(DiagonalTensor)
+@_drop_double_transpose.register(Reciprocal)
 def _drop_double_transpose_distributive(expr, self):
     """Distribute into the children of the expression. """
     return type(expr)(*map(self, expr.children))
@@ -202,6 +269,8 @@ def _push_mul_tensor(expr, self, state):
 
 
 @_push_mul.register(AssembledVector)
+@_push_mul.register(DiagonalTensor)
+@_push_mul.register(Reciprocal)
 def _push_mul_vector(expr, self, state):
     """Do not push into AssembledVectors."""
     return expr
@@ -220,8 +289,12 @@ def _push_mul_inverse(expr, self, state):
     with a coefficient into a Solve via A.inv*b = A.solve(b)
     or b*A^{-1}= (A.T.inv*b.T).T = A.T.solve(b.T).T ."""
     child, = expr.children
-    return (Solve(child, state.coeff) if state.pick_op
-            else Transpose(Solve(Transpose(child), Transpose(state.coeff))))
+    if expr.diagonal:
+        # Don't optimise further so that the translation to gem at a later can just spill ]1/a_ii[
+        return expr * state.coeff if state.pick_op else state.coeff * expr
+    else:
+        return (Solve(child, state.coeff) if state.pick_op
+                else Transpose(Solve(Transpose(child), Transpose(state.coeff))))
 
 
 @_push_mul.register(Transpose)
diff --git a/firedrake/slate/slac/tsfc_driver.py b/firedrake/slate/slac/tsfc_driver.py
index e94b8195c8..59272bcc23 100644
--- a/firedrake/slate/slac/tsfc_driver.py
+++ b/firedrake/slate/slac/tsfc_driver.py
@@ -60,7 +60,8 @@ def compile_terminal_form(tensor, prefix, *, tsfc_parameters=None, coffee=True):
         kernels = tsfc_compile(form,
                                subkernel_prefix,
                                parameters=tsfc_parameters,
-                               coffee=coffee, split=False)
+                               coffee=coffee, split=False, diagonal=tensor.diagonal)
+
         if kernels:
             cxt_k = ContextKernel(tensor=tensor,
                                   coefficients=form.coefficients(),
diff --git a/firedrake/slate/slac/utils.py b/firedrake/slate/slac/utils.py
index 21043bffac..f69ff64e90 100644
--- a/firedrake/slate/slac/utils.py
+++ b/firedrake/slate/slac/utils.py
@@ -6,7 +6,7 @@
 from ufl.algorithms.multifunction import MultiFunction
 
 from gem import (Literal, Sum, Product, Indexed, ComponentTensor, IndexSum,
-                 Solve, Inverse, Variable, view)
+                 Solve, Inverse, Variable, view, Delta, Index, Division)
 from gem import indices as make_indices
 from gem.node import Memoizer
 from gem.node import pre_traversal as traverse_dags
@@ -148,7 +148,7 @@ def visit_Symbol(self, o, *args, **kwargs):
         return SymbolWithFuncallIndexing(o.symbol, o.rank, o.offset)
 
 
-def slate_to_gem(expression):
+def slate_to_gem(expression, options):
     """Convert a slate expression to gem.
 
     :arg expression: A slate expression.
@@ -156,7 +156,7 @@ def slate_to_gem(expression):
         gem variables to UFL "terminal" forms.
     """
 
-    mapper, var2terminal = slate2gem(expression)
+    mapper, var2terminal = slate2gem(expression, options)
     return mapper, var2terminal
 
 
@@ -186,9 +186,37 @@ def _slate2gem_block(expr, self):
     return view(child, *(slice(idx, idx+extent) for idx, extent in zip(offsets, expr.shape)))
 
 
+@_slate2gem.register(sl.DiagonalTensor)
+def _slate2gem_diagonal(expr, self):
+    if not self.matfree:
+        A, = map(self, expr.children)
+        assert A.shape[0] == A.shape[1]
+        i, j = (Index(extent=s) for s in A.shape)
+        return ComponentTensor(Product(Indexed(A, (i, i)), Delta(i, j)), (i, j))
+    else:
+        raise NotImplementedError("Diagonals on Slate expressions are \
+                                   not implemented in a matrix-free manner yet.")
+
+
 @_slate2gem.register(sl.Inverse)
 def _slate2gem_inverse(expr, self):
-    return Inverse(*map(self, expr.children))
+    tensor, = expr.children
+    if expr.diagonal:
+        # optimise inverse on diagonal tensor by translating to
+        # matrix which contains the reciprocal values of the diagonal tensor
+        A, = map(self, expr.children)
+        i, j = (Index(extent=s) for s in A.shape)
+        return ComponentTensor(Product(Division(Literal(1), Indexed(A, (i, i))),
+                                       Delta(i, j)), (i, j))
+    else:
+        return Inverse(self(tensor))
+
+
+@_slate2gem.register(sl.Reciprocal)
+def _slate2gem_reciprocal(expr, self):
+    child, = map(self, expr.children)
+    indices = tuple(make_indices(len(child.shape)))
+    return ComponentTensor(Division(Literal(1.), Indexed(child, indices)), indices)
 
 
 @_slate2gem.register(sl.Solve)
@@ -237,9 +265,10 @@ def _slate2gem_factorization(expr, self):
     return A
 
 
-def slate2gem(expression):
+def slate2gem(expression, options):
     mapper = Memoizer(_slate2gem)
     mapper.var2terminal = OrderedDict()
+    mapper.matfree = options["replace_mul"]
     return mapper(expression), mapper.var2terminal
 
 
diff --git a/firedrake/slate/slate.py b/firedrake/slate/slate.py
index b57281f682..3fd9b85ec9 100644
--- a/firedrake/slate/slate.py
+++ b/firedrake/slate/slate.py
@@ -39,7 +39,8 @@
 
 __all__ = ['AssembledVector', 'Block', 'Factorization', 'Tensor',
            'Inverse', 'Transpose', 'Negative',
-           'Add', 'Mul', 'Solve', 'BlockAssembledVector']
+           'Add', 'Mul', 'Solve', 'BlockAssembledVector', 'DiagonalTensor',
+           'Reciprocal']
 
 
 class RemoveNegativeRestrictions(MultiFunction):
@@ -119,6 +120,7 @@ class TensorBase(object, metaclass=ABCMeta):
 
     terminal = False
     assembled = False
+    diagonal = False
 
     _id = count()
 
@@ -155,7 +157,7 @@ def expression_hash(self):
             elif isinstance(op, Factorization):
                 data = (type(op).__name__, op.decomposition, )
             elif isinstance(op, Tensor):
-                data = (op.form.signature(), )
+                data = (op.form.signature(), op.diagonal, )
             elif isinstance(op, (UnaryOp, BinaryOp)):
                 data = (type(op).__name__, )
             else:
@@ -815,14 +817,17 @@ class Tensor(TensorBase):
     operands = ()
     terminal = True
 
-    def __init__(self, form):
+    def __init__(self, form, diagonal=False):
         """Constructor for the Tensor class."""
         if not isinstance(form, Form):
             if isinstance(form, Function):
                 raise TypeError("Use AssembledVector instead of Tensor.")
             raise TypeError("Only UFL forms are acceptable inputs.")
 
-        r = len(form.arguments())
+        if self.diagonal:
+            assert len(form.arguments()) > 1, "Diagonal option only makes sense on rank-2 tensors."
+
+        r = len(form.arguments()) - diagonal
         if r not in (0, 1, 2):
             raise NotImplementedError("No support for tensors of rank %d." % r)
 
@@ -832,6 +837,7 @@ def __init__(self, form):
         super(Tensor, self).__init__()
 
         self.form = form
+        self.diagonal = diagonal
 
     @cached_property
     def arg_function_spaces(self):
@@ -842,7 +848,8 @@ def arg_function_spaces(self):
 
     def arguments(self):
         """Returns a tuple of arguments associated with the tensor."""
-        return self.form.arguments()
+        r = len(self.form.arguments()) - self.diagonal
+        return self.form.arguments()[0:r]
 
     def coefficients(self):
         """Returns a tuple of coefficients associated with the tensor."""
@@ -871,7 +878,7 @@ def __repr__(self):
     @cached_property
     def _key(self):
         """Returns a key for hash and equality."""
-        return (type(self), self.form)
+        return (type(self), self.form, self.diagonal)
 
 
 class TensorOp(TensorBase):
@@ -941,6 +948,37 @@ def __repr__(self):
         return "%s(%r)" % (type(self).__name__, tensor)
 
 
+class Reciprocal(UnaryOp):
+    """An abstract Slate class representing the reciprocal of a vector.
+    """
+
+    def __init__(self, A):
+        """Constructor for the Inverse class."""
+        assert A.rank == 1, "The tensor must be rank 1."
+
+        super(Reciprocal, self).__init__(A)
+
+    @cached_property
+    def arg_function_spaces(self):
+        """Returns a tuple of function spaces that the tensor
+        is defined on.
+        """
+        tensor, = self.operands
+        return tensor.arg_function_spaces
+
+    def arguments(self):
+        """Returns the expected arguments of the resulting tensor of
+        performing a specific unary operation on a tensor.
+        """
+        tensor, = self.operands
+        return tensor.arguments()
+
+    def _output_string(self, prec=None):
+        """Creates a string representation of the inverse of a tensor."""
+        tensor, = self.operands
+        return "(%s).reciprocal" % tensor
+
+
 class Inverse(UnaryOp):
     """An abstract Slate class representing the inverse of a tensor.
 
@@ -955,8 +993,9 @@ def __init__(self, A):
         assert A.shape[0] == A.shape[1], (
             "The inverse can only be computed on square tensors."
         )
+        self.diagonal = A.diagonal
 
-        if A.shape > (4, 4) and not isinstance(A, Factorization):
+        if A.shape > (4, 4) and not isinstance(A, Factorization) and not self.diagonal:
             A = Factorization(A, decomposition="PartialPivLU")
 
         super(Inverse, self).__init__(A)
@@ -1199,7 +1238,7 @@ def __init__(self, A, B, decomposition=None):
         decomposition = decomposition or "PartialPivLU"
 
         # Create a matrix factorization
-        A_factored = Factorization(A, decomposition=decomposition)
+        A_factored = Factorization(A, decomposition=decomposition) if not A.diagonal else A
 
         super(Solve, self).__init__(A_factored, B)
 
@@ -1220,6 +1259,43 @@ def arguments(self):
         return self._args
 
 
+class DiagonalTensor(UnaryOp):
+    """An abstract Slate class representing the diagonal of a tensor.
+
+    .. warning::
+
+       This class will raise an error if the tensor is not square.
+    """
+    diagonal = True
+
+    def __init__(self, A):
+        """Constructor for the Diagonal class."""
+        assert A.rank == 2, "The tensor must be rank 2."
+        assert A.shape[0] == A.shape[1], (
+            "The diagonal can only be computed on square tensors."
+        )
+
+        super(DiagonalTensor, self).__init__(A)
+
+    @cached_property
+    def arg_function_spaces(self):
+        """Returns a tuple of function spaces that the tensor
+        is defined on.
+        """
+        tensor, = self.operands
+        return tuple(arg.function_space() for arg in tensor.arguments())
+
+    def arguments(self):
+        """Returns a tuple of arguments associated with the tensor."""
+        tensor, = self.operands
+        return tensor.arguments()
+
+    def _output_string(self, prec=None):
+        """Creates a string representation of the diagonal of a tensor."""
+        tensor, = self.operands
+        return "(%s).diag" % tensor
+
+
 def space_equivalence(A, B):
     """Checks that two function spaces are equivalent.
 
@@ -1235,7 +1311,7 @@ def space_equivalence(A, B):
 
 # Establishes levels of precedence for Slate tensors
 precedences = [
-    [AssembledVector, Block, Factorization, Tensor],
+    [AssembledVector, Block, Factorization, Tensor, DiagonalTensor, Reciprocal],
     [Add],
     [Mul],
     [Solve],
diff --git a/firedrake/slate/static_condensation/hybridization.py b/firedrake/slate/static_condensation/hybridization.py
index 80037abaae..10c372e9dc 100644
--- a/firedrake/slate/static_condensation/hybridization.py
+++ b/firedrake/slate/static_condensation/hybridization.py
@@ -9,11 +9,13 @@
 from firedrake.matrix_free.operators import ImplicitMatrixContext
 from firedrake.petsc import PETSc
 from firedrake.parloops import par_loop, READ, INC
-from firedrake.slate.slate import Tensor, AssembledVector
+from firedrake.slate.slate import DiagonalTensor, Tensor, AssembledVector
 from pyop2.utils import as_tuple
+from firedrake.formmanipulation import split_form
 
+from firedrake.parameters import parameters
 
-__all__ = ['HybridizationPC']
+__all__ = ['HybridizationPC', 'SchurComplementBuilder']
 
 
 class HybridizationPC(SCBase):
@@ -41,7 +43,6 @@ def initialize(self, pc):
                                TrialFunction, TrialFunctions, TestFunction,
                                DirichletBC, assemble)
         from firedrake.assemble import allocate_matrix
-        from firedrake.formmanipulation import split_form
         from ufl.algorithms.replace import replace
 
         # Extract the problem context
@@ -199,24 +200,9 @@ def initialize(self, pc):
         # Make a SLATE tensor from Kform
         K = Tensor(Kform)
 
-        # Split forms
-        split_mixed_op = dict(split_form(Atilde.form))
-        id0, id1 = (self.vidx, self.pidx)
-        A00 = Tensor(split_mixed_op[(id0, id0)])
-        A01 = Tensor(split_mixed_op[(id0, id1)])
-        A10 = Tensor(split_mixed_op[(id1, id0)])
-        A11 = Tensor(split_mixed_op[(id1, id1)])
-        list_split_mixed_ops = [A00, A01, A10, A11]
-
-        split_trace_op = dict(split_form(K.form))
-        K0 = Tensor(split_trace_op[(0, id0)])
-        K1 = Tensor(split_trace_op[(0, id1)])
-        list_split_trace_ops = [K0, K1]
-
         # Build schur complement operator and right hand side
-        nested = PETSc.Options(prefix).getBool("nested_schur", False)
-        schur_rhs, schur_comp = self.build_schur(Atilde, K, list_split_mixed_ops,
-                                                 list_split_trace_ops, nested=nested)
+        self.schur_builder = SchurComplementBuilder(prefix, Atilde, K, pc, self.vidx, self.pidx)
+        schur_rhs, schur_comp = self.schur_builder.build_schur(self.broken_residual)
 
         # Assemble the Schur complement operator and right-hand side
         self.schur_rhs = Function(TraceSpace)
@@ -285,16 +271,11 @@ def initialize(self, pc):
             trace_ksp.setFromOptions()
 
         # Generate reconstruction calls
-        self._reconstruction_calls(list_split_mixed_ops, list_split_trace_ops)
+        self._reconstruction_calls()
 
-    def _reconstruction_calls(self, list_split_mixed_ops, list_split_trace_ops):
+    def _reconstruction_calls(self):
         """This generates the reconstruction calls for the unknowns using the
         Lagrange multipliers.
-
-        :arg split_mixed_op: a ``dict`` of split forms that make up the broken
-                             mixed operator from the original problem.
-        :arg split_trace_op: a ``dict`` of split forms that make up the trace
-                             contribution in the hybridized mixed system.
         """
         from firedrake import assemble
 
@@ -304,8 +285,12 @@ def _reconstruction_calls(self, list_split_mixed_ops, list_split_trace_ops):
         # TODO: When PyOP2 is able to write into mixed dats,
         # the reconstruction expressions can simplify into
         # one clean expression.
-        A, B, C, D = list_split_mixed_ops
-        K_0, K_1 = list_split_trace_ops
+
+        # reuse work from trace operator build
+        A, B, C, _ = self.schur_builder.list_split_mixed_ops
+        K_0, K_1 = self.schur_builder.list_split_trace_ops
+        Ahat = self.schur_builder.A00_inv_hat
+        S = self.schur_builder.inner_S
 
         # Split functions and reconstruct each bit separately
         split_residual = self.broken_residual.split()
@@ -316,10 +301,17 @@ def _reconstruction_calls(self, list_split_mixed_ops, list_split_trace_ops):
         u = split_sol[id1]
         lambdar = AssembledVector(self.trace_solution)
 
-        M = D - C * A.inv * B
-        R = K_1.T - C * A.inv * K_0.T
-        u_rec = M.solve(f - C * A.inv * g - R * lambdar,
-                        decomposition="PartialPivLU")
+        R = K_1.T - C * Ahat * K_0.T
+        rhs = f - C * Ahat * g - R * lambdar
+        if self.schur_builder.schur_approx or self.schur_builder.jacobi_S:
+            Shat = self.schur_builder.inner_S_approx_inv_hat
+            if self.schur_builder.preonly_S:
+                S = Shat
+            else:
+                S = Shat * S
+                rhs = Shat * rhs
+
+        u_rec = S.solve(rhs, decomposition="PartialPivLU")
         self._sub_unknown = functools.partial(assemble,
                                               u_rec,
                                               tensor=u,
@@ -334,49 +326,6 @@ def _reconstruction_calls(self, list_split_mixed_ops, list_split_trace_ops):
                                                form_compiler_parameters=self.ctx.fc_params,
                                                assembly_type="residual")
 
-    def build_schur(self, Atilde, K, list_split_mixed_ops, list_split_trace_ops, nested=False):
-        """The Schur complement in the operators of the trace solve contains
-        the inverse on a mixed system.  Users may want this inverse to be treated
-        with another schur complement.
-
-        Let the mixed matrix Atilde be called A here,
-        then the "nested" options rewrites with a Schur decomposition
-        as the following.
-
-        .. code-block:: text
-
-                A.inv = [[I, -A00.inv * A01]    *   [[A00.inv, 0    ]   *   [[I,             0]
-                        [0,  I             ]]       [0,        S.inv]]      [-A10* A00.inv,  I]]
-                        --------------------        -----------------       -------------------
-                        block1                      block2                  block3
-
-        with the (inner) schur complement S = A11 - A10 * A00.inv * A01
-        """
-
-        if nested:
-            A00, A01, A10, A11 = list_split_mixed_ops
-            K0, K1 = list_split_trace_ops
-            broken_residual = self.broken_residual.split()
-            split_broken_res = [AssembledVector(broken_residual[self.vidx]),
-                                AssembledVector(broken_residual[self.pidx])]
-
-            # inner schur complement
-            S = (A11 - A10 * A00.inv * A01)
-            # K * block1
-            K_Ainv_block1 = [K0, -K0 * A00.inv * A01 + K1]
-            # K * block1 * block2
-            K_Ainv_block2 = [K_Ainv_block1[0] * A00.inv, K_Ainv_block1[1] * S.inv]
-            # K * block1 * block2 * block3
-            K_Ainv_block3 = [K_Ainv_block2[0] - K_Ainv_block2[1] * A10 * A00.inv, K_Ainv_block2[1]]
-            # K * block1 * block2 * block3 * broken residual
-            schur_rhs = (K_Ainv_block3[0] * split_broken_res[0] + K_Ainv_block3[1] * split_broken_res[1])
-            # K * block1 * block2 * block3 * K.T
-            schur_comp = K_Ainv_block3[0] * K0.T + K_Ainv_block3[1] * K1.T
-        else:
-            schur_rhs = K * Atilde.inv * AssembledVector(self.broken_residual)
-            schur_comp = K * Atilde.inv * K.T
-        return schur_rhs, schur_comp
-
     @PETSc.Log.EventDecorator("HybridUpdate")
     def update(self, pc):
         """Update by assembling into the operator. No need to
@@ -452,7 +401,8 @@ def backward_substitution(self, pc, y):
 
         # We assemble the unknown which is an expression
         # of the first eliminated variable.
-        self._sub_unknown()
+        with PETSc.Log.Event("RecoverFirstElim"):
+            self._sub_unknown()
         # Recover the eliminated unknown
         self._elim_unknown()
 
@@ -486,3 +436,267 @@ def view(self, pc, viewer=None):
             self.trace_ksp.view(viewer)
             viewer.printfASCII("Locally reconstructing solutions.\n")
             viewer.printfASCII("Projecting broken flux into HDiv space.\n")
+
+    def getSchurComplementBuilder(self):
+        return self.schur_builder
+
+
+class SchurComplementBuilder(object):
+
+    """A Slate-based Schur complement expression builder. The expression is
+    used in the trace system solve and parts of it in the reconstruction
+    calls of the other two variables of the hybridised system.
+    How the Schur complement if constructed, and in particular how the local inverse of the
+    mixed matrix is built, is controlled with PETSc options. All corresponding PETSc options
+    start with ``hybridization_localsolve``.
+    The following option sets are valid together with the usual set of hybridisation options:
+
+    .. code-block:: text
+
+        {'localsolve': {'ksp_type': 'preonly',
+                        'pc_type': 'fieldsplit',
+                        'pc_fieldsplit_type': 'schur'}}
+
+    A Schur complement is requested for the mixed matrix inverse which appears inside the
+    Schur complement of the trace system solve. The Schur complements are then nested.
+    For details see defition of :meth:`build_schur`. No fieldsplit options are set so all
+    local inverses are calculated explicitly.
+
+    .. code-block:: text
+
+        'localsolve': {'ksp_type': 'preonly',
+                       'pc_type': 'fieldsplit',
+                       'pc_fieldsplit_type': 'schur',
+                       'fieldsplit_1': {'ksp_type': 'default',
+                                        'pc_type': 'python',
+                                        'pc_python_type': __name__ + '.DGLaplacian'}}
+
+    The inverse of the Schur complement inside the Schur decomposition of the mixed matrix inverse
+    is approximated by a default solver (LU in the matrix-explicit case) which is preconditioned
+    by a user-defined operator, e.g. a DG Laplacian, see :meth:`build_inner_S_inv`.
+    So :math:`P_S * S * x = P_S * b`.
+
+    .. code-block:: text
+
+        'localsolve': {'ksp_type': 'preonly',
+                        'pc_type': 'fieldsplit',
+                        'pc_fieldsplit_type': 'schur',
+                        'fieldsplit_1': {'ksp_type': 'default',
+                                        'pc_type': 'python',
+                                        'pc_python_type': __name__ + '.DGLaplacian',
+                                        'aux_ksp_type': 'preonly'}
+                                        'aux_pc_type': 'jacobi'}}}}
+
+    The inverse of the Schur complement inside the Schur decomposition of the mixed matrix inverse
+    is approximated by a default solver (LU in the matrix-explicit case) which is preconditioned
+    by a user-defined operator, e.g. a DG Laplacian. The inverse of the preconditioning matrix is
+    approximated through the inverse of only the diagonal of the provided operator, see
+    :meth:`build_Sapprox_inv`. So :math:`diag(P_S).inv * S * x = diag(P_S).inv * b`.
+
+    .. code-block:: text
+
+        'localsolve': {'ksp_type': 'preonly',
+                       'pc_type': 'fieldsplit',
+                       'pc_fieldsplit_type': 'schur',
+                       'fieldsplit_0': {'ksp_type': 'default',
+                                        'pc_type': 'jacobi'}
+
+    The inverse of the :math:`A_{00}` block of the mixed matrix is approximated by a default solver
+    (LU in the matrix-explicit case) which is preconditioned by the diagonal matrix of :math:`A_{00},
+    see :meth:`build_A00_inv`. So :math:`diag(A_{00}).inv * A_{00} * x = diag(A_{00}).inv * b`.
+
+    .. code-block:: text
+
+        'localsolve': {'ksp_type': 'preonly',
+                       'pc_type': 'fieldsplit',
+                       'pc_fieldsplit_type': 'None',
+                       'fieldsplit_0':  ...
+                       'fieldsplit_1':  ...
+
+    All the options for ``fieldsplit_`` are still valid if ``'pc_fieldsplit_type': 'None'.`` In this case
+    the mixed matrix inverse which appears inside the Schur complement of the trace system solve
+    is calculated explicitly, but the local inverses of :math:`A_{00}` and the Schur complement
+    in the reconstructions calls are still treated according to the options in ``fieldsplit_``.
+
+    """
+
+    def __init__(self, prefix, Atilde, K, pc, vidx, pidx):
+        # set options, operators and order of sub-operators
+        self.Atilde = Atilde
+        self.K = K
+        self.vidx = vidx
+        self.pidx = pidx
+        self._split_mixed_operator()
+        self.prefix = prefix + "localsolve_"
+
+        # prefixes
+        self._retrieve_options(pc)
+
+        # build all required inverses
+        self.A00_inv_hat = self.build_A00_inv()
+        self.inner_S = self.build_inner_S()
+        self.inner_S_approx_inv_hat = self.build_Sapprox_inv()
+        self.inner_S_inv_hat = self.build_inner_S_inv()
+
+    def _split_mixed_operator(self):
+        split_mixed_op = dict(split_form(self.Atilde.form))
+        id0, id1 = (self.vidx, self.pidx)
+        A00 = Tensor(split_mixed_op[(id0, id0)])
+        A01 = Tensor(split_mixed_op[(id0, id1)])
+        A10 = Tensor(split_mixed_op[(id1, id0)])
+        A11 = Tensor(split_mixed_op[(id1, id1)])
+        self.list_split_mixed_ops = [A00, A01, A10, A11]
+
+        split_trace_op = dict(split_form(self.K.form))
+        K0 = Tensor(split_trace_op[(0, id0)])
+        K1 = Tensor(split_trace_op[(0, id1)])
+        self.list_split_trace_ops = [K0, K1]
+
+    def _check_options(self, valid):
+        default = object()
+        opts = PETSc.Options(self.prefix)
+        for key, supported in valid:
+            value = opts.getString(key, default=default)
+            if value is not default and value not in supported:
+                raise ValueError(f"Unsupported value ({value}) for '{self.prefix + key}'. "
+                                 f"Should be one of {supported}")
+
+    def _retrieve_options(self, pc):
+        get_option = lambda key: PETSc.Options(self.prefix).getString(key, default="")
+
+        # Get options for Schur complement decomposition
+        self._check_options([("ksp_type", {"preonly"}), ("pc_type", {"fieldsplit"}), ("pc_fieldsplit_type", {"schur"})])
+        self.nested = (get_option("ksp_type") == "preonly"
+                       and get_option("pc_type") == "fieldsplit"
+                       and get_option("pc_fieldsplit_type") == "schur")
+
+        # Get preconditioning options for A00
+        fs0, fs1 = ("fieldsplit_"+str(idx) for idx in (self.vidx, self.pidx))
+        self._check_options([(fs0+"ksp_type", {"preonly", "default"}), (fs0+"pc_type", {"jacobi"})])
+        self.preonly_A00 = get_option(fs0+"_ksp_type") == "preonly"
+        self.jacobi_A00 = get_option(fs0+"_pc_type") == "jacobi"
+
+        # Get preconditioning options for the Schur complement
+        self._check_options([(fs1+"ksp_type", {"preonly", "default"}), (fs1+"pc_type", {"jacobi", "python"})])
+        self.preonly_S = get_option(fs1+"_ksp_type") == "preonly"
+        self.jacobi_S = get_option(fs1+"_pc_type") == "jacobi"
+
+        # Get user supplied operator and its options
+        self.schur_approx = (self.retrieve_user_S_approx(pc, get_option(fs1+"_pc_python_type"))
+                             if get_option(fs1+"_pc_type") == "python"
+                             else None)
+        self._check_options([(fs1+"aux_ksp_type", {"preonly", "default"}), (fs1+"aux_pc_type", {"jacobi"})])
+        self.preonly_Shat = get_option(fs1+"_aux_ksp_type") == "preonly"
+        self.jacobi_Shat = get_option(fs1+"_aux_pc_type") == "jacobi"
+
+        if self.jacobi_Shat or self.jacobi_A00:
+            assert parameters["slate_compiler"]["optimise"], ("Local systems should only get preconditioned with "
+                                                              "a preconditioning matrix if the Slate optimiser replaces "
+                                                              "inverses by solves.")
+
+    def build_inner_S(self):
+        """Build the inner Schur complement."""
+        _, A01, A10, A11 = self.list_split_mixed_ops
+        return A11 - A10 * self.A00_inv_hat * A01
+
+    def inv(self, A, P, prec, preonly=False):
+        """ Calculates the inverse of an operator A.
+            The inverse is potentially approximated through a solve
+            which is potentially preconditioned with the preconditioner P
+            if prec is True.
+            The inverse of A may be just approximated with the inverse of P
+            if prec and replace.
+        """
+        return (P if prec and preonly else
+                (P*A).inv * P if prec else
+                A.inv)
+
+    def build_inner_S_inv(self):
+        """ Calculates the inverse of the schur complement.
+            The inverse is potentially approximated through a solve
+            which is potentially preconditioned with the preconditioner P.
+        """
+        A = self.inner_S
+        P = self.inner_S_approx_inv_hat
+        prec = bool(self.schur_approx) or self.jacobi_S
+        return self.inv(A, P, prec, self.preonly_S)
+
+    def build_Sapprox_inv(self):
+        """ Calculates the inverse of preconditioner to the Schur complement,
+            which can be either the schur complement approximation provided by the user
+            or jacobi.
+            The inverse is potentially approximated through a solve
+            which is potentially preconditioned with jacobi.
+        """
+        prec = (bool(self.schur_approx) and self.jacobi_Shat) or self.jacobi_S
+        A = self.schur_approx if self.schur_approx else self.inner_S
+        P = DiagonalTensor(A).inv
+        preonly = self.preonly_Shat if self.schur_approx else True
+        return self.inv(A, P, prec, preonly)
+
+    def build_A00_inv(self):
+        """ Calculates the inverse of :math:`A_{00}`, the (0,0)-block of the mixed matrix Atilde.
+            The inverse is potentially approximated through a solve
+            which is potentially preconditioned with jacobi.
+        """
+        A, _, _, _ = self.list_split_mixed_ops
+        P = DiagonalTensor(A).inv
+        return self.inv(A, P, self.jacobi_A00, self.preonly_A00)
+
+    def retrieve_user_S_approx(self, pc, usercode):
+        """Retrieve a user-defined :class:firedrake.preconditioners.AuxiliaryOperator from the PETSc Options,
+        which is an approximation to the Schur complement and its inverse is used
+        to precondition the local solve in the reconstruction calls (e.g.).
+        """
+        _, _, _, A11 = self.list_split_mixed_ops
+        test, trial = A11.arguments()
+        if usercode != "":
+            (modname, funname) = usercode.rsplit('.', 1)
+            mod = __import__(modname)
+            fun = getattr(mod, funname)
+            if isinstance(fun, type):
+                fun = fun()
+            return Tensor(fun.form(pc, test, trial)[0])
+        else:
+            return None
+
+    def build_schur(self, rhs):
+        """The Schur complement in the operators of the trace solve contains
+        the inverse on a mixed system.  Users may want this inverse to be treated
+        with another Schur complement.
+
+        Let the mixed matrix Atilde be called A here.
+        Then, if a nested schur complement is requested, the inverse of Atilde
+        is rewritten with help of a a Schur decomposition as follows.
+
+        .. code-block:: text
+
+                A.inv = [[I, -A00.inv * A01]    *   [[A00.inv, 0    ]   *   [[I,             0]
+                        [0,  I             ]]       [0,        S.inv]]      [-A10* A00.inv,  I]]
+                        --------------------        -----------------       -------------------
+                        block1                      block2                  block3
+                with the (inner) schur complement S = A11 - A10 * A00.inv * A01
+        """
+
+        if self.nested:
+            _, A01, A10, _ = self.list_split_mixed_ops
+            K0, K1 = self.list_split_trace_ops
+            broken_residual = rhs.split()
+            R = [AssembledVector(broken_residual[self.vidx]),
+                 AssembledVector(broken_residual[self.pidx])]
+            # K * block1
+            K_Ainv_block1 = [K0, -K0 * self.A00_inv_hat * A01 + K1]
+            # K * block1 * block2
+            K_Ainv_block2 = [K_Ainv_block1[0] * self.A00_inv_hat,
+                             K_Ainv_block1[1] * self.inner_S_inv_hat]
+            # K * block1 * block2 * block3
+            K_Ainv_block3 = [K_Ainv_block2[0] - K_Ainv_block2[1] * A10 * self.A00_inv_hat,
+                             K_Ainv_block2[1]]
+            # K * block1 * block2 * block3 * broken residual
+            schur_rhs = (K_Ainv_block3[0] * R[0] + K_Ainv_block3[1] * R[1])
+            # K * block1 * block2 * block3 * K.T
+            schur_comp = K_Ainv_block3[0] * K0.T + K_Ainv_block3[1] * K1.T
+        else:
+            schur_rhs = self.K * self.Atilde.inv * AssembledVector(rhs)
+            schur_comp = self.K * self.Atilde.inv * self.K.T
+        return schur_rhs, schur_comp
diff --git a/tests/slate/test_assemble_tensors.py b/tests/slate/test_assemble_tensors.py
index fd8c807151..2891d91251 100644
--- a/tests/slate/test_assemble_tensors.py
+++ b/tests/slate/test_assemble_tensors.py
@@ -81,6 +81,22 @@ def mass(function_space):
     return inner(u, v) * dx
 
 
+@pytest.fixture
+def matrix_mixed_facet():
+    mesh = UnitSquareMesh(2, 2)
+    U = FunctionSpace(mesh, "RT", 1)
+    V = FunctionSpace(mesh, "DG", 0)
+    T = FunctionSpace(mesh, "HDiv Trace", 0)
+    n = FacetNormal(mesh)
+    W = U * V * T
+    u, p, lambdar = TrialFunctions(W)
+    w, q, gammar = TestFunctions(W)
+
+    return (inner(u, w)*dx + p*q*dx - div(w)*p*dx + q*div(u)*dx
+            + lambdar('+')*jump(w, n=n)*dS + gammar('+')*jump(u, n=n)*dS
+            + lambdar*gammar*ds)
+
+
 @pytest.fixture
 def rank_one_tensor(mass, f):
     return Tensor(action(mass, f))
@@ -270,3 +286,56 @@ def test_matrix_subblocks(mesh):
     for tensor, form in items:
         ref = assemble(form).M.values
         assert np.allclose(assemble(tensor).M.values, ref, rtol=1e-14)
+
+
+def test_diagonal(mass, matrix_mixed_facet):
+    n, _ = Tensor(mass).shape
+
+    # test vector built from diagonal
+    res = assemble(Tensor(mass, diagonal=True)).dat.data
+    ref = assemble(mass, diagonal=True).dat.data
+    for r, d in zip(ref, res):
+        assert np.allclose(r, d, rtol=1e-14)
+
+    # test matrix built from diagonal
+    res = assemble(DiagonalTensor(Tensor(mass))).M.values
+    ref = assemble(mass, diagonal=True).dat.data
+    ref = np.concatenate(ref) if len(np.shape(ref)) > 1 else ref  # vectorspace
+    ref = np.concatenate(ref) if len(np.shape(ref)) > 1 else ref  # tensorspace
+    for r, d in zip(ref, np.diag(res)):
+        assert np.allclose(r, d, rtol=1e-14)
+
+    # test matrix built from diagonal for non mass matrix
+    res2 = assemble(DiagonalTensor(Tensor(matrix_mixed_facet))).M.values
+    ref2 = np.concatenate(assemble(matrix_mixed_facet, diagonal=True).dat.data)
+    for r, d in zip(res2, np.diag(ref2)):
+        assert np.allclose(r, d, rtol=1e-14)
+
+    # test matrix built from diagonal
+    # for a Slate expression on a non mass matrix
+    A = Tensor(matrix_mixed_facet)
+    res3 = assemble(DiagonalTensor(A+A)).M.values
+    ref3 = np.concatenate(assemble(matrix_mixed_facet+matrix_mixed_facet, diagonal=True).dat.data)
+    for r, d in zip(res3, np.diag(ref3)):
+        assert np.allclose(r, d, rtol=1e-14)
+
+
+@pytest.mark.parametrize("function_space", ["dg0"], indirect=True)
+def test_reciprocal(function_space):
+    # test reciprocal of vector built from diagonal
+    # note: reciprocal does not commute with addition so one can only test DG
+    u = TrialFunction(function_space)
+    v = TestFunction(function_space)
+    mass = inner(u, v) * dx
+    res = assemble(Reciprocal(Tensor(mass, diagonal=True))).dat.data
+    ref = assemble(mass, diagonal=True).dat.data
+    for r, d in zip([1./d for d in ref], res):
+        assert np.allclose(r, d, rtol=1e-14)
+
+    # test inverse of matrix built from diagonal
+    # for a Slate expression on a mass matrix
+    A = Tensor(mass)
+    res = assemble(DiagonalTensor(A+A).inv).M.values
+    ref = [1./d for d in assemble(mass + mass, diagonal=True).dat.data]
+    for r, d in zip(res, np.diag(ref)):
+        assert np.allclose(r, d, rtol=1e-14)
diff --git a/tests/slate/test_optimise.py b/tests/slate/test_optimise.py
index 7d48649740..b8dc5b2f01 100644
--- a/tests/slate/test_optimise.py
+++ b/tests/slate/test_optimise.py
@@ -329,6 +329,28 @@ def test_drop_transposes(TC_non_symm):
     compare_slate_tensors(expressions, opt_expressions)
 
 
+#######################################
+# Test diagonal optimisation pass
+#######################################
+def test_push_diagonal(TC_non_symm):
+    """Test Optimisers's ability to push DiagonalTensors inside expressions."""
+    A, C = TC_non_symm
+
+    expressions = [DiagonalTensor(A), DiagonalTensor(A+A),
+                   DiagonalTensor(-A), DiagonalTensor(A*A),
+                   DiagonalTensor(A).inv]
+    opt_expressions = [DiagonalTensor(A), DiagonalTensor(A)+DiagonalTensor(A),
+                       -DiagonalTensor(A), DiagonalTensor(A*A),
+                       DiagonalTensor(A).inv]
+    compare_tensor_expressions(expressions)
+    compare_slate_tensors(expressions, opt_expressions)
+
+    expressions = [DiagonalTensor(A+A).solve(C)]
+    opt_expressions = [(DiagonalTensor(A)+DiagonalTensor(A)).solve(C)]
+    compare_vector_expressions(expressions)
+    compare_slate_tensors(expressions, opt_expressions)
+
+
 #######################################
 # Helper functions
 #######################################
diff --git a/tests/slate/test_slate_hybridization.py b/tests/slate/test_slate_hybridization.py
index d11ce0fb00..6f14c8145d 100644
--- a/tests/slate/test_slate_hybridization.py
+++ b/tests/slate/test_slate_hybridization.py
@@ -27,7 +27,8 @@
 
 
 def setup_poisson():
-    mesh = UnitSquareMesh(1, 1)
+    n = 2
+    mesh = UnitSquareMesh(n, n)
     U = FunctionSpace(mesh, "RT", 4)
     V = FunctionSpace(mesh, "DG", 3)
     W = U * V
@@ -46,6 +47,36 @@ def setup_poisson():
     return a, L, W
 
 
+def setup_poisson_3D():
+    p = 3
+    n = 2
+    mesh = SquareMesh(n, n, 1, quadrilateral=True)
+    mesh = ExtrudedMesh(mesh, n)
+    RT = FiniteElement("RTCF", quadrilateral, p+1)
+    DG_v = FiniteElement("DG", interval, p)
+    DG_h = FiniteElement("DQ", quadrilateral, p)
+    CG = FiniteElement("CG", interval, p+1)
+    HDiv_ele = EnrichedElement(HDiv(TensorProductElement(RT, DG_v)),
+                               HDiv(TensorProductElement(DG_h, CG)))
+    V = FunctionSpace(mesh, HDiv_ele)
+    U = FunctionSpace(mesh, "DQ", p)
+    W = V * U
+    sigma, u = TrialFunctions(W)
+    tau, v = TestFunctions(W)
+    V, U = W.split()
+    f = Function(U)
+    x, y, z = SpatialCoordinate(mesh)
+    expr = (1+12*pi*pi)*cos(100*pi*x)*cos(100*pi*y)*cos(100*pi*z)
+    f.interpolate(expr)
+    a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(degree=8)
+    L = -f*v*dx(degree=8)
+    return a, L, W
+
+
+def options_check(builder, expected):
+    return all(bool(getattr(builder, k)) == v for k, v in expected.items())
+
+
 @pytest.mark.parametrize(("degree", "hdiv_family", "quadrilateral"),
                          [(1, "RT", False), (1, "RTCF", True),
                           (2, "RT", False), (2, "RTCF", True)])
@@ -101,6 +132,36 @@ def test_slate_hybridization(degree, hdiv_family, quadrilateral):
     assert u_err < 1e-11
 
 
+def test_slate_hybridization_wrong_option():
+    a, L, W = setup_poisson()
+
+    w = Function(W)
+    params = {'mat_type': 'matfree',
+              'ksp_type': 'preonly',
+              'pc_type': 'python',
+              'pc_python_type': 'firedrake.HybridizationPC',
+              'hybridization': {'ksp_type': 'preonly',
+                                'pc_type': 'lu',
+                                'localsolve': {'ksp_type': 'preonly',
+                                               'pc_type': 'fieldsplit',
+                                               'pc_fieldsplit_type': 'frog'}}}
+    eq = a == L
+    problem = LinearVariationalProblem(eq.lhs, eq.rhs, w)
+    solver = LinearVariationalSolver(problem, solver_parameters=params)
+    with pytest.raises(ValueError):
+        # HybridizationPC isn't called directly from the Python interpreter,
+        # it's a callback that PETSc calls. This means that the call stack from pytest
+        # down to HybridizationPC goes via PETSc C code, which interferes with the exception
+        # before it is observed outside. Hence removing PETSc's error handler
+        # makes the problem go away, because PETSc stops interfering.
+        # We need to repush the error handler because popErrorHandler globally changes
+        # the system state for all future tests.
+        from firedrake.petsc import PETSc
+        PETSc.Sys.pushErrorHandler("ignore")
+        solver.solve()
+        PETSc.Sys.popErrorHandler("ignore")
+
+
 def test_slate_hybridization_nested_schur():
     a, L, W = setup_poisson()
 
@@ -111,8 +172,20 @@ def test_slate_hybridization_nested_schur():
               'pc_python_type': 'firedrake.HybridizationPC',
               'hybridization': {'ksp_type': 'preonly',
                                 'pc_type': 'lu',
-                                'nested_schur': 'true'}}
-    solve(a == L, w, solver_parameters=params)
+                                'localsolve': {'ksp_type': 'preonly',
+                                               'pc_type': 'fieldsplit',
+                                               'pc_fieldsplit_type': 'schur'}}}
+
+    eq = a == L
+    problem = LinearVariationalProblem(eq.lhs, eq.rhs, w)
+    solver = LinearVariationalSolver(problem, solver_parameters=params)
+    solver.solve()
+    expected = {'nested': True,
+                'preonly_A00': False, 'jacobi_A00': False,
+                'schur_approx': False,
+                'preonly_Shat': False, 'jacobi_Shat': False}
+    builder = solver.snes.ksp.pc.getPythonContext().getSchurComplementBuilder()
+    assert options_check(builder, expected), "Some solver options have not ended up in the PC as wanted."
     sigma_h, u_h = w.split()
 
     w2 = Function(W)
@@ -130,3 +203,305 @@ def test_slate_hybridization_nested_schur():
 
     assert sigma_err < 1e-11
     assert u_err < 1e-11
+
+
+class DGLaplacian(AuxiliaryOperatorPC):
+    def form(self, pc, u, v):
+        W = u.function_space()
+        n = FacetNormal(W.mesh())
+        alpha = Constant(3**3)
+        gamma = Constant(4**3)
+        h = CellSize(W.mesh())
+        h_avg = (h('+') + h('-'))/2
+        a_dg = -(inner(grad(u), grad(v))*dx
+                 - inner(jump(u, n), avg(grad(v)))*dS
+                 - inner(avg(grad(u)), jump(v, n), )*dS
+                 + alpha/h_avg * inner(jump(u, n), jump(v, n))*dS
+                 - inner(u*n, grad(v))*ds
+                 - inner(grad(u), v*n)*ds
+                 + (gamma/h)*inner(u, v)*ds)
+        bcs = None
+        return (a_dg, bcs)
+
+
+class DGLaplacian3D(AuxiliaryOperatorPC):
+    def form(self, pc, u, v):
+        W = u.function_space()
+        n = FacetNormal(W.mesh())
+        gamma = Constant(4.**3)
+        h = CellVolume(W.mesh())/FacetArea(W.mesh())
+
+        a_dg = -(dot(grad(v), grad(u))*dx(degree=8)
+                 - dot(grad(v), (u)*n)*ds_v(degree=8)
+                 - dot(v*n, grad(u))*ds_v(degree=8)
+                 + gamma/h*dot(v, u)*ds_v(degree=8)
+                 - dot(grad(v), (u)*n)*ds_t(degree=8)
+                 - dot(v*n, grad(u))*ds_t(degree=8)
+                 + gamma/h*dot(v, u)*ds_t(degree=8)
+                 - dot(grad(v), (u)*n)*ds_b(degree=8)
+                 - dot(v*n, grad(u))*ds_b(degree=8)
+                 + gamma/h*dot(v, u)*ds_b(degree=8))
+
+        bcs = []
+        return (a_dg, bcs)
+
+
+def test_mixed_poisson_approximated_schur():
+    """A test, which compares a solution to a 2D mixed Poisson problem solved
+    globally matrixfree with a HybridizationPC and CG on the trace system to
+    a solution with uses a user supplied operator as preconditioner to the
+    Schur solver.
+
+    NOTE: With the setup in this test, using the approximated schur complemement
+    defined as DGLaplacian as a preconditioner to the schur complement,
+    reduces the condition number of the local solve from 16.77 to 6.06.
+    """
+    # setup FEM
+    a, L, W = setup_poisson()
+
+    # setup first solver
+    w = Function(W)
+    params = {'ksp_type': 'preonly',
+              'pc_type': 'python',
+              'mat_type': 'matfree',
+              'pc_python_type': 'firedrake.HybridizationPC',
+              'hybridization': {'ksp_type': 'cg',
+                                'pc_type': 'none',
+                                'ksp_rtol': 1e-8,
+                                'mat_type': 'matfree',
+                                'localsolve': {'ksp_type': 'preonly',
+                                               'pc_type': 'fieldsplit',
+                                               'pc_fieldsplit_type': 'schur',
+                                               'fieldsplit_1': {'ksp_type': 'default',
+                                                                'pc_type': 'python',
+                                                                'pc_python_type': __name__ + '.DGLaplacian'}}}}
+
+    eq = a == L
+    problem = LinearVariationalProblem(eq.lhs, eq.rhs, w)
+    solver = LinearVariationalSolver(problem, solver_parameters=params)
+    solver.solve()
+
+    # double-check options are set as expected
+    expected = {'nested': True,
+                'preonly_A00': False, 'jacobi_A00': False,
+                'schur_approx': True,
+                'preonly_Shat': False, 'jacobi_Shat': False}
+    builder = solver.snes.ksp.pc.getPythonContext().getSchurComplementBuilder()
+    assert options_check(builder, expected), "Some solver options have not ended up in the PC as wanted."
+
+    sigma_h, u_h = w.split()
+
+    # setup second solver
+    w2 = Function(W)
+    aij_params = {'ksp_type': 'preonly',
+                  'pc_type': 'python',
+                  'mat_type': 'matfree',
+                  'pc_python_type': 'firedrake.HybridizationPC',
+                  'hybridization': {'ksp_type': 'cg',
+                                    'pc_type': 'none',
+                                    'ksp_rtol': 1e-8,
+                                    'mat_type': 'matfree'}}
+    solve(a == L, w2, solver_parameters=aij_params)
+    _sigma, _u = w2.split()
+
+    # Return the L2 error
+    sigma_err = errornorm(sigma_h, _sigma)
+    u_err = errornorm(u_h, _u)
+
+    assert sigma_err < 1e-8
+    assert u_err < 1e-8
+
+
+def test_slate_hybridization_jacobi_prec_A00():
+    """A test, which compares a solution to a 3D mixed Poisson problem solved
+    globally matrixfree with a HybridizationPC and CG on the trace system to
+    a solution with the same solver but which has a nested schur complement
+    in the trace solve operator and a jacobi preconditioner on the A00 block.
+
+    NOTE: With the setup in this test, using jacobi as a preconditioner to the
+    schur complement matrix, the condition number of the matrix of the local solve
+    P.inv * A.solve(...) is reduced from 36.59 to 3.06.
+    """
+    # setup FEM
+    a, L, W = setup_poisson_3D()
+
+    # setup first solver
+    w = Function(W)
+    params = {'mat_type': 'matfree',
+              'ksp_type': 'preonly',
+              'pc_type': 'python',
+              'pc_python_type': 'firedrake.HybridizationPC',
+              'hybridization': {'ksp_type': 'cg',
+                                'pc_type': 'none',
+                                'ksp_rtol': 1e-12,
+                                'mat_type': 'matfree',
+                                'localsolve': {'ksp_type': 'preonly',
+                                               'pc_type': 'fieldsplit',
+                                               'pc_fieldsplit_type': 'schur',
+                                               'fieldsplit_0': {'ksp_type': 'default',
+                                                                'pc_type': 'jacobi'}}}}
+    eq = a == L
+    problem = LinearVariationalProblem(eq.lhs, eq.rhs, w)
+    solver = LinearVariationalSolver(problem, solver_parameters=params)
+    solver.solve()
+
+    # double-check options are set as expected
+    expected = {'nested': True,
+                'preonly_A00': False, 'jacobi_A00': True,
+                'schur_approx': False,
+                'preonly_Shat': False, 'jacobi_Shat': False,
+                'preonly_S': False, 'jacobi_S': False}
+    builder = solver.snes.ksp.pc.getPythonContext().getSchurComplementBuilder()
+    assert options_check(builder, expected), "Some solver options have not ended up in the PC as wanted."
+
+    sigma_h, u_h = w.split()
+
+    # setup second solver
+    w2 = Function(W)
+    solve(a == L, w2,
+          solver_parameters={'mat_type': 'matfree',
+                             'ksp_type': 'preonly',
+                             'pc_type': 'python',
+                             'pc_python_type': 'firedrake.HybridizationPC',
+                             'hybridization': {'ksp_type': 'cg',
+                                               'pc_type': 'none',
+                                               'ksp_rtol': 1e-8,
+                                               'mat_type': 'matfree'}})
+    nh_sigma, nh_u = w2.split()
+
+    # Return the L2 error
+    sigma_err = errornorm(sigma_h, nh_sigma)
+    u_err = errornorm(u_h, nh_u)
+    assert sigma_err < 1e-8
+    assert u_err < 1e-8
+
+
+def test_slate_hybridization_jacobi_prec_schur():
+    """A test, which compares a solution to a 3D mixed Poisson problem solved
+    globally matrixfree with a HybridizationPC and CG on the trace system to
+    a solution with the same solver but which has a nested schur complement
+    in the trace solve operator and a jacobi preconditioner on the schur
+    complement.
+
+    NOTE With the setup in this test, using jacobi as apreconditioner to the
+    schur complement matrix the condition number of the matrix of the local solve
+    P.inv * A.solve(...) is reduced from 17.13 to 16.71
+    """
+    # setup FEM
+    a, L, W = setup_poisson_3D()
+
+    # setup first solver
+    w = Function(W)
+    params = {'mat_type': 'matfree',
+              'ksp_type': 'preonly',
+              'pc_type': 'python',
+              'pc_python_type': 'firedrake.HybridizationPC',
+              'hybridization': {'ksp_type': 'cg',
+                                'pc_type': 'none',
+                                'ksp_rtol': 1e-12,
+                                'mat_type': 'matfree',
+                                'localsolve': {'ksp_type': 'preonly',
+                                               'pc_type': 'fieldsplit',
+                                               'pc_fieldsplit_type': 'schur',
+                                               'fieldsplit_1': {'ksp_type': 'default',
+                                                                'pc_type': 'jacobi'}}}}
+    eq = a == L
+    problem = LinearVariationalProblem(eq.lhs, eq.rhs, w)
+    solver = LinearVariationalSolver(problem, solver_parameters=params)
+    solver.solve()
+
+    # double-check options are set as expected
+    expected = {'nested': True,
+                'preonly_A00': False, 'jacobi_A00': False,
+                'schur_approx': False,
+                'preonly_S': False, 'jacobi_S': True}
+    builder = solver.snes.ksp.pc.getPythonContext().getSchurComplementBuilder()
+    assert options_check(builder, expected), "Some solver options have not ended up in the PC as wanted."
+
+    sigma_h, u_h = w.split()
+
+    # setup second solver
+    w2 = Function(W)
+    solve(a == L, w2,
+          solver_parameters={'mat_type': 'matfree',
+                             'ksp_type': 'preonly',
+                             'pc_type': 'python',
+                             'pc_python_type': 'firedrake.HybridizationPC',
+                             'hybridization': {'ksp_type': 'cg',
+                                               'pc_type': 'none',
+                                               'ksp_rtol': 1e-8,
+                                               'mat_type': 'matfree'}})
+    nh_sigma, nh_u = w2.split()
+
+    # Return the L2 error
+    sigma_err = errornorm(sigma_h, nh_sigma)
+    u_err = errornorm(u_h, nh_u)
+
+    assert sigma_err < 1e-8
+    assert u_err < 1e-8
+
+
+def test_mixed_poisson_approximated_schur_jacobi_prec():
+    """A test, which compares a solution to a 2D mixed Poisson problem solved
+    globally matrixfree with a HybridizationPC and CG on the trace system where
+    the a user supplied operator is used as preconditioner to the
+    Schur solver to a solution where the user supplied operator is replaced
+    with the jacobi preconditioning operator.
+    """
+    # setup FEM
+    a, L, W = setup_poisson()
+
+    # setup first solver
+    w = Function(W)
+    params = {'ksp_type': 'preonly',
+              'pc_type': 'python',
+              'mat_type': 'matfree',
+              'pc_python_type': 'firedrake.HybridizationPC',
+              'hybridization': {'ksp_type': 'cg',
+                                'pc_type': 'none',
+                                'ksp_rtol': 1e-8,
+                                'mat_type': 'matfree',
+                                'localsolve': {'ksp_type': 'preonly',
+                                               'pc_type': 'fieldsplit',
+                                               'pc_fieldsplit_type': 'schur',
+                                               'fieldsplit_1': {'ksp_type': 'default',
+                                                                'pc_type': 'python',
+                                                                'pc_python_type': __name__ + '.DGLaplacian',
+                                                                'aux_ksp_type': 'preonly',
+                                                                'aux_pc_type': 'jacobi'}}}}
+
+    eq = a == L
+    problem = LinearVariationalProblem(eq.lhs, eq.rhs, w)
+    solver = LinearVariationalSolver(problem, solver_parameters=params)
+    solver.solve()
+
+    # double-check options are set as expected
+    expected = {'nested': True,
+                'preonly_A00': False, 'jacobi_A00': False,
+                'schur_approx': True,
+                'preonly_S': False, 'jacobi_S': False,
+                'preonly_Shat': True, 'jacobi_Shat': True}
+    builder = solver.snes.ksp.pc.getPythonContext().getSchurComplementBuilder()
+    assert options_check(builder, expected), "Some solver options have not ended up in the PC as wanted."
+
+    sigma_h, u_h = w.split()
+
+    # setup second solver
+    w2 = Function(W)
+    aij_params = {'ksp_type': 'preonly',
+                  'pc_type': 'python',
+                  'mat_type': 'matfree',
+                  'pc_python_type': 'firedrake.HybridizationPC',
+                  'hybridization': {'ksp_type': 'cg',
+                                    'pc_type': 'none',
+                                    'ksp_rtol': 1e-8,
+                                    'mat_type': 'matfree'}}
+    solve(a == L, w2, solver_parameters=aij_params)
+    _sigma, _u = w2.split()
+
+    # Return the L2 error
+    sigma_err = errornorm(sigma_h, _sigma)
+    u_err = errornorm(u_h, _u)
+
+    assert sigma_err < 1e-8
+    assert u_err < 1e-8