diff --git a/dune/istl/bcrsmatrix.hh b/dune/istl/bcrsmatrix.hh
index e36a20df34d90d9fca41d019ef9d4cf45a54be17..bc8bb0fac2a9149138bbfa75e2c2a42498cfe45f 100644
--- a/dune/istl/bcrsmatrix.hh
+++ b/dune/istl/bcrsmatrix.hh
@@ -963,6 +963,29 @@ namespace Dune {
 
       return *this;
     }
+
+    /*! \brief Add the scaled entries of another matrix to this one.
+     *
+     * Matrix axpy operation: *this += alpha * b
+     *
+     * \param alpha Scaling factor.
+     * \param b     The matrix to add to this one. Its sparsity pattern has to
+     *              be subset of the sparsity pattern of this matrix.
+     */
+    BCRSMatrix& axpy(field_type alpha, const BCRSMatrix& b)
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if(N()!=b.N() || M() != b.M())
+        DUNE_THROW(RangeError, "Matrix sizes do not match!");
+#endif
+      RowIterator endi=end();
+      ConstRowIterator j=b.begin();
+      for(RowIterator i=begin(); i!=endi; ++i, ++j)
+        i->axpy(alpha, *j);
+
+      return *this;
+    }
+
     //===== linear maps
 
     //! y = A x