diff --git a/dune/istl/multitypeblockmatrix.hh b/dune/istl/multitypeblockmatrix.hh
index 84d267a309396832f253c21d293e2432bbecc43e..7fc8f3d734fb47ee1ab68af47ffe39358fa94979 100644
--- a/dune/istl/multitypeblockmatrix.hh
+++ b/dune/istl/multitypeblockmatrix.hh
@@ -7,6 +7,8 @@
#include <iostream>
#include <tuple>
+#include <dune/common/hybridutilities.hh>
+
#include "istlexception.hh"
// forward declaration
@@ -31,163 +33,6 @@ namespace Dune {
- /**
- @brief prints out a matrix (type MultiTypeBlockMatrix)
-
- The parameters "crow" and "ccol" are the indices of
- the element to be printed out. Via recursive calling
- all other elements are printed, too. This is internally
- called when a MultiTypeBlockMatrix object is passed to an output stream.
- */
-
- template<int crow, int remain_rows, int ccol, int remain_cols,
- typename TMatrix>
- class MultiTypeBlockMatrix_Print {
- public:
-
- /**
- * print out a matrix block and all following
- */
- static void print(const TMatrix& m) {
- std::cout << "\t(" << crow << ", " << ccol << "): \n" << std::get<ccol>( std::get<crow>(m));
- MultiTypeBlockMatrix_Print<crow,remain_rows,ccol+1,remain_cols-1,TMatrix>::print(m); //next column
- }
- };
- template<int crow, int remain_rows, int ccol, typename TMatrix> //specialization for remain_cols=0
- class MultiTypeBlockMatrix_Print<crow,remain_rows,ccol,0,TMatrix> {
- public: static void print(const TMatrix& m) {
- MultiTypeBlockMatrix_Print<crow+1,remain_rows-1,0,TMatrix::N(),TMatrix>::print(m); //next row
- }
- };
-
- template<int crow, int ccol, int remain_cols, typename TMatrix> //recursion end: specialization for remain_rows=0
- class MultiTypeBlockMatrix_Print<crow,0,ccol,remain_cols,TMatrix> {
- public:
- static void print(const TMatrix& m)
- {std::cout << std::endl;}
- };
-
-
-
- //make MultiTypeBlockVector_Ident known (for MultiTypeBlockMatrix_Ident)
- template<int count, typename T1, typename T2>
- class MultiTypeBlockVector_Ident;
-
-
- /**
- @brief Set a MultiTypeBlockMatrix to some specific scalar value
-
- This class is used by the MultiTypeBlockMatrix class' internal assignment operator.
- Whenever a vector is assigned a scalar value, each block element
- has to provide operator= for the right side. Example:
- \code
- typedef MultiTypeBlockVector<int,int,int> CVect;
- MultiTypeBlockMatrix<CVect,CVect> CMat;
- CMat = 3; //sets all 3x2 integer elements to 3
- \endcode
- */
- template<int rowcount, typename T1, typename T2>
- class MultiTypeBlockMatrix_Ident {
- public:
-
- /**
- * equalize two matrix' element
- * note: uses MultiTypeBlockVector_Ident to equalize each row (which is of MultiTypeBlockVector type)
- */
- static void equalize(T1& a, const T2& b) {
- MultiTypeBlockVector_Ident<T1::M(),T1,T2>::equalize(a,b); //rows are cvectors
- MultiTypeBlockMatrix_Ident<rowcount-1,T1,T2>::equalize(a,b); //iterate over rows
- }
- };
-
- //recursion end for rowcount=0
- template<typename T1, typename T2>
- class MultiTypeBlockMatrix_Ident<0,T1,T2> {
- public:
- static void equalize (T1&, const T2&)
- {}
- };
-
- /**
- @brief Matrix-vector multiplication
-
- This class implements matrix vector multiplication for MultiTypeBlockMatrix/MultiTypeBlockVector types
- */
- template<int crow, int remain_rows, int ccol, int remain_cols,
- typename TVecY, typename TMatrix, typename TVecX>
- class MultiTypeBlockMatrix_VectMul {
- public:
-
- /** \brief y += A x
- */
- static void umv(TVecY& y, const TMatrix& A, const TVecX& x) {
- std::get<ccol>( std::get<crow>(A) ).umv( std::get<ccol>(x), std::get<crow>(y) );
- MultiTypeBlockMatrix_VectMul<crow,remain_rows,ccol+1,remain_cols-1,TVecY,TMatrix,TVecX>::umv(y, A, x);
- }
-
- /** \brief y -= A x
- */
- static void mmv(TVecY& y, const TMatrix& A, const TVecX& x) {
- std::get<ccol>( std::get<crow>(A) ).mmv( std::get<ccol>(x), std::get<crow>(y) );
- MultiTypeBlockMatrix_VectMul<crow,remain_rows,ccol+1,remain_cols-1,TVecY,TMatrix,TVecX>::mmv(y, A, x);
- }
-
- /** \brief y += alpha A x
- * \tparam AlphaType Type used for the scalar factor 'alpha'
- */
- template<typename AlphaType>
- static void usmv(const AlphaType& alpha, TVecY& y, const TMatrix& A, const TVecX& x) {
- std::get<ccol>( std::get<crow>(A) ).usmv(alpha, std::get<ccol>(x), std::get<crow>(y) );
- MultiTypeBlockMatrix_VectMul<crow,remain_rows,ccol+1,remain_cols-1,TVecY,TMatrix,TVecX>::usmv(alpha,y, A, x);
- }
-
-
- };
-
- //specialization for remain_cols = 0
- template<int crow, int remain_rows,int ccol, typename TVecY,
- typename TMatrix, typename TVecX>
- class MultiTypeBlockMatrix_VectMul<crow,remain_rows,ccol,0,TVecY,TMatrix,TVecX> { //start iteration over next row
-
- public:
- /**
- * do y += A x in next row
- */
- static void umv(TVecY& y, const TMatrix& A, const TVecX& x) {
- MultiTypeBlockMatrix_VectMul<crow+1,remain_rows-1,0,TMatrix::M(),TVecY,TMatrix,TVecX>::umv(y, A, x);
- }
-
- /**
- * do y -= A x in next row
- */
- static void mmv(TVecY& y, const TMatrix& A, const TVecX& x) {
- MultiTypeBlockMatrix_VectMul<crow+1,remain_rows-1,0,TMatrix::M(),TVecY,TMatrix,TVecX>::mmv(y, A, x);
- }
-
- template <typename AlphaType>
- static void usmv(const AlphaType& alpha, TVecY& y, const TMatrix& A, const TVecX& x) {
- MultiTypeBlockMatrix_VectMul<crow+1,remain_rows-1,0,TMatrix::M(),TVecY,TMatrix,TVecX>::usmv(alpha,y, A, x);
- }
- };
-
- //specialization for remain_rows = 0
- template<int crow, int ccol, int remain_cols, typename TVecY,
- typename TMatrix, typename TVecX>
- class MultiTypeBlockMatrix_VectMul<crow,0,ccol,remain_cols,TVecY,TMatrix,TVecX> {
- //end recursion
- public:
- static void umv(TVecY&, const TMatrix&, const TVecX&) {}
- static void mmv(TVecY&, const TMatrix&, const TVecX&) {}
-
- template<typename AlphaType>
- static void usmv(const AlphaType&, TVecY&, const TMatrix&, const TVecX&) {}
- };
-
-
-
-
-
-
/**
@brief A Matrix class to support different block types
@@ -259,7 +104,16 @@ namespace Dune {
* assignment operator
*/
template<typename T>
- void operator= (const T& newval) {MultiTypeBlockMatrix_Ident<N(),type,T>::equalize(*this, newval); }
+ void operator= (const T& newval) {
+ using namespace Dune::Hybrid;
+ auto size = index_constant<1+sizeof...(Args)>();
+ // Since Dune::Hybrid::size(MultiTypeBlockMatrix) is not implemented,
+ // we cannot use a plain forEach(*this, ...). This could be achieved,
+ // e.g., by implementing a static size() method.
+ forEach(integralRange(size), [&](auto&& i) {
+ (*this)[i] = newval;
+ });
+ }
/** \brief y = A x
*/
@@ -268,7 +122,7 @@ namespace Dune {
static_assert(X::size() == M(), "length of x does not match row length");
static_assert(Y::size() == N(), "length of y does not match row count");
y = 0; //reset y (for mv uses umv)
- MultiTypeBlockMatrix_VectMul<0,N(),0,M(),Y,type,X>::umv(y, *this, x); //iterate over all matrix elements
+ umv(x,y);
}
/** \brief y += A x
@@ -277,7 +131,12 @@ namespace Dune {
void umv (const X& x, Y& y) const {
static_assert(X::size() == M(), "length of x does not match row length");
static_assert(Y::size() == N(), "length of y does not match row count");
- MultiTypeBlockMatrix_VectMul<0,N(),0,M(),Y,type,X>::umv(y, *this, x); //iterate over all matrix elements
+ using namespace Dune::Hybrid;
+ forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
+ forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
+ (*this)[i][j].umv(x[j], y[i]);
+ });
+ });
}
/** \brief y -= A x
@@ -286,7 +145,12 @@ namespace Dune {
void mmv (const X& x, Y& y) const {
static_assert(X::size() == M(), "length of x does not match row length");
static_assert(Y::size() == N(), "length of y does not match row count");
- MultiTypeBlockMatrix_VectMul<0,N(),0,M(),Y,type,X>::mmv(y, *this, x); //iterate over all matrix elements
+ using namespace Dune::Hybrid;
+ forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
+ forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
+ (*this)[i][j].mmv(x[j], y[i]);
+ });
+ });
}
/** \brief y += alpha A x
@@ -295,7 +159,12 @@ namespace Dune {
void usmv (const AlphaType& alpha, const X& x, Y& y) const {
static_assert(X::size() == M(), "length of x does not match row length");
static_assert(Y::size() == N(), "length of y does not match row count");
- MultiTypeBlockMatrix_VectMul<0,N(),0,M(),Y,type,X>::usmv(alpha,y, *this, x); //iterate over all matrix elements
+ using namespace Dune::Hybrid;
+ forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
+ forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
+ (*this)[i][j].usmv(alpha, x[j], y[i]);
+ });
+ });
}
};
@@ -309,9 +178,15 @@ namespace Dune {
*/
template<typename T1, typename... Args>
std::ostream& operator<< (std::ostream& s, const MultiTypeBlockMatrix<T1,Args...>& m) {
- static const int N = MultiTypeBlockMatrix<T1,Args...>::N();
- static const int M = MultiTypeBlockMatrix<T1,Args...>::M();
- MultiTypeBlockMatrix_Print<0,N,0,M,MultiTypeBlockMatrix<T1,Args...> >::print(m);
+ auto N = index_constant<MultiTypeBlockMatrix<T1,Args...>::N()>();
+ auto M = index_constant<MultiTypeBlockMatrix<T1,Args...>::M()>();
+ using namespace Dune::Hybrid;
+ forEach(integralRange(N), [&](auto&& i) {
+ forEach(integralRange(M), [&](auto&& j) {
+ s << "\t(" << i << ", " << j << "): \n" << m[i][j];
+ });
+ });
+ s << std::endl;
return s;
}
diff --git a/dune/istl/multitypeblockvector.hh b/dune/istl/multitypeblockvector.hh
index d718cb3f5cde44f8bcd668c290e541976f8e54ab..0865072ac98b37d9dedf525920637dfc6ee1341f 100644
--- a/dune/istl/multitypeblockvector.hh
+++ b/dune/istl/multitypeblockvector.hh
@@ -9,6 +9,7 @@
#include <dune/common/dotproduct.hh>
#include <dune/common/ftraits.hh>
+#include <dune/common/hybridutilities.hh>
#include "istlexception.hh"
@@ -30,266 +31,6 @@ namespace Dune {
- /**
- @brief prints out a vector (type MultiTypeBlockVector)
-
- The parameter "current_element" is the index of
- the element to be printed out. Via recursive calling
- all other elements (number is "remaining_elements")
- are printed, too. This is internally called when
- a MultiTypeBlockVector object is passed to an output stream.
- Example:
- \code
- MultiTypeBlockVector<int,int> CVect;
- std::cout << CVect;
- \endcode
- */
-
- template<int current_element, int remaining_elements, typename TVec>
- class MultiTypeBlockVector_Print {
- public:
- /**
- * print out the current vector element and all following
- */
- static void print(const TVec& v) {
- std::cout << "\t(" << current_element << "):\n" << std::get<current_element>(v) << "\n";
- MultiTypeBlockVector_Print<current_element+1,remaining_elements-1,TVec>::print(v); //next element
- }
- };
- template<int current_element, typename TVec> //recursion end (remaining_elements=0)
- class MultiTypeBlockVector_Print<current_element,0,TVec> {
- public:
- static void print(const TVec&) {std::cout << "\n";}
- };
-
-
-
- /**
- @brief Set a MultiTypeBlockVector to some specific value
-
- This class is used by the MultiTypeBlockVector class' internal = operator.
- Whenever a vector is assigned to a value, each vector element
- has to provide the = operator for the right side. Example:
- \code
- MultiTypeBlockVector<int,int,int> CVect;
- CVect = 3; //sets all integer elements to 3
- \endcode
- */
- template<int count, typename T1, typename T2>
- class MultiTypeBlockVector_Ident {
- public:
-
- /**
- * equalize two vectors' element (index is template parameter count)
- * note: each MultiTypeBlockVector element has to provide the = operator with type T2
- */
- static void equalize(T1& a, const T2& b) {
- std::get<count-1>(a) = b; //equalize current elements
- MultiTypeBlockVector_Ident<count-1,T1,T2>::equalize(a,b); //next elements
- }
- };
- template<typename T1, typename T2> //recursion end (count=0)
- class MultiTypeBlockVector_Ident<0,T1,T2> {public: static void equalize (T1&, const T2&) {} };
-
-
-
-
-
-
- /**
- @brief Add/subtract second vector to/from the first (v1 += v2)
-
- This class implements vector addition/subtraction for any MultiTypeBlockVector-Class type.
- */
- template<int count, typename T>
- class MultiTypeBlockVector_Add {
- public:
-
- /**
- * add vector to vector
- */
- static void add (T& a, const T& b) { //add vector elements
- std::get<(count-1)>(a) += std::get<(count-1)>(b);
- MultiTypeBlockVector_Add<count-1,T>::add(a,b);
- }
-
- /**
- * Subtract vector from vector
- */
- static void sub (T& a, const T& b) { //sub vector elements
- std::get<(count-1)>(a) -= std::get<(count-1)>(b);
- MultiTypeBlockVector_Add<count-1,T>::sub(a,b);
- }
- };
- template<typename T> //recursion end; specialization for count=0
- class MultiTypeBlockVector_Add<0,T> {public: static void add (T&, const T&) {} static void sub (T&, const T&) {} };
-
-
-
- /**
- @brief AXPY operation on vectors
-
- calculates x += a * y
- */
- template<int count, typename TVec, typename Ta>
- class MultiTypeBlockVector_AXPY {
- public:
-
- /**
- * calculates x += a * y
- */
- static void axpy(TVec& x, const Ta& a, const TVec& y) {
- std::get<(count-1)>(x).axpy(a,std::get<(count-1)>(y));
- MultiTypeBlockVector_AXPY<count-1,TVec,Ta>::axpy(x,a,y);
- }
- };
- template<typename TVec, typename Ta> //specialization for count=0
- class MultiTypeBlockVector_AXPY<0,TVec,Ta> {public: static void axpy (TVec&, const Ta&, const TVec&) {} };
-
-
- /** @brief In-place multiplication with a scalar
- *
- * Calculates v *= a for each element of the given vector.
- */
- template<int count, typename TVec, typename Ta>
- class MultiTypeBlockVector_Mulscal {
- public:
-
- /**
- * calculates x *= a
- */
- static void mul(TVec& x, const Ta& a) {
- std::get<(count-1)>(x) *= a;
- MultiTypeBlockVector_Mulscal<count-1,TVec,Ta>::mul(x,a);
- }
- };
- template<typename TVec, typename Ta> //specialization for count=0
- class MultiTypeBlockVector_Mulscal<0,TVec,Ta> {public: static void mul (TVec&, const Ta&) {} };
-
-
-
- /** @brief Scalar products
- *
- * multiplies the current elements of x and y pairwise, and sum up the results.
- * Provides two variants:
- * 1) 'mul' computes the indefinite inner product and
- * 2) 'dot' provides an inner product by conjugating the first argument
- */
- template<int count, typename TVec>
- class MultiTypeBlockVector_Mul {
- public:
- static typename TVec::field_type mul(const TVec& x, const TVec& y)
- {
- return (std::get<count-1>(x) * std::get<count-1>(y)) + MultiTypeBlockVector_Mul<count-1,TVec>::mul(x,y);
- }
-
- static typename TVec::field_type dot(const TVec& x, const TVec& y)
- {
- return Dune::dot(std::get<count-1>(x),std::get<count-1>(y)) + MultiTypeBlockVector_Mul<count-1,TVec>::dot(x,y);
- }
- };
-
- template<typename TVec>
- class MultiTypeBlockVector_Mul<0,TVec> {
- public:
- static typename TVec::field_type mul(const TVec&, const TVec&) {return 0;}
- static typename TVec::field_type dot(const TVec&, const TVec&) {return 0;}
- };
-
-
-
-
-
- /** \brief Calculate the 2-norm
-
- Each element of the vector has to provide the method "two_norm2()"
- in order to calculate the whole vector's 2-norm.
- */
- template<int count, typename T>
- class MultiTypeBlockVector_Norm {
- public:
- typedef typename T::field_type field_type;
- typedef typename FieldTraits<field_type>::real_type real_type;
-
- /**
- * sum up all elements' 2-norms
- */
- static real_type result (const T& a) { //result = sum of all elements' 2-norms
- return std::get<count-1>(a).two_norm2() + MultiTypeBlockVector_Norm<count-1,T>::result(a);
- }
- };
-
- template<typename T> //recursion end: no more vector elements to add...
- class MultiTypeBlockVector_Norm<0,T> {
- public:
- typedef typename T::field_type field_type;
- typedef typename FieldTraits<field_type>::real_type real_type;
- static real_type result (const T&) {return 0.0;}
- };
-
- /** \brief Calculate the \infty-norm
-
- Each element of the vector has to provide the method "infinity_norm()"
- in order to calculate the whole vector's norm.
- */
- template<int count, typename T,
- bool hasNaN = has_nan<typename T::field_type>::value>
- class MultiTypeBlockVector_InfinityNorm {
- public:
- typedef typename T::field_type field_type;
- typedef typename FieldTraits<field_type>::real_type real_type;
-
- /**
- * Take the maximum over all elements' norms
- */
- static real_type result (const T& a)
- {
- std::pair<real_type, real_type> const ret = resultHelper(a);
- return ret.first * (ret.second / ret.second);
- }
-
- // Computes the maximum while keeping track of NaN values
- static std::pair<real_type, real_type> resultHelper(const T& a) {
- using std::max;
- auto const rest =
- MultiTypeBlockVector_InfinityNorm<count - 1, T>::resultHelper(a);
- real_type const norm = std::get<count - 1>(a).infinity_norm();
- return {max(norm, rest.first), norm + rest.second};
- }
- };
-
- template <int count, typename T>
- class MultiTypeBlockVector_InfinityNorm<count, T, false> {
- public:
- typedef typename T::field_type field_type;
- typedef typename FieldTraits<field_type>::real_type real_type;
-
- /**
- * Take the maximum over all elements' norms
- */
- static real_type result (const T& a)
- {
- using std::max;
- return max(std::get<count-1>(a).infinity_norm(), MultiTypeBlockVector_InfinityNorm<count-1,T>::result(a));
- }
- };
-
- template<typename T, bool hasNaN> //recursion end
- class MultiTypeBlockVector_InfinityNorm<0,T, hasNaN> {
- public:
- typedef typename T::field_type field_type;
- typedef typename FieldTraits<field_type>::real_type real_type;
- static real_type result (const T&)
- {
- return 0.0;
- }
-
- static std::pair<real_type,real_type> resultHelper (const T&)
- {
- return {0.0, 1.0};
- }
- };
-
/** @addtogroup DenseMatVec
@{
*/
@@ -386,28 +127,83 @@ namespace Dune {
/** \brief Assignment operator
*/
template<typename T>
- void operator= (const T& newval) {MultiTypeBlockVector_Ident<sizeof...(Args),type,T>::equalize(*this, newval); }
+ void operator= (const T& newval) {
+ Dune::Hybrid::forEach(*this, [&](auto&& entry) {
+ entry = newval;
+ });
+ }
/**
* operator for MultiTypeBlockVector += MultiTypeBlockVector operations
*/
- void operator+= (const type& newv) {MultiTypeBlockVector_Add<sizeof...(Args),type>::add(*this,newv);}
+ void operator+= (const type& newv) {
+ using namespace Dune::Hybrid;
+ forEach(integralRange(Hybrid::size(*this)), [&](auto&& i) {
+ (*this)[i] += newv[i];
+ });
+ }
/**
* operator for MultiTypeBlockVector -= MultiTypeBlockVector operations
*/
- void operator-= (const type& newv) {MultiTypeBlockVector_Add<sizeof...(Args),type>::sub(*this,newv);}
+ void operator-= (const type& newv) {
+ using namespace Dune::Hybrid;
+ forEach(integralRange(Hybrid::size(*this)), [&](auto&& i) {
+ (*this)[i] -= newv[i];
+ });
+ }
+
+ // Once we have the IsNumber traits class the following
+ // three implementations could be replaced by:
+ //
+ // template<class T,
+ // std::enable_if_t< IsNumber<T>::value, int> = 0>
+ // void operator*= (const T& w) {
+ // Dune::Hybrid::forEach(*this, [&](auto&& entry) {
+ // entry *= w;
+ // });
+ // }
+
+ void operator*= (const int& w) {
+ Dune::Hybrid::forEach(*this, [&](auto&& entry) {
+ entry *= w;
+ });
+ }
- void operator*= (const int& w) {MultiTypeBlockVector_Mulscal<sizeof...(Args),type,const int>::mul(*this,w);}
- void operator*= (const float& w) {MultiTypeBlockVector_Mulscal<sizeof...(Args),type,const float>::mul(*this,w);}
- void operator*= (const double& w) {MultiTypeBlockVector_Mulscal<sizeof...(Args),type,const double>::mul(*this,w);}
+ void operator*= (const float& w) {
+ Dune::Hybrid::forEach(*this, [&](auto&& entry) {
+ entry *= w;
+ });
+ }
+
+ void operator*= (const double& w) {
+ Dune::Hybrid::forEach(*this, [&](auto&& entry) {
+ entry *= w;
+ });
+ }
+
+ field_type operator* (const type& newv) const {
+ using namespace Dune::Hybrid;
+ return accumulate(integralRange(Hybrid::size(*this)), field_type(0), [&](auto&& a, auto&& i) {
+ return a + (*this)[i]*newv[i];
+ });
+ }
- field_type operator* (const type& newv) const {return MultiTypeBlockVector_Mul<sizeof...(Args),type>::mul(*this,newv);}
- field_type dot (const type& newv) const {return MultiTypeBlockVector_Mul<sizeof...(Args),type>::dot(*this,newv);}
+ field_type dot (const type& newv) const {
+ using namespace Dune::Hybrid;
+ return accumulate(integralRange(Hybrid::size(*this)), field_type(0), [&](auto&& a, auto&& i) {
+ return a + (*this)[i].dot(newv[i]);
+ });
+ }
/** \brief Compute the squared Euclidean norm
*/
- typename FieldTraits<field_type>::real_type two_norm2() const {return MultiTypeBlockVector_Norm<sizeof...(Args),type>::result(*this);}
+ typename FieldTraits<field_type>::real_type two_norm2() const {
+ using namespace Dune::Hybrid;
+ return accumulate(*this, typename FieldTraits<field_type>::real_type(0), [&](auto&& a, auto&& entry) {
+ return a + entry.two_norm2();
+ });
+ }
/** \brief Compute the Euclidean norm
*/
@@ -417,7 +213,29 @@ namespace Dune {
*/
typename FieldTraits<field_type>::real_type infinity_norm() const
{
- return MultiTypeBlockVector_InfinityNorm<sizeof...(Args),type>::result(*this);
+ using namespace Dune::Hybrid;
+ using std::max;
+ using real_type = typename FieldTraits<field_type>::real_type;
+
+ real_type result = 0.0;
+ // Compute max norm tracking appearing nan values
+ // if the field type supports nan.
+ ifElse(has_nan<field_type>(), [&](auto&& id) {
+ // This variable will preserve any nan value
+ real_type nanTracker = 1.0;
+ forEach(*this, [&](auto&& entry) {
+ real_type entryNorm = entry.infinity_norm();
+ result = max(entryNorm, result);
+ nanTracker += entryNorm;
+ });
+ // Incorporate possible nan value into result
+ result *= (nanTracker / nanTracker);
+ }, [&](auto&& id) {
+ forEach(*this, [&](auto&& entry) {
+ result = max(entry.infinity_norm(), result);
+ });
+ });
+ return result;
}
/** \brief Axpy operation on this vector (*this += a * y)
@@ -426,7 +244,10 @@ namespace Dune {
*/
template<typename Ta>
void axpy (const Ta& a, const type& y) {
- MultiTypeBlockVector_AXPY<sizeof...(Args),type,Ta>::axpy(*this,a,y);
+ using namespace Dune::Hybrid;
+ forEach(integralRange(Hybrid::size(*this)), [&](auto&& i) {
+ (*this)[i].axpy(a, y[i]);
+ });
}
};
@@ -437,7 +258,10 @@ namespace Dune {
*/
template <typename... Args>
std::ostream& operator<< (std::ostream& s, const MultiTypeBlockVector<Args...>& v) {
- MultiTypeBlockVector_Print<0,sizeof...(Args),MultiTypeBlockVector<Args...> >::print(v);
+ using namespace Dune::Hybrid;
+ forEach(integralRange(size(v)), [&](auto&& i) {
+ s << "\t(" << i << "):\n" << v[i] << "\n";
+ });
return s;
}
diff --git a/dune/istl/test/multitypeblockmatrixtest.cc b/dune/istl/test/multitypeblockmatrixtest.cc
index b1d296182365d1f369b9ac03b49e0ea1c714302c..359393d46507b2d7eb73fba8a8a5534ff3e24cd8 100644
--- a/dune/istl/test/multitypeblockmatrixtest.cc
+++ b/dune/istl/test/multitypeblockmatrixtest.cc
@@ -71,10 +71,22 @@ int main(int argc, char** argv) try
result[_0].resize(3);
result[_1].resize(2);
+ multiMatrix[_0][_0] = 4200;
+ multiMatrix[_0][_1] = 4201;
+ multiMatrix[_1][_0] = 4210;
+ multiMatrix[_1][_1] = 4211;
+
multiMatrix.mv(multiVector,result);
+ std::cout << "mv result" << std::endl << result << std::endl;
+
multiMatrix.umv(multiVector,result);
+ std::cout << "umv result" << std::endl << result << std::endl;
+
multiMatrix.mmv(multiVector,result);
+ std::cout << "mmv result" << std::endl << result << std::endl;
+
multiMatrix.usmv(3.14,multiVector,result);
+ std::cout << "usmv result" << std::endl << result << std::endl;
///////////////////////////////////////////////////////////////////////////////////////////////////////////
// Next test: Set up a linear system with a matrix consisting of 2x2 sparse scalar matrices.