diff --git a/dune/matrix-vector/CMakeLists.txt b/dune/matrix-vector/CMakeLists.txt index 090d12d36b758e655f6127fbaaa269dc4d099b30..f654ddffaf282935624bec5649f08f9f18de5772 100644 --- a/dune/matrix-vector/CMakeLists.txt +++ b/dune/matrix-vector/CMakeLists.txt @@ -2,6 +2,7 @@ install(FILES addtodiagonal.hh axpy.hh + axy.hh componentwisematrixmap.hh crossproduct.hh genericvectortools.hh diff --git a/dune/matrix-vector/axy.hh b/dune/matrix-vector/axy.hh new file mode 100644 index 0000000000000000000000000000000000000000..c378c89b86c49b5a3f5e3d28cb356b52e673d42c --- /dev/null +++ b/dune/matrix-vector/axy.hh @@ -0,0 +1,108 @@ +#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