Skip to content
Snippets Groups Projects
Commit a0445bea authored by Elias Pipping's avatar Elias Pipping
Browse files

*Axy moved to dune-matrix-vector

parent ed3a7bbf
No related branches found
No related tags found
No related merge requests found
...@@ -14,6 +14,7 @@ ...@@ -14,6 +14,7 @@
#include <dune/istl/scaledidmatrix.hh> #include <dune/istl/scaledidmatrix.hh>
#include <dune/matrix-vector/axpy.hh> #include <dune/matrix-vector/axpy.hh>
#include <dune/matrix-vector/axy.hh>
#include <dune/matrix-vector/crossproduct.hh> #include <dune/matrix-vector/crossproduct.hh>
#include <dune/matrix-vector/transpose.hh> #include <dune/matrix-vector/transpose.hh>
...@@ -233,87 +234,6 @@ namespace Arithmetic ...@@ -233,87 +234,6 @@ namespace Arithmetic
Dune::MatrixVector::subtractProduct(a,b,c,d); 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$ //! Compute \f$(Ax,y)\f$
template <class OperatorType, class VectorType, class VectorType2> template <class OperatorType, class VectorType, class VectorType2>
typename VectorType::field_type typename VectorType::field_type
...@@ -321,8 +241,7 @@ namespace Arithmetic ...@@ -321,8 +241,7 @@ namespace Arithmetic
const VectorType &x, const VectorType &x,
const VectorType2 &y) const VectorType2 &y)
{ {
return OperatorHelper<OperatorType, return Dune::MatrixVector::Axy(A, x, y);
MatrixTraits<OperatorType>::isMatrix>::Axy(A, x, y);
} }
//! Compute \f$(b-Ax,y)\f$ //! Compute \f$(b-Ax,y)\f$
...@@ -333,8 +252,7 @@ namespace Arithmetic ...@@ -333,8 +252,7 @@ namespace Arithmetic
const VectorType &x, const VectorType &x,
const VectorType2 &y) const VectorType2 &y)
{ {
return OperatorHelper<OperatorType, return Dune::MatrixVector::bmAxy(A, b, x, y);
MatrixTraits<OperatorType>::isMatrix>::bmAxy(A, b, x, y);
} }
} }
......
...@@ -3,12 +3,12 @@ ...@@ -3,12 +3,12 @@
#ifndef DUNE_SOLVERS_COMPUTEENERGY_HH #ifndef DUNE_SOLVERS_COMPUTEENERGY_HH
#define DUNE_SOLVERS_COMPUTEENERGY_HH #define DUNE_SOLVERS_COMPUTEENERGY_HH
#include "common/arithmetic.hh" #include <dune/matrix-vector/axy.hh>
template <class M, class V> template <class M, class V>
double computeEnergy(const M& mat, const V& x) 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> template <class M, class V>
......
...@@ -5,8 +5,9 @@ ...@@ -5,8 +5,9 @@
#include <cmath> #include <cmath>
#include <dune/matrix-vector/axy.hh>
#include <dune/solvers/norms/norm.hh> #include <dune/solvers/norms/norm.hh>
#include <dune/solvers/common/arithmetic.hh>
#include <dune/solvers/iterationsteps/lineariterationstep.hh> #include <dune/solvers/iterationsteps/lineariterationstep.hh>
/** \brief Vector norm induced by linear operator /** \brief Vector norm induced by linear operator
...@@ -78,7 +79,7 @@ ...@@ -78,7 +79,7 @@
? *(iterationStep_->getMatrix()) ? *(iterationStep_->getMatrix())
: *matrix_; : *matrix_;
const field_type ret = Arithmetic::Axy(A, f, f); const field_type ret = Dune::MatrixVector::Axy(A, f, f);
if (ret < 0) if (ret < 0)
{ {
...@@ -95,7 +96,7 @@ ...@@ -95,7 +96,7 @@
const MatrixType& A, const MatrixType& A,
const double tol=1e-10) 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) if (ret < 0)
{ {
......
...@@ -5,7 +5,8 @@ ...@@ -5,7 +5,8 @@
#include <cmath> #include <cmath>
#include <dune/fufem/arithmetic.hh> #include <dune/matrix-vector/axy.hh>
#include <dune/solvers/solvers/iterativesolver.hh> #include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/iterationsteps/lineariterationstep.hh> #include <dune/solvers/iterationsteps/lineariterationstep.hh>
#include <dune/solvers/norms/norm.hh> #include <dune/solvers/norms/norm.hh>
...@@ -41,7 +42,7 @@ class TruncatedCGSolver : public IterativeSolver<VectorType> ...@@ -41,7 +42,7 @@ class TruncatedCGSolver : public IterativeSolver<VectorType>
const VectorType& b) const { const VectorType& b) const {
if (trustRegionNormMatrix_) if (trustRegionNormMatrix_)
return Arithmetic::Axy(*trustRegionNormMatrix_, b, a); return Dune::MatrixVector::Axy(*trustRegionNormMatrix_, b, a);
return a*b; return a*b;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment