diff --git a/patches/dune-common/0001-Add-lots-of-exadune-stuff-without-thinking-about-it.patch b/patches/dune-common/0001-Add-lots-of-exadune-stuff-without-thinking-about-it.patch deleted file mode 100644 index 6566c2b13c15f23847a194b0da0dd2f1c29850ea..0000000000000000000000000000000000000000 --- a/patches/dune-common/0001-Add-lots-of-exadune-stuff-without-thinking-about-it.patch +++ /dev/null @@ -1,597 +0,0 @@ -From 5db32867b4889e25b43747dfbd888164ade291c3 Mon Sep 17 00:00:00 2001 -From: =?UTF-8?q?Ren=C3=A9=20He=C3=9F?= <rene.hess@iwr.uni-heidelberg.de> -Date: Thu, 22 Dec 2016 16:55:13 +0100 -Subject: [PATCH] Add lots of exadune stuff without thinking about it - -This commit is potentially very damaging since it was done without -checking if the changes make sense. ---- - dune/common/densematrix.hh | 251 ++++++++++++--------------------------------- - dune/common/dynmatrix.hh | 25 +---- - dune/common/fmatrix.hh | 51 +++++---- - 3 files changed, 94 insertions(+), 233 deletions(-) - -diff --git a/dune/common/densematrix.hh b/dune/common/densematrix.hh -index db774fc9..7daf4a3a 100644 ---- a/dune/common/densematrix.hh -+++ b/dune/common/densematrix.hh -@@ -16,8 +16,8 @@ - #include <dune/common/fvector.hh> - #include <dune/common/math.hh> - #include <dune/common/precision.hh> --#include <dune/common/simd.hh> --#include <dune/common/typetraits.hh> -+#include <dune/common/classname.hh> -+#include <dune/common/math.hh> - #include <dune/common/unused.hh> - - namespace Dune -@@ -75,15 +75,14 @@ namespace Dune - struct DenseMatrixAssigner; - - #ifndef DOXYGEN -- namespace Impl -+ namespace - { -- -- template< class DenseMatrix, class RHS, class = void > -- class DenseMatrixAssigner -- {}; -+ template< class DenseMatrix, class RHS, -+ bool primitive = std::is_convertible< RHS, typename DenseMatrix::field_type >::value > -+ class DenseMatrixAssignerImplementation; - - template< class DenseMatrix, class RHS > -- class DenseMatrixAssigner< DenseMatrix, RHS, std::enable_if_t< Dune::IsNumber< RHS >::value > > -+ class DenseMatrixAssignerImplementation< DenseMatrix, RHS, true > - { - public: - static void apply ( DenseMatrix &denseMatrix, const RHS &rhs ) -@@ -94,45 +93,28 @@ namespace Dune - }; - - template< class DenseMatrix, class RHS > -- class DenseMatrixAssigner< DenseMatrix, RHS, std::enable_if_t< !std::is_same< typename RHS::const_iterator, void >::value > > -+ class DenseMatrixAssignerImplementation< DenseMatrix, RHS, false > - { - public: - static void apply ( DenseMatrix &denseMatrix, const RHS &rhs ) - { -- DUNE_ASSERT_BOUNDS(rhs.N() == denseMatrix.N()); -- DUNE_ASSERT_BOUNDS(rhs.M() == denseMatrix.M()); -- typename DenseMatrix::iterator tIt = std::begin(denseMatrix); -- typename RHS::const_iterator sIt = std::begin(rhs); -- for(; sIt != std::end(rhs); ++tIt, ++sIt) -- std::copy(std::begin(*sIt), std::end(*sIt), std::begin(*tIt)); -+ static_assert( (std::is_convertible< const RHS, const DenseMatrix >::value), -+ "No template specialization of DenseMatrixAssigner found" ); -+ denseMatrix = static_cast< const DenseMatrix & >( rhs ); - } - }; -- -- } // namespace Impl -+ } - - - - template< class DenseMatrix, class RHS > - struct DenseMatrixAssigner -- : public Impl::DenseMatrixAssigner< DenseMatrix, RHS > -- {}; -- -- -- namespace Impl - { -- -- template< class DenseMatrix, class RHS > -- std::true_type hasDenseMatrixAssigner ( DenseMatrix &, const RHS &, decltype( Dune::DenseMatrixAssigner< DenseMatrix, RHS >::apply( std::declval< DenseMatrix & >(), std::declval< const RHS & >() ) ) * = nullptr ); -- -- std::false_type hasDenseMatrixAssigner ( ... ); -- -- } // namespace Impl -- -- template< class DenseMatrix, class RHS > -- struct HasDenseMatrixAssigner -- : public decltype( Impl::hasDenseMatrixAssigner( std::declval< DenseMatrix & >(), std::declval< const RHS & >() ) ) -- {}; -- -+ static void apply ( DenseMatrix &denseMatrix, const RHS &rhs ) -+ { -+ DenseMatrixAssignerImplementation< DenseMatrix, RHS >::apply( denseMatrix, rhs ); -+ } -+ }; - #endif // #ifndef DOXYGEN - - -@@ -192,15 +174,6 @@ namespace Dune - blocklevel = 1 - }; - -- private: -- //! \brief if value_type is a simd vector, then this is a simd vector of -- //! the same length that can be used for indices. -- /** -- * Just pray that the fundamental type is actually large enough... -- */ -- using simd_index_type = SimdIndex<value_type>; -- -- public: - //===== access to components - - //! random access -@@ -293,7 +266,7 @@ namespace Dune - - //===== assignment - -- template< class RHS, class = std::enable_if_t< HasDenseMatrixAssigner< MAT, RHS >::value > > -+ template< class RHS > - DenseMatrix &operator= ( const RHS &rhs ) - { - DenseMatrixAssigner< MAT, RHS >::apply( asImp(), rhs ); -@@ -718,15 +691,15 @@ namespace Dune - #ifndef DOXYGEN - struct ElimPivot - { -- ElimPivot(std::vector<simd_index_type> & pivot); -+ ElimPivot(std::vector<size_type> & pivot); - -- void swap(std::size_t i, simd_index_type j); -+ void swap(int i, int j); - - template<typename T> - void operator()(const T&, int, int) - {} - -- std::vector<simd_index_type> & pivot_; -+ std::vector<size_type> & pivot_; - }; - - template<typename V> -@@ -734,7 +707,7 @@ namespace Dune - { - Elim(V& rhs); - -- void swap(std::size_t i, simd_index_type j); -+ void swap(int i, int j); - - void operator()(const typename V::field_type& factor, int k, int i); - -@@ -746,10 +719,8 @@ namespace Dune - ElimDet(field_type& sign) : sign_(sign) - { sign_ = 1; } - -- void swap(std::size_t i, simd_index_type j) -- { -- sign_ *= cond(simd_index_type(i) == j, field_type(1), field_type(-1)); -- } -+ void swap(int, int) -+ { sign_ *= -1; } - - void operator()(const field_type&, int, int) - {} -@@ -758,52 +729,13 @@ namespace Dune - }; - #endif // DOXYGEN - -- //! do an LU-Decomposition on matrix A -- /** -- * \param A The matrix to decompose, and to store the -- * result in. -- * \param func Functor used for swapping lanes and to conduct -- * the elimination. Depending on the functor, \c -- * luDecomposition() can be used for solving, for -- * inverting, or to compute the determinant. -- * \param nonsingularLanes SimdMask of lanes that are nonsingular. -- * \param throwEarly Whether to throw an \c FMatrixError immediately -- * as soon as one lane is discovered to be -- * singular. If \c false, do not throw, instead -- * continue until finished or all lanes are -- * singular, and exit via return in both cases. -- * -- * There are two modes of operation: -- * <ul> -- * <li>Terminate as soon as one lane is discovered to be singular. Early -- * termination is done by throwing an \c FMatrixError. On entry, \c -- * all_true(nonsingularLanes) and \c throwEarly==true should hold. -- * After early termination, the contents of \c A should be considered -- * bogus, and \c nonsingularLanes has the lane(s) that triggered the -- * early termination unset. There may be more singular lanes than the -- * one reported in \c nonsingularLanes, which just havent been -- * discovered yet; so the value of \c nonsingularLanes is mostly -- * useful for diagnostics.</li> -- * <li>Terminate only when all lanes are discovered to be singular. Use -- * this when you want to apply special postprocessing in singular -- * lines (e.g. setting the determinant of singular lanes to 0 in \c -- * determinant()). On entry, \c nonsingularLanes may have any value -- * and \c throwEarly==false should hold. The function will not throw -- * an exception if some lanes are discovered to be singular, instead -- * it will continue running until all lanes are singular or until -- * finished, and terminate only via normal return. On exit, \c -- * nonsingularLanes contains the map of lanes that are valid in \c -- * A.</li> -- * </ul> -- */ -- template<class Func, class Mask> -- void luDecomposition(DenseMatrix<MAT>& A, Func func, -- Mask &nonsingularLanes, bool throwEarly) const; -+ template<class Func> -+ void luDecomposition(DenseMatrix<MAT>& A, Func func) const; - }; - - #ifndef DOXYGEN - template<typename MAT> -- DenseMatrix<MAT>::ElimPivot::ElimPivot(std::vector<simd_index_type> & pivot) -+ DenseMatrix<MAT>::ElimPivot::ElimPivot(std::vector<size_type> & pivot) - : pivot_(pivot) - { - typedef typename std::vector<size_type>::size_type size_type; -@@ -811,9 +743,9 @@ namespace Dune - } - - template<typename MAT> -- void DenseMatrix<MAT>::ElimPivot::swap(std::size_t i, simd_index_type j) -+ void DenseMatrix<MAT>::ElimPivot::swap(int i, int j) - { -- assign(pivot_[i], j, SimdScalar<simd_index_type>(i) != j); -+ pivot_[i]=j; - } - - template<typename MAT> -@@ -824,14 +756,9 @@ namespace Dune - - template<typename MAT> - template<typename V> -- void DenseMatrix<MAT>::Elim<V>::swap(std::size_t i, simd_index_type j) -+ void DenseMatrix<MAT>::Elim<V>::swap(int i, int j) - { -- using std::swap; -- -- // see the comment in luDecomposition() -- for(std::size_t l = 0; l < lanes(j); ++l) -- swap(lane(l, (*rhs_)[ i ]), -- lane(l, (*rhs_)[lane(l, j)])); -+ std::swap((*rhs_)[i], (*rhs_)[j]); - } - - template<typename MAT> -@@ -841,78 +768,47 @@ namespace Dune - { - (*rhs_)[k] -= factor*(*rhs_)[i]; - } -- - template<typename MAT> -- template<typename Func, class Mask> -- inline void DenseMatrix<MAT>:: -- luDecomposition(DenseMatrix<MAT>& A, Func func, Mask &nonsingularLanes, -- bool throwEarly) const -+ template<typename Func> -+ inline void DenseMatrix<MAT>::luDecomposition(DenseMatrix<MAT>& A, Func func) const - { -- using std::swap; -- -- typedef typename FieldTraits<value_type>::real_type real_type; -- -+ typedef typename FieldTraits<value_type>::real_type -+ real_type; - real_type norm = A.infinity_norm_real(); // for relative thresholds -- real_type pivthres = -- std::max( FMatrixPrecision< real_type >::absolute_limit(), -- norm * FMatrixPrecision< real_type >::pivoting_limit() ); -- real_type singthres = -- std::max( FMatrixPrecision< real_type >::absolute_limit(), -- norm * FMatrixPrecision< real_type >::singular_limit() ); -+ real_type pivthres = std::max( FMatrixPrecision< real_type >::absolute_limit(), norm * FMatrixPrecision< real_type >::pivoting_limit() ); -+ real_type singthres = std::max( FMatrixPrecision< real_type >::absolute_limit(), norm * FMatrixPrecision< real_type >::singular_limit() ); - - // LU decomposition of A in A - for (size_type i=0; i<rows(); i++) // loop over all rows - { -- real_type pivmax = fvmeta::absreal(A[i][i]); -- auto do_pivot = pivmax<pivthres; -+ typename FieldTraits<value_type>::real_type pivmax=fvmeta::absreal(A[i][i]); - - // pivoting ? -- if (any_true(do_pivot && nonsingularLanes)) -+ if (pivmax<pivthres) - { - // compute maximum of column -- simd_index_type imax=i; -+ size_type imax=i; -+ typename FieldTraits<value_type>::real_type abs(0.0); - for (size_type k=i+1; k<rows(); k++) -- { -- auto abs = fvmeta::absreal(A[k][i]); -- auto mask = abs > pivmax && do_pivot; -- pivmax = cond(mask, abs, pivmax); -- imax = cond(mask, simd_index_type(k), imax); -- } -- // swap rows -- if (any_true(imax != i && nonsingularLanes)) { -- for (size_type j=0; j<rows(); j++) -+ if ((abs=fvmeta::absreal(A[k][i]))>pivmax) - { -- // This is a swap operation where the second operand is scattered, -- // and on top of that is also extracted from deep within a -- // moderately complicated data structure (a DenseMatrix), where we -- // can't assume much on the memory layout. On intel processors, -- // the only instruction that might help us here is vgather, but it -- // is unclear whether that is even faster than a software -- // implementation, and we would also need vscatter which does not -- // exist. So break vectorization here and do it manually. -- for(std::size_t l = 0; l < lanes(A[i][j]); ++l) -- swap(lane(l, A[i][j]), lane(l, A[lane(l, imax)][j])); -+ pivmax = abs; imax = k; - } -+ // swap rows -+ if (imax!=i) { -+ for (size_type j=0; j<rows(); j++) -+ std::swap(A[i][j],A[imax][j]); - func.swap(i, imax); // swap the pivot or rhs - } - } - - // singular ? -- nonsingularLanes = nonsingularLanes && !(pivmax<singthres); -- if (throwEarly) { -- if(!all_true(nonsingularLanes)) -- DUNE_THROW(FMatrixError, "matrix is singular"); -- } -- else { // !throwEarly -- if(!any_true(nonsingularLanes)) -- return; -- } -+ if (pivmax<singthres) -+ DUNE_THROW(FMatrixError,"matrix is singular"); - - // eliminate - for (size_type k=i+1; k<rows(); k++) - { -- // in the simd case, A[i][i] may be close to zero in some lanes. Pray -- // that the result is no worse than a quiet NaN. - field_type factor = A[k][i]/A[i][i]; - A[k][i] = factor; - for (size_type j=i+1; j<rows(); j++) -@@ -933,8 +829,7 @@ namespace Dune - if (rows()==1) { - - #ifdef DUNE_FMatrix_WITH_CHECKING -- if (any_true(fvmeta::absreal((*this)[0][0]) -- < FMatrixPrecision<>::absolute_limit())) -+ if (fvmeta::absreal((*this)[0][0])<FMatrixPrecision<>::absolute_limit()) - DUNE_THROW(FMatrixError,"matrix is singular"); - #endif - x[0] = b[0]/(*this)[0][0]; -@@ -944,8 +839,7 @@ namespace Dune - - field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0]; - #ifdef DUNE_FMatrix_WITH_CHECKING -- if (any_true(fvmeta::absreal(detinv) -- < FMatrixPrecision<>::absolute_limit())) -+ if (fvmeta::absreal(detinv)<FMatrixPrecision<>::absolute_limit()) - DUNE_THROW(FMatrixError,"matrix is singular"); - #endif - detinv = 1.0/detinv; -@@ -958,7 +852,7 @@ namespace Dune - - field_type d = determinant(); - #ifdef DUNE_FMatrix_WITH_CHECKING -- if (any_true(fvmeta::absreal(d) < FMatrixPrecision<>::absolute_limit())) -+ if (fvmeta::absreal(d)<FMatrixPrecision<>::absolute_limit()) - DUNE_THROW(FMatrixError,"matrix is singular"); - #endif - -@@ -981,10 +875,8 @@ namespace Dune - rhs = b; // copy data - Elim<V> elim(rhs); - MAT A(asImp()); -- SimdMask<typename FieldTraits<value_type>::real_type> -- nonsingularLanes(true); - -- luDecomposition(A, elim, nonsingularLanes, true); -+ luDecomposition(A, elim); - - // backsolve - for(int i=rows()-1; i>=0; i--) { -@@ -1005,8 +897,7 @@ namespace Dune - if (rows()==1) { - - #ifdef DUNE_FMatrix_WITH_CHECKING -- if (any_true(fvmeta::absreal((*this)[0][0]) -- < FMatrixPrecision<>::absolute_limit())) -+ if (fvmeta::absreal((*this)[0][0])<FMatrixPrecision<>::absolute_limit()) - DUNE_THROW(FMatrixError,"matrix is singular"); - #endif - (*this)[0][0] = field_type( 1 ) / (*this)[0][0]; -@@ -1016,8 +907,7 @@ namespace Dune - - field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0]; - #ifdef DUNE_FMatrix_WITH_CHECKING -- if (any_true(fvmeta::absreal(detinv) -- < FMatrixPrecision<>::absolute_limit())) -+ if (fvmeta::absreal(detinv)<FMatrixPrecision<>::absolute_limit()) - DUNE_THROW(FMatrixError,"matrix is singular"); - #endif - detinv = field_type( 1 ) / detinv; -@@ -1032,10 +922,8 @@ namespace Dune - else { - - MAT A(asImp()); -- std::vector<simd_index_type> pivot(rows()); -- SimdMask<typename FieldTraits<value_type>::real_type> -- nonsingularLanes(true); -- luDecomposition(A, ElimPivot(pivot), nonsingularLanes, true); -+ std::vector<size_type> pivot(rows()); -+ luDecomposition(A, ElimPivot(pivot)); - DenseMatrix<MAT>& L=A; - DenseMatrix<MAT>& U=A; - -@@ -1063,14 +951,9 @@ namespace Dune - - for(size_type i=rows(); i>0; ) { - --i; -- for(std::size_t l = 0; l < lanes((*this)[0][0]); ++l) -- { -- std::size_t pi = lane(l, pivot[i]); -- if(i!=pi) -- for(size_type j=0; j<rows(); ++j) -- std::swap(lane(l, (*this)[j][pi]), -- lane(l, (*this)[j][ i])); -- } -+ if(i!=pivot[i]) -+ for(size_type j=0; j<rows(); ++j) -+ std::swap((*this)[j][pivot[i]], (*this)[j][i]); - } - } - } -@@ -1106,12 +989,14 @@ namespace Dune - - MAT A(asImp()); - field_type det; -- SimdMask<typename FieldTraits<value_type>::real_type> -- nonsingularLanes(true); -- -- luDecomposition(A, ElimDet(det), nonsingularLanes, false); -- assign(det, field_type(0), !nonsingularLanes); -- -+ try -+ { -+ luDecomposition(A, ElimDet(det)); -+ } -+ catch (FMatrixError&) -+ { -+ return 0; -+ } - for (size_type i = 0; i < rows(); ++i) - det *= A[i][i]; - return det; -diff --git a/dune/common/dynmatrix.hh b/dune/common/dynmatrix.hh -index ba78118a..8c7ed8b1 100644 ---- a/dune/common/dynmatrix.hh -+++ b/dune/common/dynmatrix.hh -@@ -80,13 +80,6 @@ namespace Dune - {} - - -- template <class T, -- typename = std::enable_if_t<!Dune::IsNumber<T>::value && HasDenseMatrixAssigner<DynamicMatrix, T>::value>> -- DynamicMatrix(T const& rhs) -- { -- *this = rhs; -- } -- - //==== resize related methods - /** - * \brief resize matrix to <code>r × c</code> -@@ -108,23 +101,7 @@ namespace Dune - } - - //===== assignment -- // General assignment with resizing -- template <typename T, -- typename = std::enable_if_t<!Dune::IsNumber<T>::value>> -- DynamicMatrix& operator=(T const& rhs) { -- _data.resize(rhs.N()); -- std::fill(_data.begin(), _data.end(), row_type(rhs.M(), K(0))); -- Base::operator=(rhs); -- return *this; -- } -- -- // Specialisation: scalar assignment (no resizing) -- template <typename T, -- typename = std::enable_if_t<Dune::IsNumber<T>::value>> -- DynamicMatrix& operator=(T scalar) { -- std::fill(_data.begin(), _data.end(), scalar); -- return *this; -- } -+ using Base::operator=; - - // make this thing a matrix - size_type mat_rows() const { return _data.size(); } -diff --git a/dune/common/fmatrix.hh b/dune/common/fmatrix.hh -index d6b60d81..3291b25f 100644 ---- a/dune/common/fmatrix.hh -+++ b/dune/common/fmatrix.hh -@@ -89,6 +89,24 @@ namespace Dune - */ - FieldMatrix () {} - -+ /** \brief Constructor initializing the whole matrix with a scalar -+ */ -+ template< class Other > -+ FieldMatrix ( const Other &other ) -+ { -+ DenseMatrixAssigner< FieldMatrix< K, ROWS, COLS >, Other >::apply( *this, other ); -+ } -+ -+ // /** \brief Constructor initializing the matrix from a list of lists of scalars -+ // */ -+ // FieldMatrix (std::initializer_list<std::initializer_list<K> > const &ll) -+ // { -+ // assert(ll.size() == rows); // Actually, this is not needed any more! -+ // std::copy_n(ll.begin(), std::min(static_cast<std::size_t>(ROWS), -+ // ll.size()), -+ // _data.begin()); -+ // } -+ - /** \brief Constructor initializing the matrix from a list of vector - */ - FieldMatrix(std::initializer_list<Dune::FieldVector<K, cols> > const &l) { -@@ -98,25 +116,9 @@ namespace Dune - _data.begin()); - } - -- template <class T, -- typename = std::enable_if_t<HasDenseMatrixAssigner<FieldMatrix, T>::value>> -- FieldMatrix(T const& rhs) -- { -- *this = rhs; -- } -- -+ //===== assignment - using Base::operator=; - -- // Specialisation: FieldMatrix assignment (compile-time bounds checking) -- template <typename T, int rows, int cols> -- FieldMatrix& operator=(FieldMatrix<T,rows,cols> const &rhs) -- { -- static_assert(rows == ROWS, "Size mismatch in matrix assignment (rows)"); -- static_assert(cols == COLS, "Size mismatch in matrix assignment (columns)"); -- _data = rhs._data; -- return *this; -- } -- - //! Multiplies M from the left to this matrix, this matrix is not modified - template<int l> - FieldMatrix<K,l,cols> leftmultiplyany (const FieldMatrix<K,l,rows>& M) const -@@ -228,22 +230,19 @@ namespace Dune - */ - FieldMatrix () {} - -- /** \brief Constructor initializing the matrix from a list of vector -+ /** \brief Constructor initializing the whole matrix with a scalar - */ -- FieldMatrix(std::initializer_list<Dune::FieldVector<K, 1>> const &l) -+ FieldMatrix (const K& k) - { -- std::copy_n(l.begin(), std::min(static_cast< std::size_t >( 1 ), l.size()), &_data); -+ _data[0] = k; - } - -- template <class T, -- typename = std::enable_if_t<HasDenseMatrixAssigner<FieldMatrix, T>::value>> -- FieldMatrix(T const& rhs) -+ template< class Other > -+ FieldMatrix ( const Other &other ) - { -- *this = rhs; -+ DenseMatrixAssigner< FieldMatrix< K, 1, 1 >, Other >::apply( *this, other ); - } - -- using Base::operator=; -- - //===== solve - - //! Multiplies M from the left to this matrix, this matrix is not modified --- -2.11.0 -