diff --git a/dune/solvers/common/arithmetic.hh b/dune/solvers/common/arithmetic.hh
index b1d682be9c7314a3a3809c1f51d0e226670734fb..149324d4cc77b93f262e3cded659a2b01bd98234 100644
--- a/dune/solvers/common/arithmetic.hh
+++ b/dune/solvers/common/arithmetic.hh
@@ -14,6 +14,7 @@
#include <dune/istl/scaledidmatrix.hh>
#include <dune/matrix-vector/axpy.hh>
+#include <dune/matrix-vector/axy.hh>
#include <dune/matrix-vector/crossproduct.hh>
#include <dune/matrix-vector/transpose.hh>
@@ -233,87 +234,6 @@ namespace Arithmetic
Dune::MatrixVector::subtractProduct(a,b,c,d);
}
- /** \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;
- Dune::Solvers::resizeInitializeZero(tmp,y);
- 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
@@ -321,8 +241,7 @@ namespace Arithmetic
const VectorType &x,
const VectorType2 &y)
{
- return OperatorHelper<OperatorType,
- MatrixTraits<OperatorType>::isMatrix>::Axy(A, x, y);
+ return Dune::MatrixVector::Axy(A, x, y);
}
//! Compute \f$(b-Ax,y)\f$
@@ -333,8 +252,7 @@ namespace Arithmetic
const VectorType &x,
const VectorType2 &y)
{
- return OperatorHelper<OperatorType,
- MatrixTraits<OperatorType>::isMatrix>::bmAxy(A, b, x, y);
+ return Dune::MatrixVector::bmAxy(A, b, x, y);
}
}
diff --git a/dune/solvers/computeenergy.hh b/dune/solvers/computeenergy.hh
index e6933f046fac7dd787a9d8249700f7529287344e..cef779ceb235153fa00970e16a0a2a746a4fd479 100644
--- a/dune/solvers/computeenergy.hh
+++ b/dune/solvers/computeenergy.hh
@@ -3,12 +3,12 @@
#ifndef DUNE_SOLVERS_COMPUTEENERGY_HH
#define DUNE_SOLVERS_COMPUTEENERGY_HH
-#include "common/arithmetic.hh"
+#include <dune/matrix-vector/axy.hh>
template <class M, class V>
double computeEnergy(const M& mat, const V& x)
{
- return 0.5*Arithmetic::Axy(mat,x,x);
+ return 0.5*Dune::MatrixVector::Axy(mat,x,x);
}
template <class M, class V>
diff --git a/dune/solvers/norms/energynorm.hh b/dune/solvers/norms/energynorm.hh
index 2ed7b736169c9a880e470b8b6fac4a17c10c642b..f58a74adacf3dffdba7fabc2384557512388bdca 100644
--- a/dune/solvers/norms/energynorm.hh
+++ b/dune/solvers/norms/energynorm.hh
@@ -5,8 +5,9 @@
#include <cmath>
+#include <dune/matrix-vector/axy.hh>
+
#include <dune/solvers/norms/norm.hh>
-#include <dune/solvers/common/arithmetic.hh>
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
/** \brief Vector norm induced by linear operator
@@ -78,7 +79,7 @@
? *(iterationStep_->getMatrix())
: *matrix_;
- const field_type ret = Arithmetic::Axy(A, f, f);
+ const field_type ret = Dune::MatrixVector::Axy(A, f, f);
if (ret < 0)
{
@@ -95,7 +96,7 @@
const MatrixType& A,
const double tol=1e-10)
{
- const field_type ret = Arithmetic::Axy(A, u, u);
+ const field_type ret = Dune::MatrixVector::Axy(A, u, u);
if (ret < 0)
{
diff --git a/dune/solvers/solvers/tcgsolver.hh b/dune/solvers/solvers/tcgsolver.hh
index 942d786bc1e26417e55b596556ad9526956b4a08..703d5db78122417840fafffab0793f9a30265a44 100644
--- a/dune/solvers/solvers/tcgsolver.hh
+++ b/dune/solvers/solvers/tcgsolver.hh
@@ -5,7 +5,8 @@
#include <cmath>
-#include <dune/fufem/arithmetic.hh>
+#include <dune/matrix-vector/axy.hh>
+
#include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
#include <dune/solvers/norms/norm.hh>
@@ -41,7 +42,7 @@ class TruncatedCGSolver : public IterativeSolver<VectorType>
const VectorType& b) const {
if (trustRegionNormMatrix_)
- return Arithmetic::Axy(*trustRegionNormMatrix_, b, a);
+ return Dune::MatrixVector::Axy(*trustRegionNormMatrix_, b, a);
return a*b;