Incorporate *Axy from Arithmetic

dune/matrix-vector/CMakeLists.txt

@@ -2,6 +2,7 @@
 install(FILES
   addtodiagonal.hh
   axpy.hh
+  axy.hh
   componentwisematrixmap.hh
   crossproduct.hh
   genericvectortools.hh

dune/matrix-vector/axy.hh

+#ifndef DUNE_MATRIX_VECTOR_AXY_HH
+#define DUNE_MATRIX_VECTOR_AXY_HH
+
+#include <cassert>
+
+#include "axpy.hh"
+#include "matrixtraits.hh"
+
+namespace Dune {
+namespace MatrixVector {
+  /** \brief Internal helper class for Matrix operations
+   *
+   */
+  template <class OperatorType, bool isMatrix>
+  struct OperatorHelper {
+    template <class VectorType, class VectorType2>
+    static typename VectorType::field_type Axy(const OperatorType &A,
+                                               const VectorType &x,
+                                               const VectorType2 &y) {
+      VectorType2 tmp = y;
+      tmp = 0;
+      addProduct(tmp, A, x);
+      return tmp * y;
+    }
+
+    template <class VectorType, class VectorType2>
+    static typename VectorType::field_type bmAxy(const OperatorType &A,
+                                                 const VectorType2 &b,
+                                                 const VectorType &x,
+                                                 const VectorType2 &y) {
+      VectorType2 tmp = b;
+      subtractProduct(tmp, A, x);
+      return tmp * y;
+    }
+  };
+
+  template <class MatrixType>
+  struct OperatorHelper<MatrixType, true> {
+    template <class VectorType, class VectorType2>
+    static typename VectorType::field_type Axy(const MatrixType &A,
+                                               const VectorType &x,
+                                               const VectorType2 &y) {
+      assert(x.N() == A.M());
+      assert(y.N() == A.N());
+
+      typename VectorType::field_type outer = 0.0;
+      typename VectorType2::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;
+    }
+
+    template <class VectorType, class VectorType2>
+    static typename VectorType::field_type bmAxy(const MatrixType &A,
+                                                 const VectorType2 &b,
+                                                 const VectorType &x,
+                                                 const VectorType2 &y) {
+      assert(x.N() == A.M());
+      assert(y.N() == A.N());
+      assert(y.N() == b.N());
+
+      typename VectorType::field_type outer = 0.0;
+      typename VectorType2::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$(Ax,y)\f$
+  template <class OperatorType, class VectorType, class VectorType2>
+  typename VectorType::field_type Axy(const OperatorType &A,
+                                      const VectorType &x,
+                                      const VectorType2 &y) {
+    return OperatorHelper<OperatorType,
+                          MatrixTraits<OperatorType>::isMatrix>::Axy(A, x, y);
+  }
+
+  //! Compute \f$(b-Ax,y)\f$
+  template <class OperatorType, class VectorType, class VectorType2>
+  typename VectorType::field_type bmAxy(const OperatorType &A,
+                                        const VectorType2 &b,
+                                        const VectorType &x,
+                                        const VectorType2 &y) {
+    return OperatorHelper<OperatorType,
+                          MatrixTraits<OperatorType>::isMatrix>::bmAxy(A, b, x,
+                                                                       y);
+  }
+}
+}
+
+#endif