diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc
index 03ca896b4aa..7190a888302 100644
--- a/source/simulator/solver.cc
+++ b/source/simulator/solver.cc
@@ -27,15 +27,19 @@
 #include <aspect/mesh_deformation/interface.h>
 
 #include <deal.II/base/signaling_nan.h>
+#include <deal.II/lac/diagonal_matrix.h>
+#include <deal.II/lac/linear_operator.h>
 #include <deal.II/lac/solver_gmres.h>
 #include <deal.II/lac/solver_bicgstab.h>
 #include <deal.II/lac/solver_cg.h>
 #include <deal.II/fe/fe_values.h>
+#include <deal.II/lac/trilinos_vector.h>
 
 namespace aspect
 {
   namespace internal
   {
+
     /**
      * Implement multiplication with Stokes part of system matrix. In essence, this
      * object represents a 2x2 block matrix that corresponds to the top left
@@ -165,6 +169,7 @@ namespace aspect
     }
 
 
+
     /**
      * Implement the block Schur preconditioner for the Stokes system.
      */
@@ -381,6 +386,32 @@ namespace aspect
 
     };
 
+
+    /**
+     * Given a diagonal matrix stored as a vector,
+     * create an operator that represents its action.
+     */
+
+    template<typename Range,
+             typename Domain,
+             typename Payload>
+    LinearOperator<Range, Domain, Payload> diag_operator(const LinearOperator<Range,Domain,Payload> &exemplar,
+                                                         const TrilinosWrappers::MPI::Vector &diagonal)
+    {
+      LinearOperator<Range, Domain, Payload> return_op;
+      return_op.reinit_range_vector  = exemplar.reinit_range_vector;
+      return_op.reinit_domain_vector = exemplar.reinit_domain_vector;
+
+      return_op.vmult = [&](Range &dest, const Domain &src)
+      {
+        dest = src;
+        dest.scale(diagonal);
+      };
+      return return_op;
+    }
+
+
+
     /**
      * This class approximates the Schur Complement inverse operator
      * by S^{-1} = (BC^{-1}B^T)^{-1}(BC^{-1}AD^{-1}B^T)(BD^{-1}B^T)^{-1},
@@ -458,7 +489,13 @@ namespace aspect
             SolverControl solver_control(5000, 1e-6 * src.l2_norm(), false, true);
             SolverCG<TrilinosWrappers::MPI::Vector> solver(solver_control);
             //Solve with Schur Complement approximation
-            solver.solve(mp_matrix,
+            const auto Op_A      = LinearOperator<TrilinosWrappers::MPI::Vector>(system_matrix.block(0,0));
+            const auto Op_BT     = LinearOperator<TrilinosWrappers::MPI::Vector>(system_matrix.block(0,1));
+            const auto Op_B      = LinearOperator<TrilinosWrappers::MPI::Vector>(system_matrix.block(1,0));
+            const auto Op_C_inv  = diag_operator(Op_A,inverse_lumped_mass_matrix);
+            const auto BC_invBT  = Op_B*Op_C_inv*Op_BT;
+
+            solver.solve(BC_invBT,
                          ptmp,
                          src,
                          mp_preconditioner);
@@ -471,7 +508,7 @@ namespace aspect
             system_matrix.block(1,0).vmult(ptmp,wtmp);
 
             dst=0;
-            solver.solve(mp_matrix,
+            solver.solve(BC_invBT,
                          dst,
                          ptmp,
                          mp_preconditioner);