diff --git a/dune/istl/allocator.hh b/dune/istl/allocator.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c6f607fb41a56f7d0e871316ef0f15c7cc121540
--- /dev/null
+++ b/dune/istl/allocator.hh
@@ -0,0 +1,37 @@
+#ifndef DUNE_ISTL_ALLOCATOR_HH
+#define DUNE_ISTL_ALLOCATOR_HH
+
+#include <dune/common/typetraits.hh>
+#include <memory>
+
+namespace Dune {
+
+ template<typename T>
+ struct exists{
+ static const bool value = true;
+ };
+
+ template<typename T, typename = void>
+ struct DefaultAllocatorTraits
+ {
+ using type = std::allocator<T>;
+ };
+
+ template<typename T>
+ struct DefaultAllocatorTraits<T, void_t<typename T::allocator_type> >
+ {
+ using type = typename T::allocator_type;
+ };
+
+ template<typename T>
+ struct AllocatorTraits : public DefaultAllocatorTraits<T> {};
+
+ template<typename T>
+ using AllocatorType = typename AllocatorTraits<T>::type;
+
+ template<typename T, typename X>
+ using ReboundAllocatorType = typename AllocatorTraits<T>::type::template rebind<X>::other;
+
+} // end namespace Dune
+
+#endif // DUNE_ISTL_ALLOCATOR_HH
diff --git a/dune/istl/ilu.hh b/dune/istl/ilu.hh
index edf2380744031fc071477fc9984c97ef8652f7d6..4b3f19aff5a962b6ba92ac510950055bc4618e95 100644
--- a/dune/istl/ilu.hh
+++ b/dune/istl/ilu.hh
@@ -190,7 +190,8 @@ namespace Dune {
{
// we misuse the storage to store an int. If the field_type is std::complex, we have to access the real part
// starting from C++11, we can use std::real to always return a real value, even if it is double/float
- int generation = (int) std::real( firstmatrixelement(*kj) );
+ using std::real;
+ int generation = (int) real( firstmatrixelement(*kj) );
if (generation<n)
{
mapiterator ij = rowpattern.find(kj.index());
diff --git a/dune/istl/paamg/aggregates.hh b/dune/istl/paamg/aggregates.hh
index 850355e5b49e76fdd180018a4810d63cbe960fde..65340d741d8dc7caf217725083560dad5a0c37d9 100644
--- a/dune/istl/paamg/aggregates.hh
+++ b/dune/istl/paamg/aggregates.hh
@@ -732,17 +732,9 @@ namespace Dune
}
private:
/** @brief Prevent copying. */
- AggregatesMap(const AggregatesMap<V>& map)
- {
- throw "Auch!";
- }
-
+ AggregatesMap(const AggregatesMap<V>&) = delete;
/** @brief Prevent assingment. */
- AggregatesMap<V>& operator=(const AggregatesMap<V>& map)
- {
- throw "Auch!";
- return this;
- }
+ AggregatesMap<V>& operator=(const AggregatesMap<V>&) = delete;
/**
* @brief The aggregates the vertices belong to.
diff --git a/dune/istl/paamg/amg.hh b/dune/istl/paamg/amg.hh
index 47b96a99d93bcb3beeb8b8b79b64b9f8ba55f9dc..47711a7ef49ff9f972397846b6f892d9f17a979f 100644
--- a/dune/istl/paamg/amg.hh
+++ b/dune/istl/paamg/amg.hh
@@ -512,15 +512,34 @@ namespace Dune
{
if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
// We are still participating on this level
+
+ // we have to allocate these types using the rebound allocator
+ // in order to ensure that we fulfill the alignement requirements
solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
*scalarProduct_,
*coarseSmoother_, 1E-2, 1000, 0));
else
solver_.reset();
}else
- solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
- *scalarProduct_,
- *coarseSmoother_, 1E-2, 1000, 0));
+ {
+ solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
+ *scalarProduct_,
+ *coarseSmoother_, 1E-2, 1000, 0));
+ // // we have to allocate these types using the rebound allocator
+ // // in order to ensure that we fulfill the alignement requirements
+ // using Alloc = typename A::template rebind<BiCGSTABSolver<X>>::other;
+ // Alloc alloc;
+ // auto p = alloc.allocate(1);
+ // alloc.construct(p,
+ // const_cast<M&>(*matrices_->matrices().coarsest()),
+ // *scalarProduct_,
+ // *coarseSmoother_, 1E-2, 1000, 0);
+ // solver_.reset(p,[](BiCGSTABSolver<X>* p){
+ // Alloc alloc;
+ // alloc.destroy(p);
+ // alloc.deallocate(p,1);
+ // });
+ }
}
}
diff --git a/dune/istl/paamg/hierarchy.hh b/dune/istl/paamg/hierarchy.hh
index f186c713093bdb1c9a3d650becc243a1b65146c6..f3dec8cfacf77ad81726c4c3783a97c5ee3b5f2d 100644
--- a/dune/istl/paamg/hierarchy.hh
+++ b/dune/istl/paamg/hierarchy.hh
@@ -96,13 +96,6 @@ namespace Dune
MemberType* redistributed_;
};
public:
- // enum{
- // /**
- // * @brief If true only the method addCoarser will be usable
- // * otherwise only the method addFiner will be usable.
- // */
- // coarsen = b
- // };
/**
* @brief The allocator to use for the list elements.
@@ -385,8 +378,8 @@ namespace Dune
* @brief Coarsen the vector hierarchy according to the matrix hierarchy.
* @param hierarchy The vector hierarchy to coarsen.
*/
- template<class V, class TA>
- void coarsenVector(Hierarchy<BlockVector<V,TA> >& hierarchy) const;
+ template<class V, class BA, class TA>
+ void coarsenVector(Hierarchy<BlockVector<V,BA>, TA>& hierarchy) const;
/**
* @brief Coarsen the smoother hierarchy according to the matrix hierarchy.
@@ -442,8 +435,7 @@ namespace Dune
*/
const RedistributeInfoList& redistributeInformation() const;
-
- typename MatrixOperator::field_type getProlongationDampingFactor() const
+ double getProlongationDampingFactor() const
{
return prolongDamp_;
}
@@ -478,7 +470,7 @@ namespace Dune
/** @brief The maximum number of level across all processors.*/
int maxlevels_;
- typename MatrixOperator::field_type prolongDamp_;
+ double prolongDamp_;
/**
* @brief functor to print matrix statistics.
@@ -1114,8 +1106,8 @@ namespace Dune
}
template<class M, class IS, class A>
- template<class V, class TA>
- void MatrixHierarchy<M,IS,A>::coarsenVector(Hierarchy<BlockVector<V,TA> >& hierarchy) const
+ template<class V, class BA, class TA>
+ void MatrixHierarchy<M,IS,A>::coarsenVector(Hierarchy<BlockVector<V,BA>, TA>& hierarchy) const
{
assert(hierarchy.levels()==1);
typedef typename ParallelMatrixHierarchy::ConstIterator Iterator;
diff --git a/dune/istl/preconditioners.hh b/dune/istl/preconditioners.hh
index 2728fe3c565880772f3b1b54c830397d6c79b690..84acf024cc96bca563d31be12f2b23e0c224c44e 100644
--- a/dune/istl/preconditioners.hh
+++ b/dune/istl/preconditioners.hh
@@ -74,6 +74,9 @@ namespace Dune {
typedef typename O::range_type range_type;
//! \brief The field type of the preconditioner.
typedef typename range_type::field_type field_type;
+ //! \brief scalar type underlying the field_type
+ typedef SimdScalar<field_type> scalar_field_type;
+ //! \brief type of the wrapped inverse operator
typedef O InverseOperator;
/**
@@ -137,6 +140,8 @@ namespace Dune {
typedef Y range_type;
//! \brief The field type of the preconditioner.
typedef typename X::field_type field_type;
+ //! \brief scalar type underlying the field_type
+ typedef SimdScalar<field_type> scalar_field_type;
/*! \brief Constructor.
@@ -145,7 +150,7 @@ namespace Dune {
\param n The number of iterations to perform.
\param w The relaxation factor.
*/
- SeqSSOR (const M& A, int n, field_type w)
+ SeqSSOR (const M& A, int n, scalar_field_type w)
: _A_(A), _n(n), _w(w)
{
CheckIfDiagonalPresent<M,l>::check(_A_);
@@ -198,7 +203,7 @@ namespace Dune {
//! \brief The number of steps to do in apply
int _n;
//! \brief The relaxation factor to use
- field_type _w;
+ scalar_field_type _w;
};
@@ -225,6 +230,8 @@ namespace Dune {
typedef Y range_type;
//! \brief The field type of the preconditioner.
typedef typename X::field_type field_type;
+ //! \brief scalar type underlying the field_type
+ typedef SimdScalar<field_type> scalar_field_type;
/*! \brief Constructor.
@@ -233,7 +240,7 @@ namespace Dune {
\param n The number of iterations to perform.
\param w The relaxation factor.
*/
- SeqSOR (const M& A, int n, field_type w)
+ SeqSOR (const M& A, int n, scalar_field_type w)
: _A_(A), _n(n), _w(w)
{
CheckIfDiagonalPresent<M,l>::check(_A_);
@@ -303,7 +310,7 @@ namespace Dune {
//! \brief The number of steps to perform in apply.
int _n;
//! \brief The relaxation factor to use.
- field_type _w;
+ scalar_field_type _w;
};
@@ -328,6 +335,8 @@ namespace Dune {
typedef Y range_type;
//! \brief The field type of the preconditioner.
typedef typename X::field_type field_type;
+ //! \brief scalar type underlying the field_type
+ typedef SimdScalar<field_type> scalar_field_type;
/*! \brief Constructor.
@@ -336,7 +345,7 @@ namespace Dune {
\param n The number of iterations to perform.
\param w The relaxation factor.
*/
- SeqGS (const M& A, int n, field_type w)
+ SeqGS (const M& A, int n, scalar_field_type w)
: _A_(A), _n(n), _w(w)
{
CheckIfDiagonalPresent<M,l>::check(_A_);
@@ -387,7 +396,7 @@ namespace Dune {
//! \brief The number of iterations to perform in apply.
int _n;
//! \brief The relaxation factor to use.
- field_type _w;
+ scalar_field_type _w;
};
@@ -412,6 +421,8 @@ namespace Dune {
typedef Y range_type;
//! \brief The field type of the preconditioner.
typedef typename X::field_type field_type;
+ //! \brief scalar type underlying the field_type
+ typedef SimdScalar<field_type> scalar_field_type;
/*! \brief Constructor.
@@ -420,7 +431,7 @@ namespace Dune {
\param n The number of iterations to perform.
\param w The relaxation factor.
*/
- SeqJac (const M& A, int n, field_type w)
+ SeqJac (const M& A, int n, scalar_field_type w)
: _A_(A), _n(n), _w(w)
{
CheckIfDiagonalPresent<M,l>::check(_A_);
@@ -471,7 +482,7 @@ namespace Dune {
//! \brief The number of steps to perform during apply.
int _n;
//! \brief The relaxation parameter to use.
- field_type _w;
+ scalar_field_type _w;
};
@@ -498,6 +509,8 @@ namespace Dune {
typedef Y range_type;
//! \brief The field type of the preconditioner.
typedef typename X::field_type field_type;
+ //! \brief scalar type underlying the field_type
+ typedef SimdScalar<field_type> scalar_field_type;
/*! \brief Constructor.
@@ -505,10 +518,11 @@ namespace Dune {
\param A The matrix to operate on.
\param w The relaxation factor.
*/
- SeqILU0 (const M& A, field_type w)
- : ILU(A) // copy A
+ SeqILU0 (const M& A, scalar_field_type w)
+ : _w(w),
+ ILU(A) // copy A
+
{
- _w =w;
bilu0_decomposition(ILU);
}
@@ -552,7 +566,7 @@ namespace Dune {
private:
//! \brief The relaxation factor to use.
- field_type _w;
+ scalar_field_type _w;
//! \brief The ILU0 decomposition of the matrix.
matrix_type ILU;
};
@@ -582,6 +596,8 @@ namespace Dune {
typedef Y range_type;
//! \brief The field type of the preconditioner.
typedef typename X::field_type field_type;
+ //! \brief scalar type underlying the field_type
+ typedef SimdScalar<field_type> scalar_field_type;
/*! \brief Constructor.
@@ -590,11 +606,11 @@ namespace Dune {
\param n The number of iterations to perform.
\param w The relaxation factor.
*/
- SeqILUn (const M& A, int n, field_type w)
- : ILU(A.N(),A.M(),M::row_wise)
+ SeqILUn (const M& A, int n, scalar_field_type w)
+ : ILU(A.N(),A.M(),M::row_wise),
+ _n(n),
+ _w(w)
{
- _n = n;
- _w = w;
bilu_decomposition(A,n,ILU);
}
@@ -642,7 +658,7 @@ namespace Dune {
//! \brief The number of steps to perform in apply.
int _n;
//! \brief The relaxation factor to use.
- field_type _w;
+ scalar_field_type _w;
};
@@ -664,16 +680,17 @@ namespace Dune {
typedef Y range_type;
//! \brief The field type of the preconditioner.
typedef typename X::field_type field_type;
+ //! \brief scalar type underlying the field_type
+ typedef SimdScalar<field_type> scalar_field_type;
/*! \brief Constructor.
Constructor gets all parameters to operate the prec.
\param w The relaxation factor.
*/
- Richardson (field_type w=1.0)
- {
- _w = w;
- }
+ Richardson (scalar_field_type w=1.0) :
+ _w(w)
+ {}
/*!
\brief Prepare the preconditioner.
@@ -715,7 +732,7 @@ namespace Dune {
private:
//! \brief The relaxation factor to use.
- field_type _w;
+ scalar_field_type _w;
};
diff --git a/dune/istl/solver.hh b/dune/istl/solver.hh
index 672710c6a0edc8635992b77516ccaa39423eb360..774b2f5ae9a1c95be6b5f0e1cc48318b126bc761 100644
--- a/dune/istl/solver.hh
+++ b/dune/istl/solver.hh
@@ -9,6 +9,7 @@
#include <dune/common/exceptions.hh>
#include <dune/common/shared_ptr.hh>
+#include <dune/common/simd.hh>
#include "solvertype.hh"
#include "preconditioner.hh"
@@ -98,6 +99,9 @@ namespace Dune
//! \brief The real type of the field type (is the same if using real numbers, but differs for std::complex)
typedef typename FieldTraits<field_type>::real_type real_type;
+ //! \brief scalar type underlying the field_type
+ typedef SimdScalar<real_type> scalar_real_type;
+
/**
\brief Apply inverse operator,
@@ -191,6 +195,7 @@ namespace Dune
using typename InverseOperator<X,Y>::range_type;
using typename InverseOperator<X,Y>::field_type;
using typename InverseOperator<X,Y>::real_type;
+ using typename InverseOperator<X,Y>::scalar_real_type;
/*!
\brief General constructor to initialize an iterative solver.
@@ -211,7 +216,7 @@ namespace Dune
<li> 2 : print line for each iteration </li>
</ul>
*/
- IterativeSolver (LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec, real_type reduction, int maxit, int verbose) :
+ IterativeSolver (LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec, scalar_real_type reduction, int maxit, int verbose) :
_op(stackobject_to_shared_ptr(op)),
_prec(stackobject_to_shared_ptr(prec)),
_sp(new SeqScalarProduct<X>),
@@ -244,7 +249,7 @@ namespace Dune
</ul>
*/
IterativeSolver (LinearOperator<X,Y>& op, ScalarProduct<X>& sp, Preconditioner<X,Y>& prec,
- real_type reduction, int maxit, int verbose) :
+ scalar_real_type reduction, int maxit, int verbose) :
_op(stackobject_to_shared_ptr(op)),
_prec(stackobject_to_shared_ptr(prec)),
_sp(stackobject_to_shared_ptr(sp)),
@@ -270,7 +275,7 @@ namespace Dune
*/
virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
{
- real_type saved_reduction = _reduction;
+ scalar_real_type saved_reduction = _reduction;
_reduction = reduction;
static_cast<InverseOperator<X,Y>*>(this)->apply(x,b,res);
_reduction = saved_reduction;
@@ -286,7 +291,7 @@ namespace Dune
std::shared_ptr<LinearOperator<X,Y>> _op;
std::shared_ptr<Preconditioner<X,Y>> _prec;
std::shared_ptr<ScalarProduct<X>> _sp;
- real_type _reduction;
+ scalar_real_type _reduction;
int _maxit;
int _verbose;
SolverCategory::Category _category;
diff --git a/dune/istl/solvers.hh b/dune/istl/solvers.hh
index ebc0c6a0499ed409d3e4f74c1550659855ff4a7b..41e76e310faf0346012a16acf65ea8b92c78e30d 100644
--- a/dune/istl/solvers.hh
+++ b/dune/istl/solvers.hh
@@ -14,6 +14,7 @@
#include <array>
#include <type_traits>
+#include "allocator.hh"
#include "istlexception.hh"
#include "operators.hh"
#include "scalarproducts.hh"
@@ -903,13 +904,23 @@ namespace Dune {
using typename IterativeSolver<X,Y>::field_type;
using typename IterativeSolver<X,Y>::real_type;
+ private:
+ using typename IterativeSolver<X,X>::scalar_real_type;
+
+ //! \bief field_type Allocator retrieved from domain type
+ using fAlloc = ReboundAllocatorType<X,field_type>;
+ //! \bief real_type Allocator retrieved from domain type
+ using rAlloc = ReboundAllocatorType<X,real_type>;
+
+ public:
+
/*!
\brief Set up RestartedGMResSolver solver.
\copydoc LoopSolver::LoopSolver(L&,P&,double,int,int)
\param restart number of GMRes cycles before restart
*/
- RestartedGMResSolver (LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec, real_type reduction, int restart, int maxit, int verbose) :
+ RestartedGMResSolver (LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec, scalar_real_type reduction, int restart, int maxit, int verbose) :
IterativeSolver<X,Y>::IterativeSolver(op,prec,reduction,maxit,verbose),
_restart(restart)
{}
@@ -920,7 +931,7 @@ namespace Dune {
\copydoc LoopSolver::LoopSolver(L&,S&,P&,double,int,int)
\param restart number of GMRes cycles before restart
*/
- RestartedGMResSolver (LinearOperator<X,Y>& op, ScalarProduct<X>& sp, Preconditioner<X,Y>& prec, real_type reduction, int restart, int maxit, int verbose) :
+ RestartedGMResSolver (LinearOperator<X,Y>& op, ScalarProduct<X>& sp, Preconditioner<X,Y>& prec, scalar_real_type reduction, int restart, int maxit, int verbose) :
IterativeSolver<X,Y>::IterativeSolver(op,sp,prec,reduction,maxit,verbose),
_restart(restart)
{}
@@ -953,13 +964,13 @@ namespace Dune {
const int m = _restart;
real_type norm, norm_old = 0.0, norm_0;
int j = 1;
- std::vector<field_type> s(m+1), sn(m);
- std::vector<real_type> cs(m);
+ std::vector<field_type,fAlloc> s(m+1), sn(m);
+ std::vector<real_type,rAlloc> cs(m);
// need copy of rhs if GMRes has to be restarted
Y b2(b);
// helper vector
Y w(b);
- std::vector< std::vector<field_type> > H(m+1,s);
+ std::vector< std::vector<field_type,fAlloc> > H(m+1,s);
std::vector<F> v(m+1,b);
// start timer
@@ -1105,11 +1116,11 @@ namespace Dune {
}
void update(X& w, int i,
- const std::vector<std::vector<field_type> >& H,
- const std::vector<field_type>& s,
+ const std::vector<std::vector<field_type,fAlloc> >& H,
+ const std::vector<field_type,fAlloc>& s,
const std::vector<X>& v) {
// solution vector of the upper triangular system
- std::vector<field_type> y(s);
+ std::vector<field_type,fAlloc> y(s);
// backsolve
for(int a=i-1; a>=0; a--) {
@@ -1219,13 +1230,21 @@ namespace Dune {
using typename IterativeSolver<X,X>::field_type;
using typename IterativeSolver<X,X>::real_type;
+ private:
+ using typename IterativeSolver<X,X>::scalar_real_type;
+
+ //! \bief field_type Allocator retrieved from domain type
+ using fAlloc = ReboundAllocatorType<X,field_type>;
+
+ public:
+
/*!
\brief Set up nonlinear preconditioned conjugate gradient solver.
\copydoc LoopSolver::LoopSolver(L&,P&,double,int,int)
\param restart number of GMRes cycles before restart
*/
- GeneralizedPCGSolver (LinearOperator<X,X>& op, Preconditioner<X,X>& prec, real_type reduction, int maxit, int verbose, int restart = 10) :
+ GeneralizedPCGSolver (LinearOperator<X,X>& op, Preconditioner<X,X>& prec, scalar_real_type reduction, int maxit, int verbose, int restart = 10) :
IterativeSolver<X,X>::IterativeSolver(op,prec,reduction,maxit,verbose),
_restart(restart)
{}
@@ -1237,7 +1256,7 @@ namespace Dune {
\param restart When to restart the construction of
the Krylov search space.
*/
- GeneralizedPCGSolver (LinearOperator<X,X>& op, ScalarProduct<X>& sp, Preconditioner<X,X>& prec, real_type reduction, int maxit, int verbose, int restart = 10) :
+ GeneralizedPCGSolver (LinearOperator<X,X>& op, ScalarProduct<X>& sp, Preconditioner<X,X>& prec, scalar_real_type reduction, int maxit, int verbose, int restart = 10) :
IterativeSolver<X,X>::IterativeSolver(op,sp,prec,reduction,maxit,verbose),
_restart(restart)
{}
@@ -1255,7 +1274,7 @@ namespace Dune {
_op->applyscaleadd(-1,x,b); // overwrite b with defect
std::vector<std::shared_ptr<X> > p(_restart);
- std::vector<typename X::field_type> pp(_restart);
+ std::vector<field_type,fAlloc> pp(_restart);
X q(x); // a temporary vector
X prec_res(x); // a temporary vector for preconditioner output
diff --git a/dune/istl/test/CMakeLists.txt b/dune/istl/test/CMakeLists.txt
index a78ee0b83118583ad19faa45f63b7c671e5a51dd..a0df50b87f3360c4dec88f2694d076645f789976 100644
--- a/dune/istl/test/CMakeLists.txt
+++ b/dune/istl/test/CMakeLists.txt
@@ -20,10 +20,8 @@ dune_add_test(SOURCES matrixutilstest.cc)
dune_add_test(SOURCES matrixtest.cc)
-if(Vc_FOUND)
dune_add_test(SOURCES multirhstest.cc)
add_dune_vc_flags(multirhstest)
-endif(Vc_FOUND)
dune_add_test(SOURCES bvectortest.cc)
diff --git a/dune/istl/test/laplacian.hh b/dune/istl/test/laplacian.hh
index fc2ff218e495068a3f5e07583e1477632d4df276..a02ca932dba8b4f417b9ee5dfacf43ce9e612eef 100644
--- a/dune/istl/test/laplacian.hh
+++ b/dune/istl/test/laplacian.hh
@@ -5,14 +5,14 @@
#include <dune/istl/bcrsmatrix.hh>
#include <dune/common/fvector.hh>
-template<class B>
-void setupSparsityPattern(Dune::BCRSMatrix<B>& A, int N)
+template<class B, class Alloc>
+void setupSparsityPattern(Dune::BCRSMatrix<B,Alloc>& A, int N)
{
- typedef typename Dune::BCRSMatrix<B> Matrix;
+ typedef typename Dune::BCRSMatrix<B,Alloc> Matrix;
A.setSize(N*N, N*N, N*N*5);
A.setBuildMode(Matrix::row_wise);
- for (typename Dune::BCRSMatrix<B>::CreateIterator i = A.createbegin(); i != A.createend(); ++i) {
+ for (typename Dune::BCRSMatrix<B,Alloc>::CreateIterator i = A.createbegin(); i != A.createend(); ++i) {
int x = i.index()%N; // x coordinate in the 2d field
int y = i.index()/N; // y coordinate in the 2d field
@@ -36,10 +36,10 @@ void setupSparsityPattern(Dune::BCRSMatrix<B>& A, int N)
}
-template<class B>
-void setupLaplacian(Dune::BCRSMatrix<B>& A, int N)
+template<class B, class Alloc>
+void setupLaplacian(Dune::BCRSMatrix<B,Alloc>& A, int N)
{
- typedef typename Dune::BCRSMatrix<B>::field_type FieldType;
+ typedef typename Dune::BCRSMatrix<B,Alloc>::field_type FieldType;
setupSparsityPattern(A,N);
@@ -53,7 +53,7 @@ void setupLaplacian(Dune::BCRSMatrix<B>& A, int N)
b->operator[](b.index())=-1.0;
- for (typename Dune::BCRSMatrix<B>::RowIterator i = A.begin(); i != A.end(); ++i) {
+ for (typename Dune::BCRSMatrix<B,Alloc>::RowIterator i = A.begin(); i != A.end(); ++i) {
int x = i.index()%N; // x coordinate in the 2d field
int y = i.index()/N; // y coordinate in the 2d field
@@ -92,9 +92,9 @@ void setupLaplacian(Dune::BCRSMatrix<B>& A, int N)
}
}
}
-template<int BS>
-void setBoundary(Dune::BlockVector<Dune::FieldVector<double,BS> >& lhs,
- Dune::BlockVector<Dune::FieldVector<double,BS> >& rhs,
+template<int BS, class Alloc>
+void setBoundary(Dune::BlockVector<Dune::FieldVector<double,BS>, Alloc>& lhs,
+ Dune::BlockVector<Dune::FieldVector<double,BS>, Alloc>& rhs,
int N)
{
for(int i=0; i < lhs.size(); ++i) {
diff --git a/dune/istl/test/multirhstest.cc b/dune/istl/test/multirhstest.cc
index 58b28ddb713ba7fffe25f47ce19aab20f928bdb7..da3a6667046bcbaf8ec01f66df4c7be2c0bf1e79 100644
--- a/dune/istl/test/multirhstest.cc
+++ b/dune/istl/test/multirhstest.cc
@@ -16,7 +16,9 @@
#include <cmath> // Yes, we do some math here
#include <sys/times.h> // for timing measurements
+#include <dune/common/alignedallocator.hh>
#include <dune/common/classname.hh>
+#include <dune/common/debugalign.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/simd.hh>
@@ -117,8 +119,10 @@ void test_all(unsigned int Runs = 1)
typedef typename Dune::SimdScalarTypeTraits<FT>::type MT;
typedef Dune::FieldVector<FT,1> VB;
typedef Dune::FieldMatrix<MT,1,1> MB;
- typedef Dune::BlockVector<VB> Vector;
- typedef Dune::BCRSMatrix<MB> Matrix;
+ typedef Dune::AlignedAllocator<VB> AllocV;
+ typedef Dune::AlignedAllocator<MB> AllocM;
+ typedef Dune::BlockVector<VB,AllocV> Vector;
+ typedef Dune::BCRSMatrix<MB,AllocM> Matrix;
// size
unsigned int size = 100;
@@ -158,7 +162,7 @@ void test_all(unsigned int Runs = 1)
criterion.setNoPreSmoothSteps(1);
criterion.setNoPostSmoothSteps(1);
Dune::SeqScalarProduct<Vector> sp;
- typedef Dune::Amg::AMG<Operator,Vector,Smoother> AMG;
+ typedef Dune::Amg::AMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation> AMG;
Smoother smoother(A,1,1);
AMG amg(op, criterion, smootherArgs);
@@ -176,6 +180,9 @@ int main (int argc, char ** argv)
{
test_all<float>();
test_all<double>();
+
+ test_all<Dune::AlignedNumber<double> >();
+
#if HAVE_VC
test_all<Vc::float_v>();
test_all<Vc::double_v>();