diff --git a/dune/solvers/common/arithmetic.hh b/dune/solvers/common/arithmetic.hh
index 664bb5c573fb8beef687be36bbf20af43703cf94..6461e881b246f4a7465e8ec3d73600672cf5ab4a 100644
--- a/dune/solvers/common/arithmetic.hh
+++ b/dune/solvers/common/arithmetic.hh
@@ -5,6 +5,7 @@
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>
 #include <dune/common/diagonalmatrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/scaledidmatrix.hh>
 #include <dune/common/typetraits.hh>
 
@@ -112,6 +113,10 @@ namespace Arithmetic
         {
             b.umv(c, a);
         }
+        static void subtractProduct(A& a, const B& b, const C& c)
+        {
+            b.mmv(c, a);
+        }
     };
 
     // Internal helper class for scaled product operations (i.e., b is always a scalar)
@@ -122,6 +127,10 @@ namespace Arithmetic
         {
             c.usmv(b, d, a);
         }
+        static void subtractProduct(A& a, const B& b, const C& c, const D& d)
+        {
+            c.usmv(-b, d, a);
+        }
     };
 
     template<class T, int n, int m, int k>
@@ -141,6 +150,17 @@ namespace Arithmetic
                 }
             }
         }
+        static void subtractProduct(A& a, const B& b, const C& c)
+        {
+            for (size_t row = 0; row < n; ++row)
+            {
+                for (size_t col = 0 ; col < k; ++col)
+                {
+                    for (size_t i = 0; i < m; ++i)
+                        a[row][col] -= b[row][i]*c[i][col];
+                }
+            }
+        }
     };
 
     template<class T, int n, int m, int k>
@@ -158,6 +178,15 @@ namespace Arithmetic
                 }
             }
         }
+        static void subtractProduct(A& a, const T& b, const C& c, const D& d)
+        {
+            for (size_t row = 0; row < n; ++row) {
+                for (size_t col = 0 ; col < k; ++col) {
+                    for (size_t i = 0; i < m; ++i)
+                        a[row][col] -= b * c[row][i]*d[i][col];
+                }
+            }
+        }
     };
 
     template<class T, int n>
@@ -171,6 +200,11 @@ namespace Arithmetic
             for (size_t i=0; i<n; ++i)
                 a.diagonal(i) += b.diagonal(i)*c.diagonal(i);
         }
+        static void subtractProduct(A& a, const B& b, const C& c)
+        {
+            for (size_t i=0; i<n; ++i)
+                a.diagonal(i) -= b.diagonal(i)*c.diagonal(i);
+        }
     };
 
     template<class T, int n>
@@ -184,6 +218,11 @@ namespace Arithmetic
             for (size_t i=0; i<n; ++i)
                 a.diagonal(i) += b * c.diagonal(i)*d.diagonal(i);
         }
+        static void subtractProduct(A& a, const T& b, const C& c, const D& d)
+        {
+            for (size_t i=0; i<n; ++i)
+                a.diagonal(i) -= b * c.diagonal(i)*d.diagonal(i);
+        }
     };
 
     template<class T, int n>
@@ -197,6 +236,11 @@ namespace Arithmetic
             for (size_t i=0; i<n; ++i)
                 a[i][i] += b.diagonal(i)*c.diagonal(i);
         }
+        static void subtractProduct(A& a, const B& b, const C& c)
+        {
+            for (size_t i=0; i<n; ++i)
+                a[i][i] -= b.diagonal(i)*c.diagonal(i);
+        }
     };
 
     template<class T, int n>
@@ -210,6 +254,11 @@ namespace Arithmetic
             for (size_t i=0; i<n; ++i)
                 a[i][i] += b * c.diagonal(i)*d.diagonal(i);
         }
