diff --git a/genetIC/src/ic.hpp b/genetIC/src/ic.hpp
index b95edab7..6e9bcb53 100644
--- a/genetIC/src/ic.hpp
+++ b/genetIC/src/ic.hpp
@@ -1659,7 +1659,7 @@ class ICGenerator {
   }
 
   //! Splicing: fixes the flagged region, while reinitialising the exterior from a new random field
-  virtual void splice(size_t newSeed) {
+  virtual void splice(size_t newSeed, T accuracy) {
     initialiseRandomComponentIfUninitialised();
     if(outputFields.size()>1)
       throw std::runtime_error("Splicing is not yet implemented for the case of multiple transfer functions");
@@ -1680,14 +1680,16 @@ class ICGenerator {
     newGenerator.draw();
     logging::entry() << "Finished constructing new random field. Beginning splice operation." << endl;
 
-    for(size_t level=0; level<multiLevelContext.getNumLevels(); ++level) {
-      auto &originalFieldThisLevel = outputFields[0]->getFieldForLevel(level);
-      auto &newFieldThisLevel = newField.getFieldForLevel(level);
-      auto splicedFieldThisLevel = modifications::spliceOneLevel(newFieldThisLevel, originalFieldThisLevel,
-                                                             *multiLevelContext.getCovariance(level, particle::species::all));
-      splicedFieldThisLevel.toFourier();
-      originalFieldThisLevel = std::move(splicedFieldThisLevel);
-    }
+    modifications::splice(newField, *(outputFields[0]), accuracy);
+
+    // for(size_t level=0; level<multiLevelContext.getNumLevels(); ++level) {
+    //   auto &originalFieldThisLevel = outputFields[0]->getFieldForLevel(level);
+    //   auto &newFieldThisLevel = newField.getFieldForLevel(level);
+    //   auto splicedFieldThisLevel = modifications::spliceOneLevel(newFieldThisLevel, originalFieldThisLevel,
+    //                                                          *multiLevelContext.getCovariance(level, particle::species::all));
+    //   splicedFieldThisLevel.toFourier();
+    //   originalFieldThisLevel = std::move(splicedFieldThisLevel);
+    // }
   }
 
   //! Reverses the sign of the low-k modes.
diff --git a/genetIC/src/simulation/field/multilevelfield.hpp b/genetIC/src/simulation/field/multilevelfield.hpp
index d4593f5a..8f20cbdb 100644
--- a/genetIC/src/simulation/field/multilevelfield.hpp
+++ b/genetIC/src/simulation/field/multilevelfield.hpp
@@ -168,6 +168,11 @@ namespace fields {
       addScaled(other, 1.0);
     }
 
+    //! Subtracts the specified multi-level field to this one.
+    void operator-=(const MultiLevelField<DataType> &other) {
+      addScaled(other, -1.0);
+    }
+
     //! Divides the field by the specified ratio.
     void operator/=(DataType ratio) {
       using namespace tools::numerics;
@@ -178,6 +183,15 @@ namespace fields {
       }
     }
 
+    void operator*=(DataType val) {
+      using namespace tools::numerics;
+
+      for (size_t level = 0; level < getNumLevels(); level++) {
+        auto &data = getFieldForLevel(level).getDataVector();
+        data *= val;
+      }
+    }
+
     //! Flips the sign of the field.
     void reverse() {
       for (size_t level = 0; level < getNumLevels(); ++level) {
@@ -193,19 +207,17 @@ namespace fields {
     //! Add a scaled multilevel field to the current one
     void addScaled(const MultiLevelField &other, DataType scale) {
       assertContextConsistent();
-      assert(other.isFourierOnAllLevels());
+      if (other.isFourierOnAllLevels()) toFourier();
+      else if (other.isRealOnAllLevels()) toReal();
+      else throw std::runtime_error("Incompatible field types");
+      
       assert (isCompatible(other));
-      toFourier();
 
       for (size_t level = 0; level < getNumLevels(); level++) {
         if (hasFieldForLevel(level) && other.hasFieldForLevel(level)) {
           Field<DataType> &fieldThis = getFieldForLevel(level);
           const Field<DataType> &fieldOther = other.getFieldForLevel(level);
-          T kMin = fieldThis.getGrid().getFourierKmin();
-          fieldThis.forEachFourierCellInt([&fieldOther, kMin, scale]
-                                            (ComplexType currentVal, int kx, int ky, int kz) {
-            return currentVal + scale * fieldOther.getFourierCoefficient(kx, ky, kz);
-          });
+          fieldThis.addScaled(fieldOther, scale);
         }
       }
 
@@ -569,6 +581,14 @@ namespace fields {
     }
 
 
+    //! Returns a constant reference to the field on the specified grid, assuming it has been populated.
+    const Field<DataType, T> &getFieldForLevel(size_t i) const override {
+      this->assertContextConsistent();
+      assert(i < this->fieldsOnLevels.size());
+
+      return *(this->fieldsOnLevels[i]);
+    }
+
   };
 
 
diff --git a/genetIC/src/simulation/modifications/splice.hpp b/genetIC/src/simulation/modifications/splice.hpp
index 2cda17f4..5ca6f6c6 100644
--- a/genetIC/src/simulation/modifications/splice.hpp
+++ b/genetIC/src/simulation/modifications/splice.hpp
@@ -4,6 +4,7 @@
 #include <complex>
 #include <src/tools/data_types/complex.hpp>
 #include <src/tools/numerics/cg.hpp>
+#include <src/tools/numerics/minres.hpp>
 
 namespace modifications {
   template<typename T>
@@ -33,86 +34,289 @@ namespace modifications {
     return mask;
   }
 
-  //! Return the field f which satisfies f = a in flagged region while minimising (f-b).C^-1.(f-b) elsewhere
+  // Compute the operator that applies T^+ ... T *at all levels*
   template<typename DataType, typename T=tools::datatypes::strip_complex<DataType>>
-  fields::Field<DataType,T> spliceOneLevel(fields::Field<DataType,T> & a,
-                                           fields::Field<DataType,T> & b,
-                                           const fields::Field<DataType,T> & cov) {
+  fields::OutputField<DataType> T_op_T (
+      const fields::OutputField<DataType> &inputs,
+      const std::vector<fields::Field<DataType, T>> &covs,
+      const std::function<void(const int, fields::Field<DataType, T> &)> op
+  ) {
+    fields::OutputField<DataType> outputs(inputs.getContext(), inputs.getTransferType());
+    outputs.getFieldForLevel(0); // trigger allocation
+
+    fields::OutputField<DataType> inputs_delta(inputs);
+    auto filters = inputs.getFilters();
+    int Nlevel = inputs.getNumLevels();
+
+    const auto& multiLevelContext = outputs.getContext();
+
+    assert (inputs.isRealOnAllLevels());
+
+    // Apply C^0.5 to the input fields
+    for (size_t level=0; level<Nlevel; ++level) {
+      auto& field = inputs_delta.getFieldForLevel(level);
+      field.toFourier();
+      field.applyTransferFunction(covs[level], 0.5);
+    }
+
+    // Filter all levels (but the last) in their windows
+    // Notice we are doing window then filter
+    for (size_t level=0; level<Nlevel; ++level) {
+      auto& field = inputs_delta.getFieldForLevel(level);
+      const auto &f = filters.getFilterForLevel(level);
+      if (level < Nlevel - 1) {
+        auto window = multiLevelContext.getGridForLevel(level+1).getWindow();
+        field.applyFilterInWindow(f, window, false);
+      } else {
+        field.toFourier();
+        field.applyFilter(f);
+      }
+    }
+
+    // --------------------------------------------------------------------
+    // Operator applies to the field + contributions from all finer levels
+    for (size_t level=0; level<Nlevel; ++level) {
+      auto out = inputs_delta.getFieldForLevel(level).copy();
+
+      // Add contributions from all coarser levels
+      for (size_t source_level = 0; source_level < level; ++source_level) {
+        auto source_field = inputs_delta.getFieldForLevel(source_level).copy();
+        T pixel_volume_ratio = multiLevelContext.getWeightForLevel(level) /
+                               multiLevelContext.getWeightForLevel(source_level);
+        
+        source_field->toReal();
+        *source_field *= sqrt(pixel_volume_ratio);
+
+        out->addFieldFromDifferentGrid(*source_field);
+      }
+
+      // Apply operator
+      op(level, *out);
+
+      // Add contributions from all finer levels
+      for (size_t source_level = level + 1; source_level < Nlevel; ++source_level) {
+        auto source_field = inputs_delta.getFieldForLevel(source_level).copy();
+        T pixel_volume_ratio = multiLevelContext.getWeightForLevel(level) /
+                               multiLevelContext.getWeightForLevel(source_level);
+
+        source_field->toFourier();
+        *source_field *= sqrt(pixel_volume_ratio);
+        op(source_level, *source_field);
+
+        out->addFieldFromDifferentGrid(*source_field);
+      }
+
+      outputs.getFieldForLevel(level) = std::move(*out);
+    }
+
+    // Filter all levels (but the last) in their windows
+    // Notice we are doing filter then window (the opposite order to above)
+    for (size_t level=0; level<Nlevel; ++level) {
+      auto& field = outputs.getFieldForLevel(level);
+      const auto &f = filters.getFilterForLevel(level);
+      if (level < Nlevel - 1) {
+        auto window = multiLevelContext.getGridForLevel(level+1).getWindow();
+        field.applyFilterInWindow(f, window, true);
+      } else {
+        field.toFourier();
+        field.applyFilter(f);
+      }
+    }
+    
+    // Apply C^0.5 to the input fields
+    for (size_t level=0; level<Nlevel; ++level) {
+      auto& field = outputs.getFieldForLevel(level);
+      field.toFourier();
+      field.applyTransferFunction(covs[level], 0.5);
+    }
+
+    outputs.toReal();
+    return outputs;
+  };
+
+  template<typename DataType, typename T=tools::datatypes::strip_complex<DataType>>
+  fields::OutputField<DataType> Mbar_Cm1_Mbar(
+    const fields::OutputField<DataType> & inputs,
+    const auto& covs, const auto& masks, const auto& masksCompl, size_t Nlevel
+  ) {
+    auto outputs = T_op_T<DataType, T>(
+      inputs,
+      covs,
+      [&](const int level, fields::Field<DataType,T> & input) -> void
+    {
+      input.toReal();
+      input *= masksCompl[level];
+      input.toFourier();
+      input.applyTransferFunction(covs[level], -1.0);
+      input.toReal();
+      input *= masksCompl[level];
+    });
+    return outputs;
+  }
+
+  template<typename DataType, typename T=tools::datatypes::strip_complex<DataType>>
+  fields::OutputField<DataType> Mbar_Cm1_M(
+    const fields::OutputField<DataType> & inputs,
+    const auto& covs, const auto& masks, const auto& masksCompl, size_t Nlevel
+  ) {
+    auto outputs = T_op_T<DataType, T>(
+      inputs, covs,
+      [&](const int level, fields::Field<DataType,T> & input) -> void
+    {
+      input.toReal();
+      input *= masks[level];
+      input.toFourier();
+      input.applyTransferFunction(covs[level], -1.0);
+      input.toReal();
+      input *= masksCompl[level];
+    });
+    return outputs;
+  };
+
+  template<typename DataType, typename T=tools::datatypes::strip_complex<DataType>>
+  fields::OutputField<DataType> combine(
+    const fields::OutputField<DataType> &a,
+    const std::vector<fields::Field<DataType, T>> &covs
+  ) {
+    fields::OutputField<DataType> inputs(a);
+    const auto& multiLevelContext = inputs.getContext();
+    auto filters = a.getFilters();
+    int Nlevel = a.getNumLevels();
+
+    for (size_t level=0; level<Nlevel; ++level) {
+      auto& field = inputs.getFieldForLevel(level);
+      const auto &f = filters.getFilterForLevel(level);
+      field.toFourier();
+      field.applyTransferFunction(covs[level], 0.5);
+    }
+
+    // Copy inputs
+    fields::OutputField<DataType> outputs(inputs);
+
+    // Add contribution from coarser level
+    for (size_t level = 1; level<Nlevel; ++level) {
+      auto & out = outputs.getFieldForLevel(level);
+
+      // Remove low-frequency information from this level
+      out.toFourier();
+      out.applyFilter(filters.getHighPassFilterForLevel(level));
+
+      // Replace with the low-frequency information from the level below
+      out.addFieldFromDifferentGridWithFilter(
+        inputs.getFieldForLevel(level - 1),
+        filters.getLowPassFilterForLevel(level - 1)
+      );
+    }
+    outputs.toReal();
+    outputs.getContext().setLevelsAreCombined();
+    return outputs;
+  }
+
+  template<typename DataType, typename T=tools::datatypes::strip_complex<DataType>>
+  fields::OutputField<DataType> splice(fields::OutputField<DataType> & a,
+                                       fields::OutputField<DataType> & b,
+                                       T accuracy) {
+
+      assert (a.getTransferType() == particle::species::whitenoise);
+      assert (b.getTransferType() == particle::species::whitenoise);
+      assert (a.isFourierOnAllLevels());
+      assert (b.isFourierOnAllLevels());
 
       // To understand the implementation below, first read Appendix A of Cadiou et al (2021),
       // and/or look at the 1D toy implementation (in tools/toy_implementation/gene_splicing.ipynb) which
       // contains a similar derivation and near-identical implementation.
 
-      assert (&a.getGrid() == &b.getGrid());
-      assert (&a.getGrid() == &cov.getGrid());
-      auto mask = generateMaskFromFlags(a.getGrid());
-      auto maskCompl = generateMaskComplementFromFlags(a.getGrid());
-
-      a.toFourier();
-      b.toFourier();
-
-      // The preconditioner should be almost equal to the covariance.
-      // We however set the fundamental of the power spectrum to a non-null value,
-      // otherwise, the mean value in the spliced region is unconstrained.
-      fields::Field<DataType,T> preconditioner(cov);
-      preconditioner.setFourierCoefficient(0, 0, 0, 1);
-
-      fields::Field<DataType,T> delta_diff = b-a;
-      delta_diff.applyTransferFunction(preconditioner, 0.5);
-      delta_diff.toReal();
-
-
-      fields::Field<DataType,T> z(delta_diff);
-      z*=mask;
-      z.toFourier();
-      z.applyTransferFunction(preconditioner, -1.0);
-      z.toReal();
-      z*=maskCompl;
-      z.toFourier();
-      z.applyTransferFunction(preconditioner, 0.5);
-      z.toReal();
-
-      auto X = [&preconditioner, &maskCompl](const fields::Field<DataType,T> & input) -> fields::Field<DataType,T>
-      {
-        fields::Field<DataType,T> v(input);
-        assert (!input.isFourier());
-        assert (!v.isFourier());
-        v.toFourier();
-        v.applyTransferFunction(preconditioner, 0.5);
-        v.toReal();
-        v*=maskCompl;
-        v.toFourier();
-        v.applyTransferFunction(preconditioner, -1.0);
-        v.toReal();
-        v*=maskCompl;
-        v.toFourier();
-        v.applyTransferFunction(preconditioner, 0.5);
-        v.toReal();
-        return v;
+      std::vector<fields::Field<DataType,T>> covs;
+      std::vector<fields::Field<char,T>> masks;
+      std::vector<fields::Field<char,T>> masksCompl;
+
+      int Nlevel = a.getNumLevels();
+      for(size_t level=0; level<Nlevel; ++level) {
+        auto ctxt = a.getContext();
+        fields::Field<DataType,T> cov(*ctxt.getCovariance(level, particle::species::all));
+        cov.setFourierCoefficient(0, 0, 0, 1);
+        covs.push_back(cov);
+
+        masks.push_back(generateMaskFromFlags(ctxt.getGridForLevel(level)));
+        masksCompl.push_back(generateMaskComplementFromFlags(ctxt.getGridForLevel(level)));
+      }
+
+      a.toReal();
+      b.toReal();
+
+      fields::OutputField<T> delta(b);
+      delta -= a;
+
+      fields::OutputField<T> z = Mbar_Cm1_M(delta, covs, masks, masksCompl, Nlevel);
+
+      auto A = [&](const fields::OutputField<T> & inputs) {
+        return Mbar_Cm1_Mbar(inputs, covs, masks, masksCompl, Nlevel);
       };
 
 
-      fields::Field<DataType,T> alpha = tools::numerics::conjugateGradient<DataType>(X,z);
+      // fields::OutputField<DataType> alpha = tools::numerics::conjugateGradient<DataType>(A, z);
+      fields::OutputField<DataType> alpha = tools::numerics::minres<DataType>(A, z, accuracy);
+
+      // Combine fields
+      alpha = combine(alpha, covs);
+      a = combine(a, covs);
+      b = combine(b, covs);
+      a.toReal(); b.toReal(); alpha.toReal();
+
+      // output = b + M(a - b) + Mbar alpha [all in delta basis now]
+      fields::OutputField<DataType> outputs(b.getContext(), particle::species::all);
+      outputs.getFieldForLevel(0); // trigger allocation
+      outputs.toReal();
+      
+      for (size_t level = 0; level < Nlevel; ++level) {
+        const auto & a_field = a.getFieldForLevel(level);
+        const auto & b_field = b.getFieldForLevel(level);
+        const auto & alpha_field = alpha.getFieldForLevel(level);
+
+        auto & out = outputs.getFieldForLevel(level);
+        out += b_field;
+
+        {
+          auto temp = a_field.copy();
+          *temp -= b_field;
+          *temp *= masks[level];
+
+          out += *temp;
+        }
+
+        {
+          auto temp = alpha_field.copy();
+          *temp *= masksCompl[level];
+
+          out += *temp;
+        }
+        out.toReal();
+
+        a_field.dumpGridData("a_" + std::to_string(level) + ".dat");
+        b_field.dumpGridData("b_" + std::to_string(level) + ".dat");
+        alpha_field.dumpGridData("alpha_" + std::to_string(level) + ".dat");
+        out.dumpGridData("splice_level_" + std::to_string(level) + ".dat");
+      }
 
-      alpha.toFourier();
-      alpha.applyTransferFunction(preconditioner, 0.5);
-      alpha.toReal();
+      // alpha.toFourier();
+      // alpha.applyTransferFunction(preconditioner, 0.5);
+      // alpha.toReal();
 
-      fields::Field<DataType,T> bInDeltaBasis(b);
-      bInDeltaBasis.toFourier();
-      bInDeltaBasis.applyTransferFunction(preconditioner, 0.5);
-      bInDeltaBasis.toReal();
+      // fields::Field<DataType,T> bInDeltaBasis(b);
+      // bInDeltaBasis.toFourier();
+      // bInDeltaBasis.applyTransferFunction(preconditioner, 0.5);
+      // bInDeltaBasis.toReal();
 
-      alpha*=maskCompl;
-      alpha+=bInDeltaBasis;
+      // alpha*=maskCompl;
+      // alpha+=bInDeltaBasis;
 
-      delta_diff*=mask;
-      alpha-=delta_diff;
+      // delta_diff*=mask;
+      // alpha-=delta_diff;
 
-      assert(!alpha.isFourier());
-      alpha.toFourier();
-      alpha.applyTransferFunction(preconditioner, -0.5);
-      alpha.toReal();
+      // assert(!alpha.isFourier());
+      // alpha.toFourier();
+      // alpha.applyTransferFunction(preconditioner, -0.5);
+      // alpha.toReal();
 
       return alpha;
   }
diff --git a/genetIC/src/tools/numerics/cg.hpp b/genetIC/src/tools/numerics/cg.hpp
index bce1cf18..2ac765e5 100644
--- a/genetIC/src/tools/numerics/cg.hpp
+++ b/genetIC/src/tools/numerics/cg.hpp
@@ -9,45 +9,69 @@
 namespace tools {
   namespace numerics {
 
+    template<typename T>
+    T innerProduct(const fields::OutputField<T> &a, const fields::OutputField<T> &b) {
+      T result = 0;
+      for (auto ilevel = 0; ilevel < a.getNumLevels(); ++ilevel) {
+        const auto& left = a.getFieldForLevel(ilevel);
+        const auto& right = b.getFieldForLevel(ilevel);
+        result += left.innerProduct(right);
+      }
+      return result;
+    }
+
+    template<typename T>
+    double norm(const fields::OutputField<T> &a) {
+      return std::sqrt(innerProduct(a, a));
+    }
+
     //! Solve linear equation Qx = b, and return x, using conjugate gradient
     template<typename T>
-    fields::Field<T> conjugateGradient(std::function<fields::Field<T>(const fields::Field<T> &)> Q,
-                                       const fields::Field<T> &b,
+    fields::OutputField<T> conjugateGradient(std::function<fields::OutputField<T>(const fields::OutputField<T> &)> Q,
+                                       const fields::OutputField<T> &b,
                                        double rtol = 1e-6,
                                        double atol = 1e-12) {
-      fields::Field<T> residual(b);
-      fields::Field<T> direction = -residual;
-      fields::Field<T> x = fields::Field<T>(b.getGrid(), false);
+      fields::OutputField<T> residual(b);
+      fields::OutputField<T> direction(residual);
+      direction *= -1;
+      fields::OutputField<T> x = fields::OutputField<T>(b.getContext(), b.getTransferType());
+      x.getFieldForLevel(0); // trigger allocation
 
-      double scale = b.norm();
+      double scale = norm(residual);
 
       if(scale==0.0) {
         logging::entry(logging::warning) << "Conjugate gradient: result is zero!" << std::endl;
         return x;
       }
 
-      size_t dimension = b.getGrid().size3;
+      size_t dimension = 0;
+      for (auto ilevel = 0; ilevel < b.getNumLevels(); ++ilevel) {
+        const auto ctxt = residual.getContext();
+        dimension += ctxt.getGridForLevel(ilevel).size3;
+      }
 
       size_t i;
 
       for(i=0; i<dimension+1; ++i) {
 
-        fields::Field<T> Q_direction = Q(direction);
+        auto Q_direction = Q(direction);
+  
         // distance to travel in specified direction
-        double alpha = -residual.innerProduct(direction) / direction.innerProduct(Q_direction);
+        double alpha = -innerProduct(residual, direction) / innerProduct(direction, Q_direction);
+
         x.addScaled(direction, alpha);
 
         residual = Q(x);
         residual -= b;
 
-        auto norm = residual.norm();
-        if (norm < rtol * scale || norm < atol)
+        auto res_norm = norm(residual);
+        if (res_norm < rtol * scale || res_norm < atol)
           break;
 
-        logging::entry() << "Conjugate gradient iteration " << i << " residual=" << norm << std::endl;
+        logging::entry() << "Conjugate gradient iteration " << i << " residual=" << res_norm << "/" << scale << std::endl;
 
         // update direction for next cycle; must be Q-orthogonal to all previous updates
-        double beta = residual.innerProduct(Q_direction) / direction.innerProduct(Q_direction);
+        double beta = innerProduct(residual, Q_direction) / innerProduct(direction, Q_direction);
         direction*=beta;
         direction-=residual;
 
diff --git a/genetIC/src/tools/numerics/minres.hpp b/genetIC/src/tools/numerics/minres.hpp
new file mode 100644
index 00000000..8c8d325e
--- /dev/null
+++ b/genetIC/src/tools/numerics/minres.hpp
@@ -0,0 +1,80 @@
+#ifndef IC_MINRES_HPP
+#define IC_MINRES_HPP
+
+#include <functional>
+#include <src/simulation/field/field.hpp>
+#include <src/tools/data_types/complex.hpp>
+#include <src/tools/logging.hpp>
+#include <ctime>
+#include <iostream>
+#include <fstream>
+#include <sys/stat.h>
+
+#include "cg.hpp"
+
+namespace tools {
+  namespace numerics {
+
+    /* Adapted from https://stanford.edu/group/SOL/reports/SOL-2011-2R.pdf */
+    template<typename T>
+    fields::OutputField<T> minres(std::function<fields::OutputField<T>(const fields::OutputField<T> &)> A,
+                            const fields::OutputField<T> &b,
+                            double rtol = 1e-6,
+                            double atol = 1e-12 )
+    {
+      fields::OutputField<T> x(b.getContext(), b.getTransferType());
+      x.getFieldForLevel(0);  // Trigger allocation
+      fields::OutputField<T> r(b);
+      fields::OutputField<T> s(b);
+      s = A(r);
+
+      fields::OutputField<T> p(r);
+      fields::OutputField<T> q(s);
+
+      double scale = tools::numerics::norm(b);
+      double rho = tools::numerics::innerProduct(r, s);
+
+      size_t dimension = 0;
+      for (auto ilevel = 0; ilevel < b.getNumLevels(); ++ilevel) {
+        const auto ctxt = r.getContext();
+        dimension += ctxt.getGridForLevel(ilevel).size3;
+      }
+
+      size_t max_iterations = dimension * 10;
+      double old_norm = 0;
+
+      size_t iter = 0;
+      
+      for(; iter<max_iterations; ++iter) {
+        // We have q = A(p), but no need to compute it again
+        double alpha = rho / tools::numerics::innerProduct(q, q);
+        x.addScaled(p, alpha);
+        r.addScaled(q, -alpha);
+
+        double norm = tools::numerics::norm(r);
+
+        if (norm < rtol * scale || norm < atol)
+          break;
+
+        logging::entry() << "MINRES iteration " << iter << " residual=" << norm/scale << std::endl;
+
+        s = A(r);
+        double rhobar = rho;
+        rho = tools::numerics::innerProduct(r, s);
+        double beta = rho / rhobar;
+        p *= beta;
+        p += r;
+
+        q *= beta;
+        q += s;
+      }
+     
+
+      logging::entry() << "MINRES ended after " << iter << " iterations" << std::endl;
+
+      return x;
+    }
+  }
+}
+
+#endif
\ No newline at end of file