diff --git a/dune/solvers/common/arithmetic.hh b/dune/solvers/common/arithmetic.hh
index 1adca58ad9587cd103efe0d98112eacf26742f34..bee60bae9bf5ff2e1380d84b8ff8e1ef69faae35 100644
--- a/dune/solvers/common/arithmetic.hh
+++ b/dune/solvers/common/arithmetic.hh
@@ -13,6 +13,8 @@
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/scaledidmatrix.hh>
+#include <dune/matrix-vector/axpy.hh>
+
#include <dune/solvers/common/resize.hh>
/** \brief Namespace containing helper classes and functions for arithmetic operations
@@ -142,386 +144,6 @@ namespace Arithmetic
template <bool Condition, typename Type=void>
using enable_if_t = typename std::enable_if<Condition, Type>::type;
- /** \brief Internal helper class for product operations
- *
- */
- template<class A, class B, class C, bool AisScalar, bool BisScalar, bool CisScalar>
- struct ProductHelper
- {
- template <class ADummy=A, enable_if_t<!MatrixTraits<ADummy>::isMatrix,int> SFINAE_Dummy=0>
- static void addProduct(A& a, const B& b, const C& c)
- {
- b.umv(c, a);
- }
-
- template <class ADummy=A, enable_if_t<MatrixTraits<ADummy>::isMatrix,int> SFINAE_Dummy=0>
- static void addProduct(A& a, const B& b, const C& c)
- {
- typename B::ConstRowIterator bi = b.begin();
- typename B::ConstRowIterator bEnd = b.end();
- for(; bi!=bEnd; ++bi)
- {
- typename B::ConstColIterator bik = bi->begin();
- typename B::ConstColIterator biEnd= bi->end();
- for(; bik!=biEnd; ++bik)
- {
- typename C::ConstColIterator ckj =c[bik.index()].begin();
- typename C::ConstColIterator ckEnd=c[bik.index()].end();
- for(; ckj!=ckEnd; ++ckj)
- a[bi.index()][ckj.index()] += (*bik) * (*ckj);
- }
- }
- }
- };
-
- // Internal helper class for scaled product operations (i.e., b is always a scalar)
- template<class A, class Scalar, class B, class C, bool AisScalar, bool BisScalar, bool CisScalar>
- struct ScaledProductHelper
- {
- template <class ADummy=A, enable_if_t<!MatrixTraits<ADummy>::isMatrix,int> SFINAE_Dummy=0>
- static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c)
- {
- b.usmv(scalar, c, a);
- }
-
- template <class ADummy=A, enable_if_t<MatrixTraits<ADummy>::isMatrix,int> SFINAE_Dummy=0>
- static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c)
- {
- typename B::ConstRowIterator bi = b.begin();
- typename B::ConstRowIterator bEnd = b.end();
- for(; bi!=bEnd; ++bi)
- {
- typename B::ConstColIterator bik =b[bi.index()].begin();
- typename B::ConstColIterator biEnd=b[bi.index()].end();
- for(; bik!=biEnd; ++bik)
- {
- typename C::ConstColIterator ckj =c[bik.index()].begin();
- typename C::ConstColIterator ckEnd=c[bik.index()].end();
- for(; ckj!=ckEnd; ++ckj)
- a[bi.index()][ckj.index()] += scalar * (*bik) * (*ckj);
- }
- }
- }
- };
-
- template<class T, int n, int m, int k>
- struct ProductHelper<Dune::FieldMatrix<T,n,k>, Dune::FieldMatrix<T,n,m>, Dune::FieldMatrix<T,m,k>, false, false, false>
- {
- typedef Dune::FieldMatrix<T,n,k> A;
- typedef Dune::FieldMatrix<T,n,m> B;
- typedef Dune::FieldMatrix<T,m,k> C;
- static void addProduct(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>
- struct ScaledProductHelper<Dune::FieldMatrix<T,n,k>, T, Dune::FieldMatrix<T,n,m>, Dune::FieldMatrix<T,m,k>, false, false, false>
- {
- typedef Dune::FieldMatrix<T,n,k> A;
- typedef Dune::FieldMatrix<T,n,m> C;
- typedef Dune::FieldMatrix<T,m,k> D;
- static void addProduct(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>
- struct ProductHelper<Dune::DiagonalMatrix<T,n>, Dune::DiagonalMatrix<T,n>, Dune::DiagonalMatrix<T,n>, false, false, false>
- {
- typedef Dune::DiagonalMatrix<T,n> A;
- typedef A B;
- typedef B C;
- static void addProduct(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>
- struct ScaledProductHelper<Dune::DiagonalMatrix<T,n>, T, Dune::DiagonalMatrix<T,n>, Dune::DiagonalMatrix<T,n>, false, false, false>
- {
- typedef Dune::DiagonalMatrix<T,n> A;
- typedef A C;
- typedef C D;
- static void addProduct(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>
- struct ProductHelper<Dune::FieldMatrix<T,n,n>, Dune::DiagonalMatrix<T,n>, Dune::DiagonalMatrix<T,n>, false, false, false>
- {
- typedef Dune::FieldMatrix<T,n,n> A;
- typedef Dune::DiagonalMatrix<T,n> B;
- typedef B C;
- static void addProduct(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>
- struct ScaledProductHelper<Dune::FieldMatrix<T,n,n>, T, Dune::DiagonalMatrix<T,n>, Dune::DiagonalMatrix<T,n>, false, false, false>
- {
- typedef Dune::FieldMatrix<T,n,n> A;
- typedef Dune::DiagonalMatrix<T,n> C;
- typedef C D;
- static void addProduct(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);
- }
- };
-
- template<class T, int n, int m>
- struct ProductHelper<Dune::FieldMatrix<T,n,m>, Dune::DiagonalMatrix<T,n>, Dune::FieldMatrix<T,n,m>, false, false, false>
- {
- typedef Dune::FieldMatrix<T,n,m> A;
- typedef Dune::DiagonalMatrix<T,n> B;
- typedef A C;
- static void addProduct(A& a, const B& b, const C& c)
- {
- for (size_t i=0; i<n; ++i)
- a[i].axpy(b.diagonal(i),c[i]);
- }
- };
-
- template<class T, int n, int m>
- struct ScaledProductHelper<Dune::FieldMatrix<T,n,m>, T, Dune::DiagonalMatrix<T,n>, Dune::FieldMatrix<T,n,m>, false, false, false>
- {
- typedef Dune::FieldMatrix<T,n,m> A;
- typedef Dune::DiagonalMatrix<T,n> C;
- typedef A D;
- static void addProduct(A& a, const T& b, const C& c, const D& d)
- {
- for (size_t i=0; i<n; ++i)
- a[i].axpy(b*c.diagonal(i),d[i]);
- }
- };
-
- template<class T, int n, int m>
- struct ProductHelper<Dune::FieldMatrix<T,n,m>, Dune::FieldMatrix<T,n,m>, Dune::DiagonalMatrix<T,m>, false, false, false>
- {
- typedef Dune::FieldMatrix<T,n,m> A;
- typedef Dune::DiagonalMatrix<T,m> C;
- typedef A B;
-
- static void addProduct(A& a, const B& b, const C& c)
- {
- for (size_t i=0; i<n; ++i)
- for (size_t j=0; j<m; j++)
- a[i][j] += c.diagonal(j)*b[i][j];
- }
- };
-
- template<class T, int n, int m>
- struct ScaledProductHelper<Dune::FieldMatrix<T,n,m>, T, Dune::FieldMatrix<T,n,m>, Dune::DiagonalMatrix<T,m>, false, false, false>
- {
- typedef Dune::FieldMatrix<T,n,m> A;
- typedef Dune::DiagonalMatrix<T,m> D;
- typedef A C;
-
- static void addProduct(A& a, const T& b, const C& c, const D& d)
- {
- for (size_t i=0; i<n; ++i)
- for (size_t j=0; j<m; j++)
- a[i][j] += b*d.diagonal(j)*c[i][j];
- }
- };
-
- template<class T, int n>
- struct ProductHelper<Dune::ScaledIdentityMatrix<T,n>, Dune::ScaledIdentityMatrix<T,n>, Dune::ScaledIdentityMatrix<T,n>, false, false, false>
- {
- typedef Dune::ScaledIdentityMatrix<T,n> A;
- typedef A B;
- typedef B C;
-
- static void addProduct(A& a, const B& b, const C& c)
- {
- a.scalar() += b.scalar()*c.scalar();
- }
- };
-
- template<class T, class Scalar, int n>
- struct ScaledProductHelper<Dune::ScaledIdentityMatrix<T,n>, Scalar, Dune::ScaledIdentityMatrix<T,n>, Dune::ScaledIdentityMatrix<T,n>, false, false, false>
- {
- typedef Dune::ScaledIdentityMatrix<T,n> A;
- typedef A C;
- typedef C D;
-
- static void addProduct(A& a, const Scalar& b, const C& c, const D& d)
- {
- a.scalar() += b * c.scalar()*d.scalar();
- }
- };
-
- /** \brief Specialization for b being a scalar type
- */
- template<class A, class ScalarB, class C, bool AisScalar, bool CisScalar>
- struct ProductHelper<A, ScalarB, C, AisScalar, true, CisScalar>
- {
- typedef ScalarB B;
-
- template <class ADummy=A, enable_if_t<!MatrixTraits<ADummy>::isMatrix,int> SFINAE_Dummy=0>
- static void addProduct(A& a, const B& b, const C& c)
- {
- a.axpy(b, c);
- }
-
- template <class ADummy=A, enable_if_t<MatrixTraits<ADummy>::isMatrix,int> SFINAE_Dummy=0>
- static void addProduct(A& a, const B& b, const C& c)
- {
- typename C::ConstRowIterator ci = c.begin();
- typename C::ConstRowIterator cEnd = c.end();
- for(; ci!=cEnd; ++ci)
- {
- typename C::ConstColIterator cik =c[ci.index()].begin();
- typename C::ConstColIterator ciEnd=c[ci.index()].end();
- for(; cik!=ciEnd; ++cik)
- {
- a[ci.index()][cik.index()] += b*(*cik);
- }
- }
- }
- };
-
- template<class A, class Scalar, class ScalarB, class C, bool AisScalar, bool CisScalar>
- struct ScaledProductHelper<A, Scalar, ScalarB, C, AisScalar, true, CisScalar>
- {
- typedef ScalarB B;
-
- template <class ADummy=A, enable_if_t<!MatrixTraits<ADummy>::isMatrix,int> SFINAE_Dummy=0>
- static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c)
- {
- a.axpy(scalar*b, c);
- }
-
- template <class ADummy=A, enable_if_t<MatrixTraits<ADummy>::isMatrix,int> SFINAE_Dummy=0>
- static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c)
- {
- typename C::ConstRowIterator ci = c.begin();
- typename C::ConstRowIterator cEnd = c.end();
- for(; ci!=cEnd; ++ci)
- {
- typename C::ConstColIterator cik =c[ci.index()].begin();
- typename C::ConstColIterator ciEnd=c[ci.index()].end();
- for(; cik!=ciEnd; ++cik)
- {
- a[ci.index()][cik.index()] += scalar*b*(*cik);
- }
- }
- }
- };
-
- template<class ABlock, class ScalarB, class CBlock>
- struct ProductHelper<Dune::BCRSMatrix<ABlock>, ScalarB, Dune::BCRSMatrix<CBlock>, false, true, false>
- {
- typedef Dune::BCRSMatrix<ABlock> A;
- typedef ScalarB B;
- typedef Dune::BCRSMatrix<CBlock> C;
-
- static void addProduct(A& a, const B& b, const C& c)
- {
- // Workaround for the fact that Dune::BCRSMatrix assumes
- // its blocks to have an axpy() implementation, yet
- // Dune::DiagonalMatrix does not.
-
- //a.axpy(b, c);
-
- for (typename C::ConstIterator i=c.begin(); i!=c.end(); ++i) {
- const size_t iindex = i.index();
- for (typename C::row_type::ConstIterator j=i->begin(); j!=i->end(); ++j)
- ProductHelper<typename A::block_type, B, typename C::block_type, false, true, false>::addProduct(a[iindex][j.index()],b,*j);
- }
- }
- };
-
- template<class ABlock, class Scalar, class ScalarB, class CBlock>
- struct ScaledProductHelper<Dune::BCRSMatrix<ABlock>, Scalar, ScalarB, Dune::BCRSMatrix<CBlock>, false, true, false>
- {
- typedef Dune::BCRSMatrix<ABlock> A;
- typedef ScalarB B;
- typedef Dune::BCRSMatrix<CBlock> C;
-
- static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c)
- {
- // Workaround for the fact that Dune::BCRSMatrix assumes
- // its blocks to have an axpy() implementation, yet
- // Dune::DiagonalMatrix does not.
-
- //a.axpy(scalar*b, c);
-
- for (typename C::ConstIterator i=c.begin(); i!=c.end(); ++i) {
- const size_t iindex = i.index();
- for (typename C::row_type::ConstIterator j=i->begin(); j!=i->end(); ++j)
- ProductHelper<typename A::block_type, B, typename C::block_type, false, true, false>::addProduct(a[iindex][j.index()],scalar*b,*j);
- }
- }
- };
-
- template<class T, int n, class ScalarB>
- struct ProductHelper<Dune::ScaledIdentityMatrix<T,n>, ScalarB, Dune::ScaledIdentityMatrix<T,n>, false, true, false>
- {
- typedef Dune::ScaledIdentityMatrix<T,n> A;
- typedef ScalarB B;
- typedef Dune::ScaledIdentityMatrix<T,n> C;
-
- static void addProduct(A& a, const B& b, const C& c)
- {
- a.scalar() += b*c.scalar();
- }
- };
-
- template<class T, int n, class Scalar, class ScalarB>
- struct ScaledProductHelper<Dune::ScaledIdentityMatrix<T,n>, Scalar, ScalarB, Dune::ScaledIdentityMatrix<T,n>, false, true, false>
- {
- typedef Dune::ScaledIdentityMatrix<T,n> A;
- typedef ScalarB B;
- typedef Dune::ScaledIdentityMatrix<T,n> C;
-
- static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c)
- {
- a.scalar() += scalar*b*c.scalar();
- }
- };
-
- template<class A, class B, class C>
- struct ProductHelper<A, B, C, true, true, true>
- {
- static void addProduct(A& a, const B& b, const C& c)
- {
- a += b*c;
- }
- };
-
- template<class A, class B, class C, class D>
- struct ScaledProductHelper<A, B, C, D, true, true, true>
- {
- static void addProduct(A& a, const B& b, const C& c, const D& d)
- {
- a += b*c*d;
- }
- };
-
//! Helper class for computing the cross product
template <class T, int n>
struct CrossProductHelper
@@ -655,7 +277,7 @@ namespace Arithmetic
template<class A, class B, class C>
void addProduct(A& a, const B& b, const C& c)
{
- ProductHelper<A,B,C,ScalarTraits<A>::isScalar, ScalarTraits<B>::isScalar, ScalarTraits<C>::isScalar>::addProduct(a,b,c);
+ Dune::MatrixVector::addProduct(a,b,c);
}
/** \brief Subtract a product from some matrix or vector
@@ -672,7 +294,7 @@ namespace Arithmetic
template<class A, class B, class C>
void subtractProduct(A& a, const B& b, const C& c)
{
- ScaledProductHelper<A,int,B,C,ScalarTraits<A>::isScalar, ScalarTraits<B>::isScalar, ScalarTraits<C>::isScalar>::addProduct(a,-1,b,c);
+ Dune::MatrixVector::subtractProduct(a,b,c);
}
//! Compute the cross product of two vectors. Only works for n==3
@@ -704,7 +326,7 @@ namespace Arithmetic
typename std::enable_if<ScalarTraits<B>::isScalar, void>::type
addProduct(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>::addProduct(a,b,c,d);
+ Dune::MatrixVector::addProduct(a,b,c,d);
}
/** \brief Subtract a scaled product from some matrix or vector
@@ -722,7 +344,7 @@ namespace Arithmetic
typename std::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>::addProduct(a,-b,c,d);
+ Dune::MatrixVector::subtractProduct(a,b,c,d);
}
/** \brief Internal helper class for Matrix operations
diff --git a/dune/solvers/common/staticmatrixtools.hh b/dune/solvers/common/staticmatrixtools.hh
index c9dfabe0679de3f370020b670bb89dc310e1b06d..522af0eb76eb9c568ba82c7668165d9f3dfbd2a5 100644
--- a/dune/solvers/common/staticmatrixtools.hh
+++ b/dune/solvers/common/staticmatrixtools.hh
@@ -9,6 +9,9 @@
#include "dune/istl/bcrsmatrix.hh"
#include "dune/istl/matrixindexset.hh"
+#include <dune/matrix-vector/axpy.hh>
+#include <dune/matrix-vector/scalartraits.hh>
+
#include "arithmetic.hh"
namespace StaticMatrix
@@ -169,7 +172,7 @@ namespace StaticMatrix
{
static void addTransformedMatrix(MatrixA& A, const ScalarTransform& T1, const MatrixB& B, const ScalarTransform& T2)
{
- Arithmetic::addProduct(A, T1*T2, B);
+ Dune::MatrixVector::addProduct(A, T1*T2, B);
}
};
@@ -283,7 +286,7 @@ namespace StaticMatrix
void addTransformedMatrix(MatrixA& A, const TransformationMatrix& T1, const MatrixB& B, const TransformationMatrix& T2)
{
TransformMatrixHelper<MatrixA,MatrixB,TransformationMatrix,
- Arithmetic::ScalarTraits<MatrixA>::isScalar, Arithmetic::ScalarTraits<MatrixB>::isScalar, Arithmetic::ScalarTraits<TransformationMatrix>::isScalar>
+ Dune::MatrixVector::ScalarTraits<MatrixA>::isScalar, Dune::MatrixVector::ScalarTraits<MatrixB>::isScalar, Dune::MatrixVector::ScalarTraits<TransformationMatrix>::isScalar>
::addTransformedMatrix(A,T1,B,T2);
}
@@ -292,7 +295,7 @@ namespace StaticMatrix
{
A = 0;
TransformMatrixHelper<MatrixA,MatrixB,TransformationMatrix,
- Arithmetic::ScalarTraits<MatrixA>::isScalar, Arithmetic::ScalarTraits<MatrixB>::isScalar, Arithmetic::ScalarTraits<TransformationMatrix>::isScalar>
+ Dune::MatrixVector::ScalarTraits<MatrixA>::isScalar, Dune::MatrixVector::ScalarTraits<MatrixB>::isScalar, Dune::MatrixVector::ScalarTraits<TransformationMatrix>::isScalar>
::addTransformedMatrix(A,T1,B,T2);
}
diff --git a/dune/solvers/iterationsteps/blockgsstep.cc b/dune/solvers/iterationsteps/blockgsstep.cc
index d49c588d7335a68459b04c0cfb2b8ceeb54867ba..9092e4d07649e47491ec8386878b3a7342df6f57 100644
--- a/dune/solvers/iterationsteps/blockgsstep.cc
+++ b/dune/solvers/iterationsteps/blockgsstep.cc
@@ -3,7 +3,7 @@
#include <cassert>
-#include <dune/solvers/common/arithmetic.hh>
+#include <dune/matrix-vector/axpy.hh>
template<class MatrixType, class DiscFuncType, class BitVectorType>
inline
@@ -15,7 +15,7 @@ residual(int index, VectorBlock& r) const
const auto& row = (*this->mat_)[index];
for (auto cIt = row.begin(); cIt!=row.end(); ++cIt) {
// r_i -= A_ij x_j
- Arithmetic::subtractProduct(r, *cIt, (*this->x_)[cIt.index()]);
+ Dune::MatrixVector::subtractProduct(r, *cIt, (*this->x_)[cIt.index()]);
}
}
diff --git a/dune/solvers/iterationsteps/linegsstep.hh b/dune/solvers/iterationsteps/linegsstep.hh
index 45c0fc36c8fed1f23e0257524d63fbae63cbd89a..1c59163995009e35d3b2b9703d3108de4eca76f7 100755
--- a/dune/solvers/iterationsteps/linegsstep.hh
+++ b/dune/solvers/iterationsteps/linegsstep.hh
@@ -3,9 +3,9 @@
#ifndef DUNE_LINE_GAUSS_SEIDEL_STEP_HH
#define DUNE_LINE_GAUSS_SEIDEL_STEP_HH
-#include <dune/common/bitsetvector.hh>
+#include <dune/matrix-vector/axpy.hh>
-#include <dune/solvers/common/arithmetic.hh>
+#include <dune/common/bitsetvector.hh>
template<class MatrixType,
class DiscFuncType,
@@ -51,7 +51,7 @@ template<class MatrixType,
const auto& row = (*this->mat_)[index];
for (auto cIt = row.begin(); cIt!=row.end(); ++cIt) {
// r_i -= A_ij x_j
- Arithmetic::subtractProduct(r, *cIt, (*this->x_)[cIt.index()]);
+ Dune::MatrixVector::subtractProduct(r, *cIt, (*this->x_)[cIt.index()]);
}
}
};
diff --git a/dune/solvers/iterationsteps/projectedblockgsstep.hh b/dune/solvers/iterationsteps/projectedblockgsstep.hh
index 111b093078a88d1caf3644e8507838cd1ce7b69a..9a491168828f8b27f4fbd4e98b602d4865f7d881 100644
--- a/dune/solvers/iterationsteps/projectedblockgsstep.hh
+++ b/dune/solvers/iterationsteps/projectedblockgsstep.hh
@@ -5,9 +5,9 @@
#include <vector>
-#include <dune/solvers/common/boxconstraint.hh>
+#include <dune/matrix-vector/axpy.hh>
-#include <dune/solvers/common/arithmetic.hh>
+#include <dune/solvers/common/boxconstraint.hh>
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
template<class MatrixType, class VectorType>
@@ -42,7 +42,7 @@ public:
const auto& row = (*this->mat_)[index];
for (auto cIt = row.begin(); cIt!=row.end(); ++cIt) {
// r_i -= A_ij x_j
- Arithmetic::subtractProduct(r, *cIt, (*this->x_)[cIt.index()]);
+ Dune::MatrixVector::subtractProduct(r, *cIt, (*this->x_)[cIt.index()]);
}
}
diff --git a/dune/solvers/iterationsteps/truncatedblockgsstep.hh b/dune/solvers/iterationsteps/truncatedblockgsstep.hh
index c0c2083810079cf36e42e1aa97508c4be03e28dc..d84e285c5f8edbc0b15a1bad242fe9040329be5b 100644
--- a/dune/solvers/iterationsteps/truncatedblockgsstep.hh
+++ b/dune/solvers/iterationsteps/truncatedblockgsstep.hh
@@ -7,12 +7,13 @@
#include <dune/common/bitsetvector.hh>
+#include <dune/matrix-vector/axpy.hh>
+
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
#include <dune/solvers/operators/sumoperator.hh>
#include <dune/solvers/common/staticmatrixtools.hh>
#include <dune/solvers/common/arithmetic.hh>
-
namespace TruncatedBlockGSStepNS {
template<int blocklevel> struct RecursiveGSStep;
}
@@ -186,7 +187,7 @@ public:
{
size_t col = it.index();
if (col == row)
- Arithmetic::addProduct(Aii, 1.0, *it);
+ Dune::MatrixVector::addProduct(Aii, 1.0, *it);
else
it->mmv(x[col],r);
}
@@ -198,7 +199,7 @@ public:
const LRBlock& mki(lowrank_mat.lowRankFactor()[k][row]);
// Aii += m^tm[i][i]
// this should actually be \sum(lr_block^t*lr_block) if the blocks are nonsymmetric
- Arithmetic::addProduct(Aii,mki,mki);
+ Dune::MatrixVector::addProduct(Aii,mki,mki);
// r -= m^tm*x, where x[row] were zero
mki.mmv(mx[k], r);
diff --git a/dune/solvers/transferoperators/compressedmultigridtransfer.hh b/dune/solvers/transferoperators/compressedmultigridtransfer.hh
index 31a32163d8ca9624816703ccc42fa5c9a66f62b2..1cdd9ebd8891da21db9062fb4ec3fc30310dde3c 100644
--- a/dune/solvers/transferoperators/compressedmultigridtransfer.hh
+++ b/dune/solvers/transferoperators/compressedmultigridtransfer.hh
@@ -10,8 +10,8 @@
#include <dune/common/bitsetvector.hh>
#include <dune/common/shared_ptr.hh>
+#include <dune/matrix-vector/axpy.hh>
#include <dune/solvers/operators/sumoperator.hh>
-#include <dune/solvers/common/staticmatrixtools.hh>
#include "genericmultigridtransfer.hh"
#include "multigridtransfer.hh"
@@ -265,7 +265,7 @@ public:
size_t j = j_it.index();
for (size_t k=0; k<lr_fine_factor.N(); ++k)
- Arithmetic::addProduct(lr_coarse_factor[k][j], *j_it, lr_fine_factor[k][i]);
+ Dune::MatrixVector::addProduct(lr_coarse_factor[k][j], *j_it, lr_fine_factor[k][i]);
}
}
}
diff --git a/dune/solvers/transferoperators/genericmultigridtransfer.hh b/dune/solvers/transferoperators/genericmultigridtransfer.hh
index b1c3ffac3b583efcae2a26a0992caf07944d1bd2..b4a6ca4284c98a5bf5fb8af7ea7c3c5dbd34a00c 100644
--- a/dune/solvers/transferoperators/genericmultigridtransfer.hh
+++ b/dune/solvers/transferoperators/genericmultigridtransfer.hh
@@ -12,8 +12,9 @@
#include <dune/localfunctions/lagrange/pqkfactory.hh>
+#include <dune/matrix-vector/axpy.hh>
+
#include "dune/solvers/common/staticmatrixtools.hh"
-#include <dune/solvers/common/arithmetic.hh>
/** \brief Restriction and prolongation operator for standard multigrid
@@ -697,7 +698,7 @@ public:
// Compute cm = im^T * m * jm
if(TransferOperatorType::block_type::rows==1)
- Arithmetic::addProduct(cm, (*im)[0][0] * (*jm)[0][0], *m);
+ Dune::MatrixVector::addProduct(cm, (*im)[0][0] * (*jm)[0][0], *m);
else
StaticMatrix::addTransformedMatrix(cm, *im, *m, *jm);
}
@@ -764,7 +765,7 @@ public:
// Compute cm = im^T * m * jm
if(TransferOperatorType::block_type::rows==1)
- Arithmetic::addProduct(cm, (*im)[0][0] * (*jm)[0][0], *m);
+ Dune::MatrixVector::addProduct(cm, (*im)[0][0] * (*jm)[0][0], *m);
else
StaticMatrix::addTransformedMatrix(cm, *im, *m, *jm);
}