+        static void subtractProduct(A& a, const T& b, const C& c, const D& d)
+        {
+            for (size_t i=0; i<n; ++i)
+                a[i][i] -= b * c.diagonal(i)*d.diagonal(i);
+        }
     };
 
     /** \brief Specialization for b being a scalar type
@@ -221,6 +270,10 @@ namespace Arithmetic
         {
             a.axpy(b, c);
         }
+        static void subtractProduct(A& a, const B& b, const C& c)
+        {
+            a.axpy(-b, c);
+        }
     };
 
     /** \brief Specialization for c being a scalar type
@@ -232,6 +285,10 @@ namespace Arithmetic
         {
             a.axpy(b * c, d);
         }
+        static void subtractProduct(A& a, const B& b, const C& c, const D& d)
+        {
+            a.axpy(-b * c, d);
+        }
     };
 
     template<class A, class B, class C>
@@ -241,6 +298,10 @@ namespace Arithmetic
         {
             a += b*c;
         }
+        static void subtractProduct(A& a, const B& b, const C& c)
+        {
+            a -= b*c;
+        }
     };
 
     template<class A, class B, class C, class D>
@@ -250,6 +311,10 @@ namespace Arithmetic
         {
             a += b*c*d;
         }
+        static void subtractProduct(A& a, const B& b, const C& c, const D& d)
+        {
+            a -= b*c*d;
+        }
     };
 
     /** \brief Add a product to some matrix or vector
@@ -269,9 +334,26 @@ namespace Arithmetic
         ProductHelper<A,B,C,ScalarTraits<A>::isScalar, ScalarTraits<B>::isScalar, ScalarTraits<C>::isScalar>::addProduct(a,b,c);
     }
 
+    /** \brief Subtract a product from some matrix or vector
+     *
+     * This function computes a-=b*c.
+     *
+     * This function should tolerate all meaningful
+     * combinations of scalars, vectors, and matrices.
+     *
+     * a,b,c could be matrices with appropriate
+     * dimensions. But b can also always be a scalar
+     * represented by a 1-dim vector or a 1 by 1 matrix.
+     */
+    template<class A, class B, class C>
+    void subtractProduct(A& a, const B& b, const C& c)
+    {
+        ProductHelper<A,B,C,ScalarTraits<A>::isScalar, ScalarTraits<B>::isScalar, ScalarTraits<C>::isScalar>::subtractProduct(a,b,c);
+    }
+
     /** \brief Add a scaled product to some matrix or vector
      *
-     * This function computes a+=b*c.
+     * This function computes a+=b*c*d.
      *
      * This function should tolerate all meaningful
      * combinations of scalars, vectors, and matrices.
@@ -286,6 +368,155 @@ namespace Arithmetic
     {
         ScaledProductHelper<A,B,C,D,ScalarTraits<A>::isScalar, ScalarTraits<C>::isScalar, ScalarTraits<D>::isScalar>::addProduct(a,b,c,d);
     }
+
+    /** \brief Subtract a scaled product from some matrix or vector
+     *
+     * This function computes a-=b*c*d.
+     *
+     * This function should tolerate all meaningful
+     * combinations of scalars, vectors, and matrices.
+     *
+     * a,c,d could be matrices with appropriate dimensions. But b must
+     * (currently) and c can also always be a scalar represented by a
+     * 1-dim vector or a 1 by 1 matrix.
+     */
+    template<class A, class B, class C, class D>
+    typename Dune::enable_if<ScalarTraits<B>::isScalar, void>::type
+    subtractProduct(A& a, const B& b, const C& c, const D& d)
+    {
+        ScaledProductHelper<A,B,C,D,ScalarTraits<A>::isScalar, ScalarTraits<C>::isScalar, ScalarTraits<D>::isScalar>::subtractProduct(a,b,c,d);
+    }
+
+    namespace {
+        template <class MatrixType, class VectorType>
+        typename VectorType::field_type
+        AxyWithTemporary(const MatrixType &A,
+                         const VectorType &x,
+                         const VectorType &y)
+        {
+            VectorType tmp(y.size());
+            tmp = 0.0;
+            addProduct(tmp, A, x);
+            return tmp * y;
+        }
+
+        template <class MatrixType, class VectorType>
+        typename VectorType::field_type
+        AxyWithoutTemporary(const MatrixType &A,
+                            const VectorType &x,
+                            const VectorType &y)
+        {
+            assert(x.N() == A.M());
+            assert(y.N() == A.N());
+
+            typename VectorType::field_type outer = 0.0;
+            typename VectorType::block_type inner;
+            typename MatrixType::ConstIterator endi=A.end();
+            for (typename MatrixType::ConstIterator i=A.begin(); i!=endi; ++i) {
+                inner = 0.0;
+                const size_t iindex = i.index();
+                typename MatrixType::row_type::ConstIterator endj = i->end();
+                for (typename MatrixType::row_type::ConstIterator j=i->begin();
+                     j!=endj; ++j)
+                    addProduct(inner, *j, x[j.index()]);
+                outer += inner * y[iindex];
+            }
+            return outer;
+        }
+    }
+
+    //! Compute \f$(Ax,y)\f$
+    template <class MatrixType, class VectorType>
+    typename VectorType::field_type
+    Axy(const MatrixType &A,
+        const VectorType &x,
+        const VectorType &y)
+    {
+        return AxyWithTemporary(A, x, y);
+    }
+
+    //! Compute \f$(Ax,y)\f$
+    template <int n, typename field_type, class VectorType>
+    typename VectorType::field_type
+    Axy(const Dune::FieldMatrix<field_type, n, n> &A,
+        const VectorType &x,
+        const VectorType &y)
+    {
+        return AxyWithoutTemporary(A, x, y);
+    }
+
+    //! Compute \f$(Ax,y)\f$
+    template <typename block_type, class VectorType>
+    typename VectorType::field_type
+    Axy(const Dune::BCRSMatrix<block_type> &A,
+        const VectorType &x,
+        const VectorType &y)
+    {
+        return AxyWithoutTemporary(A, x, y);
+    }
+
+    namespace {
+        template <class MatrixType, class VectorType>
+        typename VectorType::field_type
+        bmAxyWithTemporary(const MatrixType &A, const VectorType &b,
+                           const VectorType &x, const VectorType &y)
+        {
+            VectorType tmp = b;
+            subtractProduct(tmp, A, x);
+            return tmp * y;
+        }
+
+        template <class MatrixType, class VectorType>
+        typename VectorType::field_type
+        bmAxyWithoutTemporary(const MatrixType &A, const VectorType &b,
+                              const VectorType &x, const VectorType &y)
+        {
+            assert(x.N() == A.M());
+            assert(y.N() == A.N());
+            assert(y.N() == b.N());
+
+            typename VectorType::field_type outer = 0.0;
+            typename VectorType::block_type inner;
+            typename MatrixType::ConstIterator endi=A.end();
+            for (typename MatrixType::ConstIterator i=A.begin(); i!=endi; ++i) {
+                const size_t iindex = i.index();
+                inner = b[iindex];
+                typename MatrixType::row_type::ConstIterator endj = i->end();
+                for (typename MatrixType::row_type::ConstIterator j=i->begin();
+                     j!=endj; ++j)
+                    subtractProduct(inner, *j, x[j.index()]);
+                outer += inner * y[iindex];
+            }
+            return outer;
+        }
+    }
+
+    //! Compute \f$(b-Ax,y)\f$
+    template <class MatrixType, class VectorType>
+    typename VectorType::field_type
+    bmAxy(const MatrixType &A, const VectorType &b,
+          const VectorType &x, const VectorType &y)
+    {
+        return bmAxyWithTemporary(A, b, x, y);
+    }
+
+    //! Compute \f$(b-Ax,y)\f$
+    template <int n, typename field_type, class VectorType>
+    typename VectorType::field_type
+    bmAxy(const Dune::FieldMatrix<field_type, n, n> &A, const VectorType &b,
+          const VectorType &x, const VectorType &y)
+    {
+        return bmAxyWithoutTemporary(A, b, x, y);
+    }
+
+    //! Compute \f$(b-Ax,y)\f$
+    template <typename block_type, class VectorType>
+    typename VectorType::field_type
+    bmAxy(const Dune::BCRSMatrix<block_type> &A, const VectorType &b,
+          const VectorType &x, const VectorType &y)
+    {
+        return bmAxyWithoutTemporary(A, b, x, y);
+    }
 };
 
 #endif