diff --git a/dune/xt/common/fmatrix-2.6.hh b/dune/xt/common/fmatrix-2.6.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4f4155f92f59f9f9eaf65eb772c59ec64ea7bbf2
--- /dev/null
+++ b/dune/xt/common/fmatrix-2.6.hh
@@ -0,0 +1,1171 @@
+// This file is part of the dune-xt project:
+// https://github.com/dune-community/dune-xt
+// Copyright 2009-2018 dune-xt developers and contributors. All rights reserved.
+// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+// or GPL-2.0+ (http://opensource.org/licenses/gpl-license)
+// with "runtime exception" (http://www.dune-project.org/license.html)
+// Authors:
+// Felix Schindler (2014, 2016 - 2018)
+// René Fritze (2015 - 2016, 2018 - 2019)
+// Tim Keil (2018)
+// Tobias Leibner (2014, 2018 - 2019)
+
+#ifndef DUNE_XT_COMMON_FMATRIX_26_HH
+#define DUNE_XT_COMMON_FMATRIX_26_HH
+
+#include <initializer_list>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/xt/common/exceptions.hh>
+#include <dune/xt/common/debug.hh>
+#include <dune/xt/common/matrix.hh>
+#include <dune/xt/common/fvector.hh>
+#include <dune/xt/common/type_traits.hh>
+
+namespace Dune {
+namespace XT {
+namespace Common {
+
+
+/**
+ * \todo We need to implement all operators from the base which return the base, to rather return ourselfes!
+ */
+template <class K, int ROWS, int COLS>
+class FieldMatrix : public Dune::FieldMatrix<K, ROWS, COLS>
+{
+ typedef Dune::FieldMatrix<K, ROWS, COLS> BaseType;
+ typedef FieldMatrix<K, ROWS, COLS> ThisType;
+
+public:
+ using typename BaseType::field_type;
+ using typename BaseType::size_type;
+ using typename BaseType::value_type;
+
+ FieldMatrix(const K& kk = suitable_default<K>::value())
+ : BaseType()
+ {
+ // This is required because BaseType(kk) does not work for std::string
+ for (size_t i = 0; i < ROWS; ++i) {
+ for (size_t j = 0; j < COLS; ++j)
+ (*this)[i][j] = kk;
+ }
+ }
+
+ FieldMatrix(const size_t DXTC_DEBUG_ONLY(rr),
+ const size_t DXTC_DEBUG_ONLY(cc),
+ const K& kk = suitable_default<K>::value())
+ : BaseType()
+ {
+#ifndef NDEBUG
+ if (rr != ROWS || cc != COLS)
+ DUNE_THROW(Exceptions::wrong_input_given,
+ "You are trying to construct a FieldMatrix< ..., " << ROWS << ", " << COLS << " > (of "
+ << "static size) with " << rr << " rows and " << cc
+ << " columns!");
+#endif // NDEBUG
+ // This is required because BaseType(kk) does not work for std::string
+ for (size_t i = 0; i < ROWS; ++i) {
+ for (size_t j = 0; j < COLS; ++j)
+ (*this)[i][j] = kk;
+ }
+
+ } // ... FieldMatrix(...)
+
+ FieldMatrix(std::initializer_list<std::initializer_list<K>> list_of_rows)
+ : BaseType()
+ {
+#ifndef NDEBUG
+ if (list_of_rows.size() != ROWS)
+ DUNE_THROW(Exceptions::wrong_input_given,
+ "You are trying to construct a FieldMatrix< ..., "
+ << ROWS << ", " << COLS << " > (of "
+ << "static size) from a list modeling a matrix with " << list_of_rows.size() << " rows!");
+#endif // NDEBUG
+ size_t rr = 0;
+ for (auto row : list_of_rows) {
+#ifndef NDEBUG
+ if (row.size() != COLS)
+ DUNE_THROW(Exceptions::wrong_input_given,
+ "You are trying to construct a FieldMatrix< ..., "
+ << ROWS << ", " << COLS << " > (of "
+ << "static size) from a list with a row of length " << row.size() << "!");
+#endif // NDEBUG
+ size_t cc = 0;
+ for (auto entry : row)
+ (*this)[rr][cc++] = entry;
+ ++rr;
+ }
+ } // FieldMatrix(std::initializer_list<std::initializer_list<K>> list_of_rows)
+
+ template <class OtherMatrixType>
+ FieldMatrix(const OtherMatrixType& other, std::enable_if_t<is_matrix<OtherMatrixType>::value, int> /*dummy*/ = 0)
+ {
+ *this = other;
+ }
+
+ template <class OtherMatrixType>
+ std::enable_if_t<is_matrix<OtherMatrixType>::value, ThisType>& operator=(const OtherMatrixType& other)
+ {
+ for (size_t rr = 0; rr < ROWS; ++rr)
+ for (size_t cc = 0; cc < COLS; ++cc)
+ (*this)[rr][cc] = MatrixAbstraction<OtherMatrixType>::get_entry(other, rr, cc);
+ return *this;
+ }
+
+ Dune::XT::Common::FieldMatrix<K, COLS, ROWS> transpose() const
+ {
+ Dune::XT::Common::FieldMatrix<K, COLS, ROWS> ret;
+ for (size_t rr = 0; rr < ROWS; ++rr)
+ for (size_t cc = 0; cc < COLS; ++cc)
+ ret[cc][rr] = (*this)[rr][cc];
+ return ret;
+ }
+
+ Dune::XT::Common::FieldVector<K, ROWS> operator*(const Dune::FieldVector<K, COLS>& vec) const
+ {
+ Dune::FieldVector<K, ROWS> ret;
+ this->mv(vec, ret);
+ return ret;
+ }
+
+ //! This op is not redundant
+ ThisType operator*(const K& scal) const
+ {
+ ThisType ret(*this);
+ ret *= scal;
+ return ret;
+ }
+
+ field_type determinant() const;
+
+ template <class Func>
+ void luDecomposition(ThisType& A, Func func) const;
+
+ void invert();
+
+ template <class V, class W>
+ void solve(V& x, const W& b) const;
+
+private:
+ // copy from dune/common/densematrix.hh, we have to copy it as it is a private member of Dune::DenseMatrix
+ struct ElimPivot
+ {
+ ElimPivot(std::vector<size_type>& pivot)
+ : pivot_(pivot)
+ {
+ typedef typename std::vector<size_type>::size_type size_type;
+ for (size_type i = 0; i < pivot_.size(); ++i)
+ pivot_[i] = i;
+ }
+
+ void swap(int i, int j)
+ {
+ pivot_[i] = j;
+ }
+
+ template <typename T>
+ void operator()(const T&, int, int)
+ {}
+
+ std::vector<size_type>& pivot_;
+ }; // struct ElimPivot
+
+ template <typename V>
+ struct Elim
+ {
+ Elim(V& rhs)
+ : rhs_(&rhs)
+ {}
+
+ void swap(int i, int j)
+ {
+ std::swap((*rhs_)[i], (*rhs_)[j]);
+ }
+
+ void operator()(const typename V::field_type& factor, int k, int i)
+ {
+ (*rhs_)[k] -= factor * (*rhs_)[i];
+ }
+
+ V* rhs_;
+ };
+
+ struct ElimDet
+ {
+ ElimDet(field_type& sign)
+ : sign_(sign)
+ {
+ sign_ = 1;
+ }
+
+ void swap(int, int)
+ {
+ sign_ *= -1;
+ }
+
+ void operator()(const field_type&, int, int) {}
+
+ field_type& sign_;
+ };
+}; // class FieldMatrix<...>
+
+// Direct copy of the luDecomposition function in dune/common/densematrix.hh
+// The only (functional) change is that this version always performs pivotization (the version in dune-common only
+// performs pivotization if the diagonal entry is below a certain threshold)
+// See dune/xt/la/test/matrixinverter_for_real_matrix_from_3d_pointsource.tpl for an example where the dune-common
+// version fails due to stability issues.
+// TODO: Fixed in dune-common master (see MR !449 in dune-common's gitlab), remove this copy once we depend on a
+// suitable version of dune-common (probably 2.7).
+template <class K, int ROWS, int COLS>
+template <typename Func>
+inline void FieldMatrix<K, ROWS, COLS>::luDecomposition(FieldMatrix<K, ROWS, COLS>& A, Func func) const
+{
+ typedef typename FieldTraits<value_type>::real_type real_type;
+ real_type norm = A.infinity_norm_real(); // for relative thresholds
+ 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
+ {
+ typename FieldTraits<value_type>::real_type pivmax = fvmeta::absreal(A[i][i]);
+
+ // compute maximum of column
+ size_type imax = i;
+ typename FieldTraits<value_type>::real_type abs(0.0);
+ for (size_type k = i + 1; k < ROWS; k++)
+ if ((abs = fvmeta::absreal(A[k][i])) > pivmax) {
+ 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]);
+ assert(imax < std::numeric_limits<int>::max() && i < std::numeric_limits<int>::max());
+ func.swap(static_cast<int>(i), static_cast<int>(imax)); // swap the pivot or rhs
+ }
+
+ // singular ?
+ if (pivmax < singthres)
+ DUNE_THROW(FMatrixError, "matrix is singular");
+
+ // eliminate
+ for (size_type k = i + 1; k < ROWS; k++) {
+ field_type factor = A[k][i] / A[i][i];
+ A[k][i] = factor;
+ for (size_type j = i + 1; j < ROWS; j++)
+ A[k][j] -= factor * A[i][j];
+ assert(k < std::numeric_limits<int>::max() && i < std::numeric_limits<int>::max());
+ func(factor, static_cast<int>(k), static_cast<int>(i));
+ }
+ }
+}
+
+// Direct copy of the invert function in dune/common/densematrix.hh
+// The only (functional) change is the replacement of the luDecomposition of DenseMatrix by our own version.
+// TODO: Fixed in dune-common master (see MR !449 in dune-common's gitlab), remove this copy once we depend on a
+// suitable version of dune-common (probably 2.7).
+template <class K, int ROWS, int COLS>
+inline void FieldMatrix<K, ROWS, COLS>::invert()
+{
+ // never mind those ifs, because they get optimized away
+ if (ROWS != COLS)
+ DUNE_THROW(Dune::FMatrixError, "Can't invert a " << ROWS << "x" << COLS << " matrix!");
+ if (ROWS <= 3) {
+ BaseType::invert();
+ } else {
+ auto A = *this;
+ std::vector<size_type> pivot(ROWS);
+ this->luDecomposition(A, ElimPivot(pivot));
+ auto& L = A;
+ auto& U = A;
+
+ // initialize inverse
+ *this = field_type();
+
+ for (size_type i = 0; i < ROWS; ++i)
+ (*this)[i][i] = 1;
+
+ // L Y = I; multiple right hand sides
+ for (size_type i = 0; i < ROWS; i++)
+ for (size_type j = 0; j < i; j++)
+ for (size_type k = 0; k < ROWS; k++)
+ (*this)[i][k] -= L[i][j] * (*this)[j][k];
+
+ // U A^{-1} = Y
+ for (size_type i = ROWS; i > 0;) {
+ --i;
+ for (size_type k = 0; k < ROWS; k++) {
+ for (size_type j = i + 1; j < ROWS; j++)
+ (*this)[i][k] -= U[i][j] * (*this)[j][k];
+ (*this)[i][k] /= U[i][i];
+ }
+ }
+
+ for (size_type i = ROWS; i > 0;) {
+ --i;
+ if (i != pivot[i])
+ for (size_type j = 0; j < ROWS; ++j)
+ std::swap((*this)[j][pivot[i]], (*this)[j][i]);
+ }
+ }
+}
+
+// Direct copy of the determinant function in dune/common/densematrix.hh
+// The only (functional) change is the replacement of the luDecomposition of DenseMatrix by our own version.
+// TODO: Fixed in dune-common master (see MR !449 in dune-common's gitlab), remove this copy once we depend on a
+// suitable version of dune-common (probably 2.7).
+template <class K, int ROWS, int COLS>
+inline typename FieldMatrix<K, ROWS, COLS>::field_type FieldMatrix<K, ROWS, COLS>::determinant() const
+{
+ // never mind those ifs, because they get optimized away
+ if (ROWS != COLS)
+ DUNE_THROW(FMatrixError, "There is no determinant for a " << ROWS << "x" << COLS << " matrix!");
+
+ if (ROWS == 1)
+ return (*this)[0][0];
+
+ if (ROWS == 2)
+ return (*this)[0][0] * (*this)[1][1] - (*this)[0][1] * (*this)[1][0];
+
+ if (ROWS == 3) {
+ // code generated by maple
+ field_type t4 = (*this)[0][0] * (*this)[1][1];
+ field_type t6 = (*this)[0][0] * (*this)[1][2];
+ field_type t8 = (*this)[0][1] * (*this)[1][0];
+ field_type t10 = (*this)[0][2] * (*this)[1][0];
+ field_type t12 = (*this)[0][1] * (*this)[2][0];
+ field_type t14 = (*this)[0][2] * (*this)[2][0];
+
+ return (t4 * (*this)[2][2] - t6 * (*this)[2][1] - t8 * (*this)[2][2] + t10 * (*this)[2][1] + t12 * (*this)[1][2]
+ - t14 * (*this)[1][1]);
+ }
+
+ auto A = *this;
+ field_type det;
+ try {
+ this->luDecomposition(A, ElimDet(det));
+ } catch (FMatrixError&) {
+ return 0;
+ }
+ for (size_type i = 0; i < ROWS; ++i)
+ det *= A[i][i];
+ return det;
+}
+
+
+// Direct copy of the solve function in dune/common/densematrix.hh
+// The only (functional) change is the replacement of the luDecomposition of DenseMatrix by our own version.
+// TODO: Fixed in dune-common master (see MR !449 in dune-common's gitlab), remove this copy once we depend on a
+// suitable version of dune-common (probably 2.7).
+template <class K, int ROWS, int COLS>
+template <class V, class W>
+inline void FieldMatrix<K, ROWS, COLS>::solve(V& x, const W& b) const
+{
+ // never mind those ifs, because they get optimized away
+ if (ROWS != COLS)
+ DUNE_THROW(FMatrixError, "Can't solve for a " << ROWS << "x" << COLS << " matrix!");
+
+ if (ROWS == 1) {
+
+#ifdef DUNE_FMatrix_WITH_CHECKING
+ if (fvmeta::absreal((*this)[0][0]) < FMatrixPrecision<>::absolute_limit())
+ DUNE_THROW(FMatrixError, "matrix is singular");
+#endif
+ x[0] = b[0] / (*this)[0][0];
+
+ } else if (ROWS == 2) {
+
+ field_type detinv = (*this)[0][0] * (*this)[1][1] - (*this)[0][1] * (*this)[1][0];
+#ifdef DUNE_FMatrix_WITH_CHECKING
+ if (fvmeta::absreal(detinv) < FMatrixPrecision<>::absolute_limit())
+ DUNE_THROW(FMatrixError, "matrix is singular");
+#endif
+ detinv = 1.0 / detinv;
+
+ x[0] = detinv * ((*this)[1][1] * b[0] - (*this)[0][1] * b[1]);
+ x[1] = detinv * ((*this)[0][0] * b[1] - (*this)[1][0] * b[0]);
+
+ } else if (ROWS == 3) {
+
+ field_type d = determinant();
+#ifdef DUNE_FMatrix_WITH_CHECKING
+ if (fvmeta::absreal(d) < FMatrixPrecision<>::absolute_limit())
+ DUNE_THROW(FMatrixError, "matrix is singular");
+#endif
+
+ x[0] = (b[0] * (*this)[1][1] * (*this)[2][2] - b[0] * (*this)[2][1] * (*this)[1][2]
+ - b[1] * (*this)[0][1] * (*this)[2][2] + b[1] * (*this)[2][1] * (*this)[0][2]
+ + b[2] * (*this)[0][1] * (*this)[1][2] - b[2] * (*this)[1][1] * (*this)[0][2])
+ / d;
+
+ x[1] = ((*this)[0][0] * b[1] * (*this)[2][2] - (*this)[0][0] * b[2] * (*this)[1][2]
+ - (*this)[1][0] * b[0] * (*this)[2][2] + (*this)[1][0] * b[2] * (*this)[0][2]
+ + (*this)[2][0] * b[0] * (*this)[1][2] - (*this)[2][0] * b[1] * (*this)[0][2])
+ / d;
+
+ x[2] = ((*this)[0][0] * (*this)[1][1] * b[2] - (*this)[0][0] * (*this)[2][1] * b[1]
+ - (*this)[1][0] * (*this)[0][1] * b[2] + (*this)[1][0] * (*this)[2][1] * b[0]
+ + (*this)[2][0] * (*this)[0][1] * b[1] - (*this)[2][0] * (*this)[1][1] * b[0])
+ / d;
+
+ } else {
+
+ V& rhs = x; // use x to store rhs
+ rhs = b; // copy data
+ Elim<V> elim(rhs);
+ auto A = *this;
+
+ this->luDecomposition(A, elim);
+
+ // backsolve
+ for (int i = ROWS - 1; i >= 0; i--) {
+ for (size_type j = i + 1; j < ROWS; j++)
+ rhs[i] -= A[i][j] * x[j];
+ x[i] = rhs[i] / A[i][i];
+ }
+ }
+}
+
+
+/**
+ * \todo We need to implement all operators from the base which return the base, to rather return ourselfes!
+ */
+template <class K>
+class FieldMatrix<K, 1, 1> : public Dune::FieldMatrix<K, 1, 1>
+{
+ static const int ROWS = 1;
+ static const int COLS = 1;
+ typedef Dune::FieldMatrix<K, ROWS, COLS> BaseType;
+ typedef FieldMatrix<K, ROWS, COLS> ThisType;
+
+public:
+ FieldMatrix(const K& kk = suitable_default<K>::value())
+ : BaseType()
+ {
+ (*this)[0][0] = kk;
+ }
+
+ FieldMatrix(std::initializer_list<std::initializer_list<K>> list_of_rows)
+ : BaseType()
+ {
+#ifndef NDEBUG
+ if (list_of_rows.size() != ROWS)
+ DUNE_THROW(Exceptions::wrong_input_given,
+ "You are trying to construct a FieldMatrix< ..., 1, 1 > (of "
+ << "static size) from a list modeling a matrix with " << list_of_rows.size() << " rows!");
+#endif // NDEBUG
+ for (auto row : list_of_rows) {
+#ifndef NDEBUG
+ if (row.size() != COLS)
+ DUNE_THROW(Exceptions::wrong_input_given,
+ "You are trying to construct a FieldMatrix< ..., 1, 1 > (of "
+ << "static size) from a list with a row of length " << row.size() << "!");
+#endif // NDEBUG
+ for (auto entry : row)
+ (*this)[0][0] = entry;
+ }
+ } // FieldMatrix(std::initializer_list<std::initializer_list<K>> list_of_rows)
+
+ FieldMatrix(const size_t DXTC_DEBUG_ONLY(rr),
+ const size_t DXTC_DEBUG_ONLY(cc),
+ const K& kk = suitable_default<K>::value())
+ : BaseType()
+ {
+#ifndef NDEBUG
+ if (rr != ROWS || cc != COLS)
+ DUNE_THROW(Exceptions::wrong_input_given,
+ "You are trying to construct a FieldMatrix< ..., " << ROWS << ", " << COLS << " > (of "
+ << "static size) with " << rr << " rows and " << cc
+ << " columns!");
+#endif // NDEBUG
+ (*this)[0][0] = kk;
+ } // ... FieldMatrix(...)
+
+ FieldMatrix(const BaseType& other)
+ : BaseType(other)
+ {}
+
+ FieldMatrix(const Dune::XT::Common::FieldVector<K, 1>& other)
+ : BaseType()
+ {
+ (*this)[0][0] = other[0];
+ }
+
+ FieldMatrix(const Dune::FieldVector<K, 1>& other)
+ : BaseType()
+ {
+ (*this)[0][0] = other[0];
+ }
+
+ Dune::XT::Common::FieldMatrix<K, COLS, ROWS> transpose() const
+ {
+ return *this;
+ }
+
+ using BaseType::operator=;
+
+ ThisType& operator=(const FieldVector<K, 1>& other)
+ {
+ BaseType::operator=(other[0]);
+ return *this;
+ }
+
+ ThisType operator*(const K& scal) const
+ {
+ ThisType ret(*this);
+ ret *= scal;
+ return ret;
+ }
+}; // class FieldMatrix
+
+
+template <class K, size_t num_blocks, size_t block_rows, size_t block_cols = block_rows>
+class BlockedFieldMatrix
+{
+ using ThisType = BlockedFieldMatrix;
+
+public:
+ static constexpr size_t num_rows = num_blocks * block_rows;
+ static constexpr size_t num_cols = num_blocks * block_cols;
+ using MatrixType = Dune::FieldMatrix<K, num_rows, num_cols>;
+ using BlockType = FieldMatrix<K, block_rows, block_cols>;
+
+ BlockedFieldMatrix(const K& val = K(0.))
+ : backend_(BlockType(val))
+ {}
+
+ BlockedFieldMatrix(const size_t DXTC_DEBUG_ONLY(rows), const size_t DXTC_DEBUG_ONLY(cols), const K& val = K(0.))
+ : backend_(BlockType(val))
+ {
+ assert(rows == num_rows && cols == num_cols && "Requested shape has to match static shape!");
+ }
+
+ template <class OtherMatrixType>
+ BlockedFieldMatrix(const OtherMatrixType& other,
+ std::enable_if_t<is_matrix<OtherMatrixType>::value, int> /*dummy*/ = 0)
+ {
+ *this = other;
+ }
+
+ template <class OtherMatrixType>
+ std::enable_if_t<is_matrix<OtherMatrixType>::value, ThisType>& operator=(const OtherMatrixType& other)
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ for (size_t ll = 0; ll < block_rows; ++ll)
+ for (size_t mm = 0; mm < block_cols; ++mm)
+ backend_[jj][ll][mm] =
+ MatrixAbstraction<OtherMatrixType>::get_entry(other, jj * block_rows + ll, jj * block_cols + mm);
+ return *this;
+ }
+
+
+ BlockedFieldMatrix(const MatrixType& other)
+ {
+ *this = other;
+ }
+
+ BlockedFieldMatrix(const BlockType& block)
+ : backend_(block)
+ {}
+
+ bool operator==(const ThisType& other) const
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ if (block(jj) != other.block(jj))
+ return false;
+ return true;
+ }
+
+ K get_entry(const size_t ii, const size_t jj) const
+ {
+ assert(ii < num_rows && jj < num_cols);
+ if (!in_pattern(ii, jj))
+ return K(0.);
+ return backend_[ii / block_rows][ii % block_rows][jj % block_cols];
+ }
+
+ K& get_entry(const size_t jj, const size_t ll, const size_t mm)
+ {
+ assert(is_valid_entry(jj, ll, mm));
+ return backend_[jj][ll][mm];
+ }
+
+ const K& get_entry(const size_t jj, const size_t ll, const size_t mm) const
+ {
+ assert(is_valid_entry(jj, ll, mm));
+ return backend_[jj][ll][mm];
+ }
+
+ void set_entry(const size_t ii, const size_t jj, const K& val)
+ {
+ assert(ii < num_rows && jj < num_cols);
+ if (in_pattern(ii, jj))
+ backend_[ii / block_rows][ii % block_rows][jj % block_cols] = val;
+ else if (XT::Common::FloatCmp::ne(val, K(0)))
+ DUNE_THROW(Dune::MathError, "Tried to modify a value that is not in the pattern!");
+ }
+
+ void set_entry(const size_t jj, const size_t ll, const size_t mm, const K& val)
+ {
+ assert(is_valid_entry(jj, ll, mm));
+ backend_[jj][ll][mm] = val;
+ }
+
+ void add_to_entry(const size_t ii, const size_t jj, const K& val)
+ {
+ assert(ii < num_rows && jj < num_cols);
+ if (!in_pattern(ii, jj) && XT::Common::FloatCmp::ne(val, K(0)))
+ DUNE_THROW(Dune::MathError, "Tried to modify a value that is not in the pattern!");
+ backend_[ii / block_rows][ii % block_rows][jj % block_cols] += val;
+ }
+
+ void add_to_entry(const size_t jj, const size_t ll, const size_t mm, const K& val)
+ {
+ assert(is_valid_entry(jj, ll, mm));
+ backend_[jj][ll][mm] += val;
+ }
+
+ BlockType& block(const size_t jj)
+ {
+ assert(jj < num_blocks);
+ return backend_[jj];
+ }
+
+ const BlockType& block(const size_t jj) const
+ {
+ assert(jj < num_blocks);
+ return backend_[jj];
+ }
+
+ void mv(const Dune::FieldVector<K, num_cols>& x, Dune::FieldVector<K, num_rows>& ret) const
+ {
+ std::fill(ret.begin(), ret.end(), 0.);
+ for (size_t jj = 0; jj < num_blocks; ++jj) {
+ const auto row_offset = block_rows * jj;
+ const auto col_offset = block_cols * jj;
+ for (size_t ll = 0; ll < block_rows; ++ll)
+ for (size_t mm = 0; mm < block_cols; ++mm)
+ ret[row_offset + ll] += backend_[jj][ll][mm] * x[col_offset + mm];
+ } // jj
+ } // void mv(...)
+
+ void mv(const BlockedFieldVector<K, num_blocks, block_cols>& x,
+ BlockedFieldVector<K, num_blocks, block_rows>& ret) const
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ backend_[jj].mv(x.block(jj), ret.block(jj));
+ } // void mv(...)
+
+ void mtv(const Dune::FieldVector<K, num_rows>& x, Dune::FieldVector<K, num_cols>& ret) const
+ {
+ std::fill(ret.begin(), ret.end(), 0.);
+ for (size_t jj = 0; jj < num_blocks; ++jj) {
+ const auto row_offset = block_rows * jj;
+ const auto col_offset = block_cols * jj;
+ for (size_t mm = 0; mm < block_cols; ++mm)
+ for (size_t ll = 0; ll < block_rows; ++ll)
+ ret[col_offset + mm] += backend_[jj][ll][mm] * x[row_offset + ll];
+ } // jj
+ } // void mtv(...)
+
+ void mtv(const BlockedFieldVector<K, num_blocks, block_rows>& x,
+ BlockedFieldVector<K, num_blocks, block_cols>& ret) const
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ backend_[jj].mtv(x.block(jj), ret.block(jj));
+ } // void mv(...)
+
+ template <size_t br, size_t bc>
+ ThisType& rightmultiply(const BlockedFieldMatrix<K, num_blocks, br, bc>& other)
+ {
+ assert((this != &other) && "Multiplying a matrix by itself gives wrong results, please copy before!");
+ static_assert(br == bc, "Cannot rightmultiply with non-square matrix");
+ static_assert(br == block_cols, "Size mismatch");
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ backend_[jj].rightmultiply(other.backend_[jj]);
+ return *this;
+ }
+
+ BlockedFieldMatrix<K, num_blocks, block_cols, block_rows> transpose()
+ {
+ BlockedFieldMatrix<K, num_blocks, block_cols, block_rows> ret;
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ ret.block(jj) = block(jj).transpose();
+ return ret;
+ }
+
+ ThisType operator*(const ThisType& other) const
+ {
+ ThisType ret(*this);
+ ret.rightmultiply(other);
+ return ret;
+ }
+
+ ThisType& operator*=(const K& val)
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ block(jj) *= val;
+ return *this;
+ }
+
+ ThisType& operator+=(const ThisType& other)
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ block(jj) += other.block(jj);
+ return *this;
+ }
+
+ ThisType operator+(const ThisType& other) const
+ {
+ ThisType ret(*this);
+ ret += other;
+ return ret;
+ }
+
+ ThisType& operator-=(const ThisType& other)
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ block(jj) -= other.block(jj);
+ return *this;
+ }
+
+ ThisType operator-(const ThisType& other) const
+ {
+ ThisType ret(*this);
+ ret -= other;
+ return ret;
+ }
+
+ static bool in_pattern(const size_t ii, const size_t jj)
+ {
+ return (ii / block_rows == jj / block_cols);
+ }
+
+ static bool is_valid_entry(const size_t jj, const size_t ll, const size_t mm)
+ {
+ return (jj < num_blocks && ll < block_rows && mm < block_cols);
+ }
+
+ DynamicMatrix<K> convert_to_dynamic_matrix() const
+ {
+ DynamicMatrix<K> ret(num_rows, num_cols, 0.);
+ for (size_t jj = 0; jj < num_blocks; ++jj) {
+ const size_t row_offset = jj * block_rows;
+ const size_t col_offset = jj * block_cols;
+ for (size_t rr = 0; rr < block_rows; ++rr)
+ for (size_t cc = 0; cc < block_cols; ++cc)
+ ret[row_offset + rr][col_offset + cc] = block(jj)[rr][cc];
+ } // jj
+ return ret;
+ }
+
+ template <class CharType, class CharTraits>
+ friend std::basic_ostream<CharType, CharTraits>& operator<<(std::basic_ostream<CharType, CharTraits>& out,
+ const ThisType& mat)
+ {
+ return output_matrix(out, mat);
+ } // ... operator<<(...)
+
+private:
+ FieldVector<BlockType, num_blocks> backend_;
+};
+
+
+template <class K, int N, int M>
+struct MatrixAbstraction<Dune::XT::Common::FieldMatrix<K, N, M>>
+{
+ typedef Dune::XT::Common::FieldMatrix<K, N, M> MatrixType;
+ typedef typename Dune::FieldTraits<K>::field_type ScalarType;
+ typedef typename Dune::FieldTraits<K>::real_type RealType;
+ typedef ScalarType S;
+ typedef RealType R;
+ template <size_t rows = N, size_t cols = M, class FieldType = K>
+ using MatrixTypeTemplate = Dune::XT::Common::FieldMatrix<FieldType, rows, cols>;
+
+ static const bool is_matrix = true;
+
+ static const bool has_static_size = true;
+
+ static const size_t static_rows = N;
+
+ static const size_t static_cols = M;
+
+ static const constexpr StorageLayout storage_layout = StorageLayout::dense_row_major;
+
+ static constexpr bool has_ostream = true;
+
+ template <class SparsityPatternType = FullPattern>
+ static inline MatrixType create(const size_t rows,
+ const size_t cols,
+ const ScalarType& val = suitable_default<ScalarType>::value(),
+ const SparsityPatternType& /*pattern*/ = SparsityPatternType())
+ {
+ return MatrixType(rows, cols, val);
+ }
+
+ template <class SparsityPatternType = FullPattern>
+ static inline std::unique_ptr<MatrixType> make_unique(const size_t rows,
+ const size_t cols,
+ const ScalarType& val = suitable_default<ScalarType>::value(),
+ const SparsityPatternType& /*pattern*/ = SparsityPatternType())
+ {
+ return std::make_unique<MatrixType>(rows, cols, val);
+ }
+
+ static constexpr size_t rows(const MatrixType& /*mat*/)
+ {
+ return N;
+ }
+
+ static constexpr size_t cols(const MatrixType& /*mat*/)
+ {
+ return M;
+ }
+
+ static inline void set_entry(MatrixType& mat, const size_t row, const size_t col, const ScalarType& val)
+ {
+ mat[row][col] = val;
+ }
+
+ static inline ScalarType get_entry(const MatrixType& mat, const size_t row, const size_t col)
+ {
+ return mat[row][col];
+ }
+
+ static inline void add_to_entry(MatrixType& mat, const size_t row, const size_t col, const ScalarType& val)
+ {
+ mat[row][col] += val;
+ }
+
+
+ static inline ScalarType* data(MatrixType& mat)
+ {
+ return &(mat[0][0]);
+ }
+
+ static inline const ScalarType* data(const MatrixType& mat)
+ {
+ return &(mat[0][0]);
+ }
+};
+
+template <class K, size_t num_blocks, size_t block_rows, size_t block_cols>
+struct MatrixAbstraction<Dune::XT::Common::BlockedFieldMatrix<K, num_blocks, block_rows, block_cols>>
+{
+ typedef Dune::XT::Common::BlockedFieldMatrix<K, num_blocks, block_rows, block_cols> MatrixType;
+ typedef typename Dune::FieldTraits<K>::field_type ScalarType;
+ typedef typename Dune::FieldTraits<K>::real_type RealType;
+ typedef ScalarType S;
+ typedef RealType R;
+
+ static const bool is_matrix = true;
+
+ static const bool has_static_size = true;
+
+ static const size_t static_rows = MatrixType::num_rows;
+
+ static const size_t static_cols = MatrixType::num_cols;
+
+ template <size_t rows = static_rows, size_t cols = static_cols, class FieldType = K>
+ using MatrixTypeTemplate = Dune::XT::Common::BlockedFieldMatrix<FieldType, rows / block_rows, block_rows, block_cols>;
+
+ static const constexpr StorageLayout storage_layout = StorageLayout::other;
+
+ static constexpr bool has_ostream = true;
+
+ template <class SparsityPatternType = FullPattern>
+ static inline MatrixType create(const size_t rows,
+ const size_t cols,
+ const ScalarType& val = suitable_default<ScalarType>::value(),
+ const SparsityPatternType& /*pattern*/ = SparsityPatternType())
+ {
+ return MatrixType(rows, cols, val);
+ }
+
+ template <class SparsityPatternType = FullPattern>
+ static inline std::unique_ptr<MatrixType> make_unique(const size_t rows,
+ const size_t cols,
+ const ScalarType& val = suitable_default<ScalarType>::value(),
+ const SparsityPatternType& /*pattern*/ = SparsityPatternType())
+ {
+ return std::make_unique<MatrixType>(rows, cols, val);
+ }
+
+ static constexpr size_t rows(const MatrixType& /*mat*/)
+ {
+ return static_rows;
+ }
+
+ static constexpr size_t cols(const MatrixType& /*mat*/)
+ {
+ return static_cols;
+ }
+
+ static inline void set_entry(MatrixType& mat, const size_t row, const size_t col, const ScalarType& val)
+ {
+ mat.set_entry(row, col, val);
+ }
+
+ static inline ScalarType get_entry(const MatrixType& mat, const size_t row, const size_t col)
+ {
+ return mat.get_entry(row, col);
+ }
+
+ static inline void add_to_entry(MatrixType& mat, const size_t row, const size_t col, const ScalarType& val)
+ {
+ mat.add_to_entry(row, col, val);
+ }
+
+ static inline ScalarType* data(MatrixType& /*mat*/)
+ {
+ DUNE_THROW(InvalidStateException, "Do not call me if storage layout is not dense!");
+ return nullptr;
+ }
+
+ static inline const ScalarType* data(const MatrixType& /*mat*/)
+ {
+ DUNE_THROW(InvalidStateException, "Do not call me if storage_layout is not dense!");
+ return nullptr;
+ }
+};
+
+
+template <class M>
+typename std::enable_if<is_matrix<M>::value && MatrixAbstraction<M>::has_static_size,
+ std::unique_ptr<FieldMatrix<typename MatrixAbstraction<M>::S,
+ MatrixAbstraction<M>::static_rows,
+ MatrixAbstraction<M>::static_cols>>>::type
+make_field_container_ptr(const M& mat)
+{
+ static const size_t rows = MatrixAbstraction<M>::static_rows;
+ static const size_t cols = MatrixAbstraction<M>::static_cols;
+ auto ret = std::make_unique<FieldMatrix<typename MatrixAbstraction<M>::S, rows, cols>>;
+ for (size_t ii = 0; ii < rows; ++ii)
+ for (size_t jj = 0; jj < cols; ++jj)
+ (*ret)[ii][jj] = get_matrix_entry(mat, ii, jj);
+ return std::move(ret);
+}
+
+
+template <class M>
+typename std::enable_if<is_matrix<M>::value && MatrixAbstraction<M>::has_static_size,
+ FieldMatrix<typename MatrixAbstraction<M>::S,
+ MatrixAbstraction<M>::static_rows,
+ MatrixAbstraction<M>::static_cols>>::type
+make_field_container(const M& mat)
+{
+ static const size_t rows = MatrixAbstraction<M>::static_rows;
+ static const size_t cols = MatrixAbstraction<M>::static_cols;
+ FieldMatrix<typename MatrixAbstraction<M>::S, rows, cols> ret;
+ for (size_t ii = 0; ii < rows; ++ii)
+ for (size_t jj = 0; jj < cols; ++jj)
+ ret[ii][jj] = get_matrix_entry(mat, ii, jj);
+ return ret;
+}
+
+
+template <class K, int ROWS, int COLS>
+FieldMatrix<K, ROWS, COLS> make_field_container(Dune::FieldMatrix<K, ROWS, COLS>&& vec)
+{
+ return std::move(vec);
+}
+
+
+template <class K, int ROWS, int COLS>
+FieldMatrix<K, ROWS, COLS> make_field_container(FieldMatrix<K, ROWS, COLS>&& vec)
+{
+ return std::move(vec);
+}
+
+
+template <class K, int ROWS, int COLS>
+typename std::enable_if<is_arithmetic<K>::value && !is_complex<K>::value, FieldMatrix<K, ROWS, COLS>>::type
+real(const Dune::FieldMatrix<K, ROWS, COLS>& real_mat)
+{
+ return real_mat;
+}
+
+template <class K, int ROWS, int COLS>
+typename std::enable_if<is_arithmetic<K>::value && !is_complex<K>::value, FieldMatrix<K, ROWS, COLS>>::type
+real(const FieldMatrix<K, ROWS, COLS>& real_mat)
+{
+ return real_mat;
+}
+
+template <class K, int ROWS, int COLS>
+typename std::enable_if<is_arithmetic<K>::value && !is_complex<K>::value, FieldMatrix<K, ROWS, COLS>>::type
+real(Dune::FieldMatrix<K, ROWS, COLS>&& real_mat)
+{
+ return std::move(real_mat);
+}
+
+template <class K, int ROWS, int COLS>
+typename std::enable_if<is_arithmetic<K>::value && !is_complex<K>::value, FieldMatrix<K, ROWS, COLS>>::type
+real(FieldMatrix<K, ROWS, COLS>&& real_mat)
+{
+ return std::move(real_mat);
+}
+
+template <class K, int ROWS, int COLS>
+typename std::enable_if<is_complex<K>::value, FieldMatrix<real_t<K>, ROWS, COLS>>::type
+real(const Dune::FieldMatrix<K, ROWS, COLS>& complex_mat)
+{
+ FieldMatrix<real_t<K>, ROWS, COLS> real_mat;
+ for (size_t ii = 0; ii < ROWS; ++ii)
+ for (size_t jj = 0; jj < COLS; ++jj)
+ real_mat[ii][jj] = complex_mat[ii][jj].real();
+ return real_mat;
+}
+
+template <class K, int ROWS, int COLS>
+typename std::enable_if<is_complex<K>::value, FieldMatrix<real_t<K>, ROWS, COLS>>::type
+real(const FieldMatrix<K, ROWS, COLS>& complex_mat)
+{
+ FieldMatrix<real_t<K>, ROWS, COLS> real_mat;
+ for (size_t ii = 0; ii < ROWS; ++ii)
+ for (size_t jj = 0; jj < COLS; ++jj)
+ real_mat[ii][jj] = complex_mat[ii][jj].real();
+ return real_mat;
+}
+
+
+template <class K, int ROWS, int COLS>
+typename std::enable_if<is_arithmetic<K>::value && !is_complex<K>::value, FieldMatrix<K, ROWS, COLS>>::type
+imag(const Dune::FieldMatrix<K, ROWS, COLS>& /*real_mat*/)
+{
+ return FieldMatrix<K, ROWS, COLS>(0);
+}
+
+template <class K, int ROWS, int COLS>
+typename std::enable_if<is_arithmetic<K>::value && !is_complex<K>::value, FieldMatrix<K, ROWS, COLS>>::type
+imag(const FieldMatrix<K, ROWS, COLS>& /*real_mat*/)
+{
+ return FieldMatrix<K, ROWS, COLS>(0);
+}
+
+template <class K, int ROWS, int COLS>
+typename std::enable_if<is_complex<K>::value, FieldMatrix<real_t<K>, ROWS, COLS>>::type
+imag(const Dune::FieldMatrix<K, ROWS, COLS>& complex_mat)
+{
+ FieldMatrix<real_t<K>, ROWS, COLS> real_mat;
+ for (size_t ii = 0; ii < ROWS; ++ii)
+ for (size_t jj = 0; jj < COLS; ++jj)
+ real_mat[ii][jj] = complex_mat[ii][jj].imag();
+ return real_mat;
+}
+
+template <class K, int ROWS, int COLS>
+typename std::enable_if<is_complex<K>::value, FieldMatrix<real_t<K>, ROWS, COLS>>::type
+imag(const FieldMatrix<K, ROWS, COLS>& complex_mat)
+{
+ FieldMatrix<real_t<K>, ROWS, COLS> real_mat;
+ for (size_t ii = 0; ii < ROWS; ++ii)
+ for (size_t jj = 0; jj < COLS; ++jj)
+ real_mat[ii][jj] = complex_mat[ii][jj].imag();
+ return real_mat;
+}
+
+
+} // namespace Common
+} // namespace XT
+
+
+template <class L, int ROWS, int COLS, class R>
+Dune::XT::Common::FieldMatrix<typename PromotionTraits<L, R>::PromotedType, ROWS, COLS>
+operator-(const Dune::FieldMatrix<L, ROWS, COLS>& left, const Dune::FieldMatrix<R, ROWS, COLS>& right)
+{
+ Dune::XT::Common::FieldMatrix<typename PromotionTraits<L, R>::PromotedType, ROWS, COLS> ret = left;
+ ret -= right;
+ return ret;
+}
+
+
+template <class L, int ROWS, int COLS, class R>
+Dune::XT::Common::FieldMatrix<typename PromotionTraits<L, R>::PromotedType, ROWS, COLS>
+operator+(const Dune::FieldMatrix<L, ROWS, COLS>& left, const Dune::FieldMatrix<R, ROWS, COLS>& right)
+{
+ Dune::XT::Common::FieldMatrix<typename PromotionTraits<L, R>::PromotedType, ROWS, COLS> ret = left;
+ ret += right;
+ return ret;
+}
+
+
+template <class K, int L_ROWS, int L_COLS, int R_COLS>
+Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS> operator*(const Dune::FieldMatrix<K, L_ROWS, L_COLS>& left,
+ const Dune::FieldMatrix<K, L_COLS, R_COLS>& right)
+{
+ return left.rightmultiplyany(right);
+}
+
+// we need this explicit overload to fix an ambiguous operator* error due to the automatic conversion from
+// FieldMatrix<K, 1, 1> to K
+template <class K, int L_ROWS>
+Dune::XT::Common::FieldMatrix<K, L_ROWS, 1> operator*(const Dune::XT::Common::FieldMatrix<K, L_ROWS, 1>& left,
+ const Dune::FieldMatrix<K, 1, 1>& right)
+{
+ return left.rightmultiplyany(right);
+}
+
+
+template <class K, int L_ROWS, int L_COLS, int R_COLS>
+void rightmultiply(Dune::FieldMatrix<K, L_ROWS, R_COLS>& ret,
+ const Dune::FieldMatrix<K, L_ROWS, L_COLS>& left,
+ const Dune::FieldMatrix<K, L_COLS, R_COLS>& right)
+{
+ for (size_t ii = 0; ii < L_ROWS; ++ii) {
+ for (size_t jj = 0; jj < R_COLS; ++jj) {
+ ret[ii][jj] = 0.;
+ for (size_t kk = 0; kk < L_COLS; ++kk)
+ ret[ii][jj] += left[ii][kk] * right[kk][jj];
+ }
+ }
+}
+
+template <class L, int L_ROWS, int L_COLS, class R, int R_COLS>
+typename std::enable_if<
+ !std::is_same<L, R>::value,
+ Dune::XT::Common::FieldMatrix<typename PromotionTraits<L, R>::PromotedType, L_ROWS, R_COLS>>::type
+operator*(const Dune::FieldMatrix<L, L_ROWS, L_COLS>& left, const Dune::FieldMatrix<R, L_COLS, R_COLS>& right)
+{
+ using Promoted = Dune::XT::Common::FieldMatrix<typename PromotionTraits<L, R>::PromotedType, L_ROWS, R_COLS>;
+ using Dune::XT::Common::convert_to;
+ return convert_to<Promoted>(left).rightmultiplyany(convert_to<Promoted>(right));
+}
+
+// versions that do not allocate matrices on the stack (for large matrices)
+template <class K, int L_ROWS, int L_COLS, int R_COLS>
+std::unique_ptr<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>
+operator*(const std::unique_ptr<Dune::FieldMatrix<K, L_ROWS, L_COLS>>& left,
+ const Dune::FieldMatrix<K, L_COLS, R_COLS>& right)
+{
+ auto ret = std::make_unique<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>();
+ rightmultiply(*ret, *left, right);
+ return ret;
+}
+
+template <class K, int L_ROWS, int L_COLS, int R_COLS>
+std::unique_ptr<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>
+operator*(const Dune::FieldMatrix<K, L_ROWS, L_COLS>& left,
+ const std::unique_ptr<Dune::FieldMatrix<K, L_COLS, R_COLS>>& right)
+{
+ auto ret = std::make_unique<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>();
+ rightmultiply(*ret, left, *right);
+ return ret;
+}
+
+template <class K, int L_ROWS, int L_COLS, int R_COLS>
+std::unique_ptr<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>
+operator*(const std::unique_ptr<Dune::FieldMatrix<K, L_ROWS, L_COLS>>& left,
+ const std::unique_ptr<Dune::FieldMatrix<K, L_COLS, R_COLS>>& right)
+{
+ return left * *right;
+}
+
+
+} // namespace Dune
+
+#endif // DUNE_XT_COMMON_FMATRIX_26_HH
diff --git a/dune/xt/common/fmatrix-2.7.hh b/dune/xt/common/fmatrix-2.7.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e495318c21fc58fa50e6f19ec6e65be41db4830d
--- /dev/null
+++ b/dune/xt/common/fmatrix-2.7.hh
@@ -0,0 +1,1008 @@
+// This file is part of the dune-xt project:
+// https://github.com/dune-community/dune-xt
+// Copyright 2009-2018 dune-xt developers and contributors. All rights reserved.
+// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+// or GPL-2.0+ (http://opensource.org/licenses/gpl-license)
+// with "runtime exception" (http://www.dune-project.org/license.html)
+// Authors:
+// Felix Schindler (2014, 2016 - 2018)
+// René Fritze (2015 - 2016, 2018 - 2019)
+// Tim Keil (2018)
+// Tobias Leibner (2014, 2018 - 2019)
+
+#ifndef DUNE_XT_COMMON_FMATRIX_27_HH
+#define DUNE_XT_COMMON_FMATRIX_27_HH
+
+#include <initializer_list>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/xt/common/exceptions.hh>
+#include <dune/xt/common/debug.hh>
+#include <dune/xt/common/matrix.hh>
+#include <dune/xt/common/fvector.hh>
+#include <dune/xt/common/type_traits.hh>
+
+namespace Dune {
+namespace XT {
+namespace Common {
+
+
+/**
+ * \todo We need to implement all operators from the base which return the base, to rather return ourselfes!
+ */
+template <class K, int ROWS, int COLS>
+class FieldMatrix : public Dune::FieldMatrix<K, ROWS, COLS>
+{
+ using BaseType = Dune::FieldMatrix<K, ROWS, COLS>;
+ using ThisType = FieldMatrix;
+
+public:
+ using typename BaseType::field_type;
+ using typename BaseType::size_type;
+ using typename BaseType::value_type;
+
+ FieldMatrix(const K& kk = suitable_default<K>::value())
+ : BaseType()
+ {
+ // This is required because BaseType(kk) does not work for std::string
+ for (size_t i = 0; i < ROWS; ++i) {
+ for (size_t j = 0; j < COLS; ++j)
+ (*this)[i][j] = kk;
+ }
+ }
+
+ FieldMatrix(const size_t DXTC_DEBUG_ONLY(rr),
+ const size_t DXTC_DEBUG_ONLY(cc),
+ const K& kk = suitable_default<K>::value())
+ : BaseType()
+ {
+#ifndef NDEBUG
+ if (rr != ROWS || cc != COLS)
+ DUNE_THROW(Exceptions::wrong_input_given,
+ "You are trying to construct a FieldMatrix< ..., " << ROWS << ", " << COLS << " > (of "
+ << "static size) with " << rr << " rows and " << cc
+ << " columns!");
+#endif // NDEBUG
+ // This is required because BaseType(kk) does not work for std::string
+ for (size_t i = 0; i < ROWS; ++i) {
+ for (size_t j = 0; j < COLS; ++j)
+ (*this)[i][j] = kk;
+ }
+
+ } // ... FieldMatrix(...)
+
+ FieldMatrix(std::initializer_list<std::initializer_list<K>> list_of_rows)
+ : BaseType()
+ {
+#ifndef NDEBUG
+ if (list_of_rows.size() != ROWS)
+ DUNE_THROW(Exceptions::wrong_input_given,
+ "You are trying to construct a FieldMatrix< ..., "
+ << ROWS << ", " << COLS << " > (of "
+ << "static size) from a list modeling a matrix with " << list_of_rows.size() << " rows!");
+#endif // NDEBUG
+ size_t rr = 0;
+ for (auto row : list_of_rows) {
+#ifndef NDEBUG
+ if (row.size() != COLS)
+ DUNE_THROW(Exceptions::wrong_input_given,
+ "You are trying to construct a FieldMatrix< ..., "
+ << ROWS << ", " << COLS << " > (of "
+ << "static size) from a list with a row of length " << row.size() << "!");
+#endif // NDEBUG
+ size_t cc = 0;
+ for (auto entry : row)
+ (*this)[rr][cc++] = entry;
+ ++rr;
+ }
+ } // FieldMatrix(std::initializer_list<std::initializer_list<K>> list_of_rows)
+
+ template <class OtherMatrixType>
+ FieldMatrix(const OtherMatrixType& other, std::enable_if_t<is_matrix<OtherMatrixType>::value, int> /*dummy*/ = 0)
+ {
+ *this = other;
+ }
+
+ template <class OtherMatrixType>
+ std::enable_if_t<is_matrix<OtherMatrixType>::value, ThisType>& operator=(const OtherMatrixType& other)
+ {
+ for (size_t rr = 0; rr < ROWS; ++rr)
+ for (size_t cc = 0; cc < COLS; ++cc)
+ (*this)[rr][cc] = MatrixAbstraction<OtherMatrixType>::get_entry(other, rr, cc);
+ return *this;
+ }
+
+ Dune::XT::Common::FieldMatrix<K, COLS, ROWS> transpose() const
+ {
+ Dune::XT::Common::FieldMatrix<K, COLS, ROWS> ret;
+ for (size_t rr = 0; rr < ROWS; ++rr)
+ for (size_t cc = 0; cc < COLS; ++cc)
+ ret[cc][rr] = (*this)[rr][cc];
+ return ret;
+ }
+
+ Dune::XT::Common::FieldVector<K, ROWS> operator*(const Dune::FieldVector<K, COLS>& vec) const
+ {
+ Dune::FieldVector<K, ROWS> ret;
+ this->mv(vec, ret);
+ return ret;
+ }
+
+ //! vector space addition -- two-argument version
+ template <class OtherScalar>
+ friend auto operator+(const ThisType& matrixA, const Dune::FieldMatrix<OtherScalar, ROWS, COLS>& matrixB)
+ {
+ FieldMatrix<typename PromotionTraits<K, OtherScalar>::PromotedType, ROWS, COLS> result;
+ for (size_type i = 0; i < ROWS; ++i)
+ for (size_type j = 0; j < COLS; ++j)
+ result[i][j] = matrixA[i][j] + matrixB[i][j];
+ return result;
+ }
+
+ //! vector space subtraction -- two-argument version
+ template <class OtherScalar>
+ friend auto operator-(const ThisType& matrixA, const Dune::FieldMatrix<OtherScalar, ROWS, COLS>& matrixB)
+ {
+ FieldMatrix<typename PromotionTraits<K, OtherScalar>::PromotedType, ROWS, COLS> result;
+ for (size_type i = 0; i < ROWS; ++i)
+ for (size_type j = 0; j < COLS; ++j)
+ result[i][j] = matrixA[i][j] - matrixB[i][j];
+ return result;
+ }
+
+ //! vector space multiplication with scalar
+ template <class Scalar, std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
+ friend auto operator*(const ThisType& matrix, Scalar scalar)
+ {
+ FieldMatrix<typename PromotionTraits<K, Scalar>::PromotedType, ROWS, COLS> result;
+ for (size_type i = 0; i < ROWS; ++i)
+ for (size_type j = 0; j < COLS; ++j)
+ result[i][j] = matrix[i][j] * scalar;
+ return result;
+ }
+
+ //! vector space multiplication with scalar
+ template <class Scalar, std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
+ friend auto operator*(Scalar scalar, const ThisType& matrix)
+ {
+ FieldMatrix<typename PromotionTraits<K, Scalar>::PromotedType, ROWS, COLS> result;
+ for (size_type i = 0; i < ROWS; ++i)
+ for (size_type j = 0; j < COLS; ++j)
+ result[i][j] = scalar * matrix[i][j];
+ return result;
+ }
+
+ //! vector space division by scalar
+ template <class Scalar, std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
+ friend auto operator/(const ThisType& matrix, Scalar scalar)
+ {
+ FieldMatrix<typename PromotionTraits<K, Scalar>::PromotedType, ROWS, COLS> result;
+ for (size_type i = 0; i < ROWS; ++i)
+ for (size_type j = 0; j < COLS; ++j)
+ result[i][j] = matrix[i][j] / scalar;
+ return result;
+ }
+
+ /** \brief Matrix-matrix multiplication
+ */
+ template <class OtherScalar, int otherCols>
+ friend auto operator*(const ThisType& matrixA, const Dune::FieldMatrix<OtherScalar, COLS, otherCols>& matrixB)
+ {
+ FieldMatrix<typename PromotionTraits<K, OtherScalar>::PromotedType, ROWS, otherCols> result;
+
+ for (size_type i = 0; i < matrixA.mat_rows(); ++i)
+ for (size_type j = 0; j < matrixB.mat_cols(); ++j) {
+ result[i][j] = 0;
+ for (size_type k = 0; k < matrixA.mat_cols(); ++k)
+ result[i][j] += matrixA[i][k] * matrixB[k][j];
+ }
+
+ return result;
+ }
+
+ template <class OtherScalar, int otherRows>
+ friend auto operator*(const Dune::FieldMatrix<OtherScalar, otherRows, ROWS>& matrixA, const ThisType& matrixB)
+ {
+ FieldMatrix<typename PromotionTraits<K, OtherScalar>::PromotedType, otherRows, COLS> result;
+ for (size_type i = 0; i < matrixA.mat_rows(); ++i)
+ for (size_type j = 0; j < matrixB.mat_cols(); ++j) {
+ result[i][j] = 0;
+ for (size_type k = 0; k < matrixA.mat_cols(); ++k)
+ result[i][j] += matrixA[i][k] * matrixB[k][j];
+ }
+ return result;
+ }
+
+ template <class OtherScalar, int otherCols>
+ friend auto operator*(const ThisType& matrixA, const FieldMatrix<OtherScalar, COLS, otherCols>& matrixB)
+ {
+ FieldMatrix<typename PromotionTraits<K, OtherScalar>::PromotedType, ROWS, otherCols> result;
+ for (size_type i = 0; i < matrixA.mat_rows(); ++i)
+ for (size_type j = 0; j < matrixB.mat_cols(); ++j) {
+ result[i][j] = 0;
+ for (size_type k = 0; k < matrixA.mat_cols(); ++k)
+ result[i][j] += matrixA[i][k] * matrixB[k][j];
+ }
+ return result;
+ }
+}; // class FieldMatrix<...>
+
+
+/**
+ * \todo We need to implement all operators from the base which return the base, to rather return ourselfes!
+ */
+template <class K>
+class FieldMatrix<K, 1, 1> : public Dune::FieldMatrix<K, 1, 1>
+{
+ static const int ROWS = 1;
+ static const int COLS = 1;
+ using BaseType = Dune::FieldMatrix<K, ROWS, COLS>;
+ using ThisType = FieldMatrix;
+
+public:
+ using typename BaseType::size_type;
+
+ FieldMatrix(const K& kk = suitable_default<K>::value())
+ : BaseType()
+ {
+ (*this)[0][0] = kk;
+ }
+
+ FieldMatrix(std::initializer_list<std::initializer_list<K>> list_of_rows)
+ : BaseType()
+ {
+#ifndef NDEBUG
+ if (list_of_rows.size() != ROWS)
+ DUNE_THROW(Exceptions::wrong_input_given,
+ "You are trying to construct a FieldMatrix< ..., 1, 1 > (of "
+ << "static size) from a list modeling a matrix with " << list_of_rows.size() << " rows!");
+#endif // NDEBUG
+ for (auto row : list_of_rows) {
+#ifndef NDEBUG
+ if (row.size() != COLS)
+ DUNE_THROW(Exceptions::wrong_input_given,
+ "You are trying to construct a FieldMatrix< ..., 1, 1 > (of "
+ << "static size) from a list with a row of length " << row.size() << "!");
+#endif // NDEBUG
+ for (auto entry : row)
+ (*this)[0][0] = entry;
+ }
+ } // FieldMatrix(std::initializer_list<std::initializer_list<K>> list_of_rows)
+
+ FieldMatrix(const size_t DXTC_DEBUG_ONLY(rr),
+ const size_t DXTC_DEBUG_ONLY(cc),
+ const K& kk = suitable_default<K>::value())
+ : BaseType()
+ {
+#ifndef NDEBUG
+ if (rr != ROWS || cc != COLS)
+ DUNE_THROW(Exceptions::wrong_input_given,
+ "You are trying to construct a FieldMatrix< ..., " << ROWS << ", " << COLS << " > (of "
+ << "static size) with " << rr << " rows and " << cc
+ << " columns!");
+#endif // NDEBUG
+ (*this)[0][0] = kk;
+ } // ... FieldMatrix(...)
+
+ FieldMatrix(const BaseType& other)
+ : BaseType(other)
+ {}
+
+ FieldMatrix(const Dune::XT::Common::FieldVector<K, 1>& other)
+ : BaseType()
+ {
+ (*this)[0][0] = other[0];
+ }
+
+ FieldMatrix(const Dune::FieldVector<K, 1>& other)
+ : BaseType()
+ {
+ (*this)[0][0] = other[0];
+ }
+
+ Dune::XT::Common::FieldMatrix<K, COLS, ROWS> transpose() const
+ {
+ return *this;
+ }
+
+ using BaseType::operator=;
+
+ ThisType& operator=(const FieldVector<K, 1>& other)
+ {
+ BaseType::operator=(other[0]);
+ return *this;
+ }
+
+ //! vector space addition -- two-argument version
+ template <class OtherScalar>
+ friend auto operator+(const ThisType& matrixA, const Dune::FieldMatrix<OtherScalar, 1, 1>& matrixB)
+ {
+ return FieldMatrix<typename PromotionTraits<K, OtherScalar>::PromotedType, 1, 1>{matrixA[0][0] + matrixB[0][0]};
+ }
+
+ //! Binary addition when treating FieldMatrix<K,1,1> like K
+ template <class Scalar, std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
+ friend auto operator+(const ThisType& matrix, const Scalar& scalar)
+ {
+ return FieldMatrix<typename PromotionTraits<K, Scalar>::PromotedType, 1, 1>{matrix[0][0] + scalar};
+ }
+
+ //! Binary addition when treating FieldMatrix<K,1,1> like K
+ template <class Scalar, std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
+ friend auto operator+(const Scalar& scalar, const ThisType& matrix)
+ {
+ return FieldMatrix<typename PromotionTraits<Scalar, K>::PromotedType, 1, 1>{scalar + matrix[0][0]};
+ }
+
+ //! vector space subtraction -- two-argument version
+ template <class OtherScalar>
+ friend auto operator-(const ThisType& matrixA, const Dune::FieldMatrix<OtherScalar, 1, 1>& matrixB)
+ {
+ return FieldMatrix<typename PromotionTraits<K, OtherScalar>::PromotedType, 1, 1>{matrixA[0][0] - matrixB[0][0]};
+ }
+
+ //! Binary subtraction when treating FieldMatrix<K,1,1> like K
+ template <class Scalar, std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
+ friend auto operator-(const ThisType& matrix, const Scalar& scalar)
+ {
+ return FieldMatrix<typename PromotionTraits<K, Scalar>::PromotedType, 1, 1>{matrix[0][0] - scalar};
+ }
+
+ //! Binary subtraction when treating FieldMatrix<K,1,1> like K
+ template <class Scalar, std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
+ friend auto operator-(const Scalar& scalar, const ThisType& matrix)
+ {
+ return FieldMatrix<typename PromotionTraits<Scalar, K>::PromotedType, 1, 1>{scalar - matrix[0][0]};
+ }
+
+ //! vector space multiplication with scalar
+ template <class Scalar, std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
+ friend auto operator*(const ThisType& matrix, Scalar scalar)
+ {
+ return FieldMatrix<typename PromotionTraits<K, Scalar>::PromotedType, 1, 1>{matrix[0][0] * scalar};
+ }
+
+ //! vector space multiplication with scalar
+ template <class Scalar, std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
+ friend auto operator*(Scalar scalar, const ThisType& matrix)
+ {
+ return FieldMatrix<typename PromotionTraits<K, Scalar>::PromotedType, 1, 1>{scalar * matrix[0][0]};
+ }
+
+ //! vector space division by scalar
+ template <class Scalar, std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
+ friend auto operator/(const ThisType& matrix, Scalar scalar)
+ {
+ return FieldMatrix<typename PromotionTraits<K, Scalar>::PromotedType, 1, 1>{matrix[0][0] / scalar};
+ }
+
+ /** \brief Matrix-matrix multiplication
+ */
+ template <class OtherScalar, int otherCols>
+ friend auto operator*(const ThisType& matrixA, const Dune::FieldMatrix<OtherScalar, 1, otherCols>& matrixB)
+ {
+ FieldMatrix<typename PromotionTraits<K, OtherScalar>::PromotedType, 1, otherCols> result;
+ for (size_type j = 0; j < matrixB.mat_cols(); ++j)
+ result[0][j] = matrixA[0][0] * matrixB[0][j];
+ return result;
+ }
+
+ template <class OtherScalar, int otherRows>
+ friend auto operator*(const Dune::FieldMatrix<OtherScalar, otherRows, 1>& matrixA, const ThisType& matrixB)
+ {
+ FieldMatrix<typename PromotionTraits<K, OtherScalar>::PromotedType, otherRows, 1> result;
+ for (size_type i = 0; i < matrixA.mat_rows(); ++i)
+ result[i][0] = matrixA[i][0] * matrixB[0][0];
+ return result;
+ }
+
+ template <class OtherScalar, int otherCols>
+ friend auto operator*(const ThisType& matrixA, const FieldMatrix<OtherScalar, 1, otherCols>& matrixB)
+ {
+ FieldMatrix<typename PromotionTraits<K, OtherScalar>::PromotedType, 1, otherCols> result;
+ for (size_type j = 0; j < matrixB.mat_cols(); ++j)
+ result[0][j] = matrixA[0][0] * matrixB[0][j];
+ return result;
+ }
+}; // class FieldMatrix
+
+
+template <class K, size_t num_blocks, size_t block_rows, size_t block_cols = block_rows>
+class BlockedFieldMatrix
+{
+ using ThisType = BlockedFieldMatrix;
+
+public:
+ static constexpr size_t num_rows = num_blocks * block_rows;
+ static constexpr size_t num_cols = num_blocks * block_cols;
+ using MatrixType = Dune::FieldMatrix<K, num_rows, num_cols>;
+ using BlockType = FieldMatrix<K, block_rows, block_cols>;
+
+ BlockedFieldMatrix(const K& val = K(0.))
+ : backend_(BlockType(val))
+ {}
+
+ BlockedFieldMatrix(const size_t DXTC_DEBUG_ONLY(rows), const size_t DXTC_DEBUG_ONLY(cols), const K& val = K(0.))
+ : backend_(BlockType(val))
+ {
+ assert(rows == num_rows && cols == num_cols && "Requested shape has to match static shape!");
+ }
+
+ template <class OtherMatrixType>
+ BlockedFieldMatrix(const OtherMatrixType& other,
+ std::enable_if_t<is_matrix<OtherMatrixType>::value, int> /*dummy*/ = 0)
+ {
+ *this = other;
+ }
+
+ template <class OtherMatrixType>
+ std::enable_if_t<is_matrix<OtherMatrixType>::value, ThisType>& operator=(const OtherMatrixType& other)
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ for (size_t ll = 0; ll < block_rows; ++ll)
+ for (size_t mm = 0; mm < block_cols; ++mm)
+ backend_[jj][ll][mm] =
+ MatrixAbstraction<OtherMatrixType>::get_entry(other, jj * block_rows + ll, jj * block_cols + mm);
+ return *this;
+ }
+
+
+ BlockedFieldMatrix(const MatrixType& other)
+ {
+ *this = other;
+ }
+
+ BlockedFieldMatrix(const BlockType& block)
+ : backend_(block)
+ {}
+
+ bool operator==(const ThisType& other) const
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ if (block(jj) != other.block(jj))
+ return false;
+ return true;
+ }
+
+ K get_entry(const size_t ii, const size_t jj) const
+ {
+ assert(ii < num_rows && jj < num_cols);
+ if (!in_pattern(ii, jj))
+ return K(0.);
+ return backend_[ii / block_rows][ii % block_rows][jj % block_cols];
+ }
+
+ K& get_entry(const size_t jj, const size_t ll, const size_t mm)
+ {
+ assert(is_valid_entry(jj, ll, mm));
+ return backend_[jj][ll][mm];
+ }
+
+ const K& get_entry(const size_t jj, const size_t ll, const size_t mm) const
+ {
+ assert(is_valid_entry(jj, ll, mm));
+ return backend_[jj][ll][mm];
+ }
+
+ void set_entry(const size_t ii, const size_t jj, const K& val)
+ {
+ assert(ii < num_rows && jj < num_cols);
+ if (in_pattern(ii, jj))
+ backend_[ii / block_rows][ii % block_rows][jj % block_cols] = val;
+ else if (XT::Common::FloatCmp::ne(val, K(0)))
+ DUNE_THROW(Dune::MathError, "Tried to modify a value that is not in the pattern!");
+ }
+
+ void set_entry(const size_t jj, const size_t ll, const size_t mm, const K& val)
+ {
+ assert(is_valid_entry(jj, ll, mm));
+ backend_[jj][ll][mm] = val;
+ }
+
+ void add_to_entry(const size_t ii, const size_t jj, const K& val)
+ {
+ assert(ii < num_rows && jj < num_cols);
+ if (!in_pattern(ii, jj) && XT::Common::FloatCmp::ne(val, K(0)))
+ DUNE_THROW(Dune::MathError, "Tried to modify a value that is not in the pattern!");
+ backend_[ii / block_rows][ii % block_rows][jj % block_cols] += val;
+ }
+
+ void add_to_entry(const size_t jj, const size_t ll, const size_t mm, const K& val)
+ {
+ assert(is_valid_entry(jj, ll, mm));
+ backend_[jj][ll][mm] += val;
+ }
+
+ BlockType& block(const size_t jj)
+ {
+ assert(jj < num_blocks);
+ return backend_[jj];
+ }
+
+ const BlockType& block(const size_t jj) const
+ {
+ assert(jj < num_blocks);
+ return backend_[jj];
+ }
+
+ void mv(const Dune::FieldVector<K, num_cols>& x, Dune::FieldVector<K, num_rows>& ret) const
+ {
+ std::fill(ret.begin(), ret.end(), 0.);
+ for (size_t jj = 0; jj < num_blocks; ++jj) {
+ const auto row_offset = block_rows * jj;
+ const auto col_offset = block_cols * jj;
+ for (size_t ll = 0; ll < block_rows; ++ll)
+ for (size_t mm = 0; mm < block_cols; ++mm)
+ ret[row_offset + ll] += backend_[jj][ll][mm] * x[col_offset + mm];
+ } // jj
+ } // void mv(...)
+
+ void mv(const BlockedFieldVector<K, num_blocks, block_cols>& x,
+ BlockedFieldVector<K, num_blocks, block_rows>& ret) const
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ backend_[jj].mv(x.block(jj), ret.block(jj));
+ } // void mv(...)
+
+ void mtv(const Dune::FieldVector<K, num_rows>& x, Dune::FieldVector<K, num_cols>& ret) const
+ {
+ std::fill(ret.begin(), ret.end(), 0.);
+ for (size_t jj = 0; jj < num_blocks; ++jj) {
+ const auto row_offset = block_rows * jj;
+ const auto col_offset = block_cols * jj;
+ for (size_t mm = 0; mm < block_cols; ++mm)
+ for (size_t ll = 0; ll < block_rows; ++ll)
+ ret[col_offset + mm] += backend_[jj][ll][mm] * x[row_offset + ll];
+ } // jj
+ } // void mtv(...)
+
+ void mtv(const BlockedFieldVector<K, num_blocks, block_rows>& x,
+ BlockedFieldVector<K, num_blocks, block_cols>& ret) const
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ backend_[jj].mtv(x.block(jj), ret.block(jj));
+ } // void mv(...)
+
+ template <size_t br, size_t bc>
+ ThisType& rightmultiply(const BlockedFieldMatrix<K, num_blocks, br, bc>& other)
+ {
+ assert((this != &other) && "Multiplying a matrix by itself gives wrong results, please copy before!");
+ static_assert(br == bc, "Cannot rightmultiply with non-square matrix");
+ static_assert(br == block_cols, "Size mismatch");
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ backend_[jj].rightmultiply(other.backend_[jj]);
+ return *this;
+ }
+
+ BlockedFieldMatrix<K, num_blocks, block_cols, block_rows> transpose()
+ {
+ BlockedFieldMatrix<K, num_blocks, block_cols, block_rows> ret;
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ ret.block(jj) = block(jj).transpose();
+ return ret;
+ }
+
+ ThisType operator*(const ThisType& other) const
+ {
+ ThisType ret(*this);
+ ret.rightmultiply(other);
+ return ret;
+ }
+
+ ThisType& operator*=(const K& val)
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ block(jj) *= val;
+ return *this;
+ }
+
+ ThisType& operator+=(const ThisType& other)
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ block(jj) += other.block(jj);
+ return *this;
+ }
+
+ ThisType operator+(const ThisType& other) const
+ {
+ ThisType ret(*this);
+ ret += other;
+ return ret;
+ }
+
+ ThisType& operator-=(const ThisType& other)
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ block(jj) -= other.block(jj);
+ return *this;
+ }
+
+ ThisType operator-(const ThisType& other) const
+ {
+ ThisType ret(*this);
+ ret -= other;
+ return ret;
+ }
+
+ static bool in_pattern(const size_t ii, const size_t jj)
+ {
+ return (ii / block_rows == jj / block_cols);
+ }
+
+ static bool is_valid_entry(const size_t jj, const size_t ll, const size_t mm)
+ {
+ return (jj < num_blocks && ll < block_rows && mm < block_cols);
+ }
+
+ DynamicMatrix<K> convert_to_dynamic_matrix() const
+ {
+ DynamicMatrix<K> ret(num_rows, num_cols, 0.);
+ for (size_t jj = 0; jj < num_blocks; ++jj) {
+ const size_t row_offset = jj * block_rows;
+ const size_t col_offset = jj * block_cols;
+ for (size_t rr = 0; rr < block_rows; ++rr)
+ for (size_t cc = 0; cc < block_cols; ++cc)
+ ret[row_offset + rr][col_offset + cc] = block(jj)[rr][cc];
+ } // jj
+ return ret;
+ }
+
+ template <class CharType, class CharTraits>
+ friend std::basic_ostream<CharType, CharTraits>& operator<<(std::basic_ostream<CharType, CharTraits>& out,
+ const ThisType& mat)
+ {
+ return output_matrix(out, mat);
+ } // ... operator<<(...)
+
+private:
+ FieldVector<BlockType, num_blocks> backend_;
+};
+
+
+template <class K, int N, int M>
+struct MatrixAbstraction<Dune::XT::Common::FieldMatrix<K, N, M>>
+{
+ typedef Dune::XT::Common::FieldMatrix<K, N, M> MatrixType;
+ typedef typename Dune::FieldTraits<K>::field_type ScalarType;
+ typedef typename Dune::FieldTraits<K>::real_type RealType;
+ typedef ScalarType S;
+ typedef RealType R;
+ template <size_t rows = N, size_t cols = M, class FieldType = K>
+ using MatrixTypeTemplate = Dune::XT::Common::FieldMatrix<FieldType, rows, cols>;
+
+ static const bool is_matrix = true;
+
+ static const bool has_static_size = true;
+
+ static const size_t static_rows = N;
+
+ static const size_t static_cols = M;
+
+ static const constexpr StorageLayout storage_layout = StorageLayout::dense_row_major;
+
+ static constexpr bool has_ostream = true;
+
+ template <class SparsityPatternType = FullPattern>
+ static inline MatrixType create(const size_t rows,
+ const size_t cols,
+ const ScalarType& val = suitable_default<ScalarType>::value(),
+ const SparsityPatternType& /*pattern*/ = SparsityPatternType())
+ {
+ return MatrixType(rows, cols, val);
+ }
+
+ template <class SparsityPatternType = FullPattern>
+ static inline std::unique_ptr<MatrixType> make_unique(const size_t rows,
+ const size_t cols,
+ const ScalarType& val = suitable_default<ScalarType>::value(),
+ const SparsityPatternType& /*pattern*/ = SparsityPatternType())
+ {
+ return std::make_unique<MatrixType>(rows, cols, val);
+ }
+
+ static constexpr size_t rows(const MatrixType& /*mat*/)
+ {
+ return N;
+ }
+
+ static constexpr size_t cols(const MatrixType& /*mat*/)
+ {
+ return M;
+ }
+
+ static inline void set_entry(MatrixType& mat, const size_t row, const size_t col, const ScalarType& val)
+ {
+ mat[row][col] = val;
+ }
+
+ static inline ScalarType get_entry(const MatrixType& mat, const size_t row, const size_t col)
+ {
+ return mat[row][col];
+ }
+
+ static inline void add_to_entry(MatrixType& mat, const size_t row, const size_t col, const ScalarType& val)
+ {
+ mat[row][col] += val;
+ }
+
+
+ static inline ScalarType* data(MatrixType& mat)
+ {
+ return &(mat[0][0]);
+ }
+
+ static inline const ScalarType* data(const MatrixType& mat)
+ {
+ return &(mat[0][0]);
+ }
+};
+
+template <class K, size_t num_blocks, size_t block_rows, size_t block_cols>
+struct MatrixAbstraction<Dune::XT::Common::BlockedFieldMatrix<K, num_blocks, block_rows, block_cols>>
+{
+ typedef Dune::XT::Common::BlockedFieldMatrix<K, num_blocks, block_rows, block_cols> MatrixType;
+ typedef typename Dune::FieldTraits<K>::field_type ScalarType;
+ typedef typename Dune::FieldTraits<K>::real_type RealType;
+ typedef ScalarType S;
+ typedef RealType R;
+
+ static const bool is_matrix = true;
+
+ static const bool has_static_size = true;
+
+ static const size_t static_rows = MatrixType::num_rows;
+
+ static const size_t static_cols = MatrixType::num_cols;
+
+ template <size_t rows = static_rows, size_t cols = static_cols, class FieldType = K>
+ using MatrixTypeTemplate = Dune::XT::Common::BlockedFieldMatrix<FieldType, rows / block_rows, block_rows, block_cols>;
+
+ static const constexpr StorageLayout storage_layout = StorageLayout::other;
+
+ static constexpr bool has_ostream = true;
+
+ template <class SparsityPatternType = FullPattern>
+ static inline MatrixType create(const size_t rows,
+ const size_t cols,
+ const ScalarType& val = suitable_default<ScalarType>::value(),
+ const SparsityPatternType& /*pattern*/ = SparsityPatternType())
+ {
+ return MatrixType(rows, cols, val);
+ }
+
+ template <class SparsityPatternType = FullPattern>
+ static inline std::unique_ptr<MatrixType> make_unique(const size_t rows,
+ const size_t cols,
+ const ScalarType& val = suitable_default<ScalarType>::value(),
+ const SparsityPatternType& /*pattern*/ = SparsityPatternType())
+ {
+ return std::make_unique<MatrixType>(rows, cols, val);
+ }
+
+ static constexpr size_t rows(const MatrixType& /*mat*/)
+ {
+ return static_rows;
+ }
+
+ static constexpr size_t cols(const MatrixType& /*mat*/)
+ {
+ return static_cols;
+ }
+
+ static inline void set_entry(MatrixType& mat, const size_t row, const size_t col, const ScalarType& val)
+ {
+ mat.set_entry(row, col, val);
+ }
+
+ static inline ScalarType get_entry(const MatrixType& mat, const size_t row, const size_t col)
+ {
+ return mat.get_entry(row, col);
+ }
+
+ static inline void add_to_entry(MatrixType& mat, const size_t row, const size_t col, const ScalarType& val)
+ {
+ mat.add_to_entry(row, col, val);
+ }
+
+ static inline ScalarType* data(MatrixType& /*mat*/)
+ {
+ DUNE_THROW(InvalidStateException, "Do not call me if storage layout is not dense!");
+ return nullptr;
+ }
+
+ static inline const ScalarType* data(const MatrixType& /*mat*/)
+ {
+ DUNE_THROW(InvalidStateException, "Do not call me if storage_layout is not dense!");
+ return nullptr;
+ }
+};
+
+
+template <class M>
+typename std::enable_if<is_matrix<M>::value && MatrixAbstraction<M>::has_static_size,
+ std::unique_ptr<FieldMatrix<typename MatrixAbstraction<M>::S,
+ MatrixAbstraction<M>::static_rows,
+ MatrixAbstraction<M>::static_cols>>>::type
+make_field_container_ptr(const M& mat)
+{
+ static const size_t rows = MatrixAbstraction<M>::static_rows;
+ static const size_t cols = MatrixAbstraction<M>::static_cols;
+ auto ret = std::make_unique<FieldMatrix<typename MatrixAbstraction<M>::S, rows, cols>>;
+ for (size_t ii = 0; ii < rows; ++ii)
+ for (size_t jj = 0; jj < cols; ++jj)
+ (*ret)[ii][jj] = get_matrix_entry(mat, ii, jj);
+ return std::move(ret);
+}
+
+
+template <class M>
+typename std::enable_if<is_matrix<M>::value && MatrixAbstraction<M>::has_static_size,
+ FieldMatrix<typename MatrixAbstraction<M>::S,
+ MatrixAbstraction<M>::static_rows,
+ MatrixAbstraction<M>::static_cols>>::type
+make_field_container(const M& mat)
+{
+ static const size_t rows = MatrixAbstraction<M>::static_rows;
+ static const size_t cols = MatrixAbstraction<M>::static_cols;
+ FieldMatrix<typename MatrixAbstraction<M>::S, rows, cols> ret;
+ for (size_t ii = 0; ii < rows; ++ii)
+ for (size_t jj = 0; jj < cols; ++jj)
+ ret[ii][jj] = get_matrix_entry(mat, ii, jj);
+ return ret;
+}
+
+
+template <class K, int ROWS, int COLS>
+FieldMatrix<K, ROWS, COLS> make_field_container(Dune::FieldMatrix<K, ROWS, COLS>&& vec)
+{
+ return std::move(vec);
+}
+
+
+template <class K, int ROWS, int COLS>
+FieldMatrix<K, ROWS, COLS> make_field_container(FieldMatrix<K, ROWS, COLS>&& vec)
+{
+ return std::move(vec);
+}
+
+
+template <class K, int ROWS, int COLS>
+typename std::enable_if<is_arithmetic<K>::value && !is_complex<K>::value, FieldMatrix<K, ROWS, COLS>>::type
+real(const Dune::FieldMatrix<K, ROWS, COLS>& real_mat)
+{
+ return real_mat;
+}
+
+template <class K, int ROWS, int COLS>
+typename std::enable_if<is_arithmetic<K>::value && !is_complex<K>::value, FieldMatrix<K, ROWS, COLS>>::type
+real(const FieldMatrix<K, ROWS, COLS>& real_mat)
+{
+ return real_mat;
+}
+
+template <class K, int ROWS, int COLS>
+typename std::enable_if<is_arithmetic<K>::value && !is_complex<K>::value, FieldMatrix<K, ROWS, COLS>>::type
+real(Dune::FieldMatrix<K, ROWS, COLS>&& real_mat)
+{
+ return std::move(real_mat);
+}
+
+template <class K, int ROWS, int COLS>
+typename std::enable_if<is_arithmetic<K>::value && !is_complex<K>::value, FieldMatrix<K, ROWS, COLS>>::type
+real(FieldMatrix<K, ROWS, COLS>&& real_mat)
+{
+ return std::move(real_mat);
+}
+
+template <class K, int ROWS, int COLS>
+typename std::enable_if<is_complex<K>::value, FieldMatrix<real_t<K>, ROWS, COLS>>::type
+real(const Dune::FieldMatrix<K, ROWS, COLS>& complex_mat)
+{
+ FieldMatrix<real_t<K>, ROWS, COLS> real_mat;
+ for (size_t ii = 0; ii < ROWS; ++ii)
+ for (size_t jj = 0; jj < COLS; ++jj)
+ real_mat[ii][jj] = complex_mat[ii][jj].real();
+ return real_mat;
+}
+
+template <class K, int ROWS, int COLS>
+typename std::enable_if<is_complex<K>::value, FieldMatrix<real_t<K>, ROWS, COLS>>::type
+real(const FieldMatrix<K, ROWS, COLS>& complex_mat)
+{
+ FieldMatrix<real_t<K>, ROWS, COLS> real_mat;
+ for (size_t ii = 0; ii < ROWS; ++ii)
+ for (size_t jj = 0; jj < COLS; ++jj)
+ real_mat[ii][jj] = complex_mat[ii][jj].real();
+ return real_mat;
+}
+
+
+template <class K, int ROWS, int COLS>
+typename std::enable_if<is_arithmetic<K>::value && !is_complex<K>::value, FieldMatrix<K, ROWS, COLS>>::type
+imag(const Dune::FieldMatrix<K, ROWS, COLS>& /*real_mat*/)
+{
+ return FieldMatrix<K, ROWS, COLS>(0);
+}
+
+template <class K, int ROWS, int COLS>
+typename std::enable_if<is_arithmetic<K>::value && !is_complex<K>::value, FieldMatrix<K, ROWS, COLS>>::type
+imag(const FieldMatrix<K, ROWS, COLS>& /*real_mat*/)
+{
+ return FieldMatrix<K, ROWS, COLS>(0);
+}
+
+template <class K, int ROWS, int COLS>
+typename std::enable_if<is_complex<K>::value, FieldMatrix<real_t<K>, ROWS, COLS>>::type
+imag(const Dune::FieldMatrix<K, ROWS, COLS>& complex_mat)
+{
+ FieldMatrix<real_t<K>, ROWS, COLS> real_mat;
+ for (size_t ii = 0; ii < ROWS; ++ii)
+ for (size_t jj = 0; jj < COLS; ++jj)
+ real_mat[ii][jj] = complex_mat[ii][jj].imag();
+ return real_mat;
+}
+
+template <class K, int ROWS, int COLS>
+typename std::enable_if<is_complex<K>::value, FieldMatrix<real_t<K>, ROWS, COLS>>::type
+imag(const FieldMatrix<K, ROWS, COLS>& complex_mat)
+{
+ FieldMatrix<real_t<K>, ROWS, COLS> real_mat;
+ for (size_t ii = 0; ii < ROWS; ++ii)
+ for (size_t jj = 0; jj < COLS; ++jj)
+ real_mat[ii][jj] = complex_mat[ii][jj].imag();
+ return real_mat;
+}
+
+
+} // namespace Common
+} // namespace XT
+
+template <class K, int L_ROWS, int L_COLS, int R_COLS>
+void rightmultiply(Dune::FieldMatrix<K, L_ROWS, R_COLS>& ret,
+ const Dune::FieldMatrix<K, L_ROWS, L_COLS>& left,
+ const Dune::FieldMatrix<K, L_COLS, R_COLS>& right)
+{
+ for (size_t ii = 0; ii < L_ROWS; ++ii) {
+ for (size_t jj = 0; jj < R_COLS; ++jj) {
+ ret[ii][jj] = 0.;
+ for (size_t kk = 0; kk < L_COLS; ++kk)
+ ret[ii][jj] += left[ii][kk] * right[kk][jj];
+ }
+ }
+}
+
+// versions that do not allocate matrices on the stack (for large matrices)
+template <class K, int L_ROWS, int L_COLS, int R_COLS>
+std::unique_ptr<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>
+operator*(const std::unique_ptr<Dune::FieldMatrix<K, L_ROWS, L_COLS>>& left,
+ const Dune::FieldMatrix<K, L_COLS, R_COLS>& right)
+{
+ auto ret = std::make_unique<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>();
+ rightmultiply(*ret, *left, right);
+ return ret;
+}
+
+template <class K, int L_ROWS, int L_COLS, int R_COLS>
+std::unique_ptr<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>
+operator*(const Dune::FieldMatrix<K, L_ROWS, L_COLS>& left,
+ const std::unique_ptr<Dune::FieldMatrix<K, L_COLS, R_COLS>>& right)
+{
+ auto ret = std::make_unique<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>();
+ rightmultiply(*ret, left, *right);
+ return ret;
+}
+
+template <class K, int L_ROWS, int L_COLS, int R_COLS>
+std::unique_ptr<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>
+operator*(const std::unique_ptr<Dune::FieldMatrix<K, L_ROWS, L_COLS>>& left,
+ const std::unique_ptr<Dune::FieldMatrix<K, L_COLS, R_COLS>>& right)
+{
+ return left * *right;
+}
+
+
+} // namespace Dune
+
+#endif // DUNE_XT_COMMON_FMATRIX_27_HH
diff --git a/dune/xt/common/fmatrix.hh b/dune/xt/common/fmatrix.hh
index a87b6bf5172d3473e2307194472c9eed5c54652f..5b311f6a63395d2b10badcf1d99e7046f2a32239 100644
--- a/dune/xt/common/fmatrix.hh
+++ b/dune/xt/common/fmatrix.hh
@@ -13,996 +13,12 @@
#ifndef DUNE_XT_COMMON_FMATRIX_HH
#define DUNE_XT_COMMON_FMATRIX_HH
-#include <initializer_list>
+#include <dune/common/version.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fvector.hh>
-
-#include <dune/xt/common/exceptions.hh>
-#include <dune/xt/common/debug.hh>
-#include <dune/xt/common/matrix.hh>
-#include <dune/xt/common/fvector.hh>
-#include <dune/xt/common/type_traits.hh>
-
-namespace Dune {
-namespace XT {
-namespace Common {
-
-
-/**
- * \todo We need to implement all operators from the base which return the base, to rather return ourselfes!
- */
-template <class K, int ROWS, int COLS>
-class FieldMatrix : public Dune::FieldMatrix<K, ROWS, COLS>
-{
- using BaseType = Dune::FieldMatrix<K, ROWS, COLS>;
- using ThisType = FieldMatrix;
-
-public:
- using typename BaseType::field_type;
- using typename BaseType::size_type;
- using typename BaseType::value_type;
-
- FieldMatrix(const K& kk = suitable_default<K>::value())
- : BaseType()
- {
- // This is required because BaseType(kk) does not work for std::string
- for (size_t i = 0; i < ROWS; ++i) {
- for (size_t j = 0; j < COLS; ++j)
- (*this)[i][j] = kk;
- }
- }
-
- FieldMatrix(const size_t DXTC_DEBUG_ONLY(rr),
- const size_t DXTC_DEBUG_ONLY(cc),
- const K& kk = suitable_default<K>::value())
- : BaseType()
- {
-#ifndef NDEBUG
- if (rr != ROWS || cc != COLS)
- DUNE_THROW(Exceptions::wrong_input_given,
- "You are trying to construct a FieldMatrix< ..., " << ROWS << ", " << COLS << " > (of "
- << "static size) with " << rr << " rows and " << cc
- << " columns!");
-#endif // NDEBUG
- // This is required because BaseType(kk) does not work for std::string
- for (size_t i = 0; i < ROWS; ++i) {
- for (size_t j = 0; j < COLS; ++j)
- (*this)[i][j] = kk;
- }
-
- } // ... FieldMatrix(...)
-
- FieldMatrix(std::initializer_list<std::initializer_list<K>> list_of_rows)
- : BaseType()
- {
-#ifndef NDEBUG
- if (list_of_rows.size() != ROWS)
- DUNE_THROW(Exceptions::wrong_input_given,
- "You are trying to construct a FieldMatrix< ..., "
- << ROWS << ", " << COLS << " > (of "
- << "static size) from a list modeling a matrix with " << list_of_rows.size() << " rows!");
-#endif // NDEBUG
- size_t rr = 0;
- for (auto row : list_of_rows) {
-#ifndef NDEBUG
- if (row.size() != COLS)
- DUNE_THROW(Exceptions::wrong_input_given,
- "You are trying to construct a FieldMatrix< ..., "
- << ROWS << ", " << COLS << " > (of "
- << "static size) from a list with a row of length " << row.size() << "!");
-#endif // NDEBUG
- size_t cc = 0;
- for (auto entry : row)
- (*this)[rr][cc++] = entry;
- ++rr;
- }
- } // FieldMatrix(std::initializer_list<std::initializer_list<K>> list_of_rows)
-
- template <class OtherMatrixType>
- FieldMatrix(const OtherMatrixType& other, std::enable_if_t<is_matrix<OtherMatrixType>::value, int> /*dummy*/ = 0)
- {
- *this = other;
- }
-
- template <class OtherMatrixType>
- std::enable_if_t<is_matrix<OtherMatrixType>::value, ThisType>& operator=(const OtherMatrixType& other)
- {
- for (size_t rr = 0; rr < ROWS; ++rr)
- for (size_t cc = 0; cc < COLS; ++cc)
- (*this)[rr][cc] = MatrixAbstraction<OtherMatrixType>::get_entry(other, rr, cc);
- return *this;
- }
-
- Dune::XT::Common::FieldMatrix<K, COLS, ROWS> transpose() const
- {
- Dune::XT::Common::FieldMatrix<K, COLS, ROWS> ret;
- for (size_t rr = 0; rr < ROWS; ++rr)
- for (size_t cc = 0; cc < COLS; ++cc)
- ret[cc][rr] = (*this)[rr][cc];
- return ret;
- }
-
- Dune::XT::Common::FieldVector<K, ROWS> operator*(const Dune::FieldVector<K, COLS>& vec) const
- {
- Dune::FieldVector<K, ROWS> ret;
- this->mv(vec, ret);
- return ret;
- }
-
- //! vector space addition -- two-argument version
- template <class OtherScalar>
- friend auto operator+(const ThisType& matrixA, const Dune::FieldMatrix<OtherScalar, ROWS, COLS>& matrixB)
- {
- FieldMatrix<typename PromotionTraits<K, OtherScalar>::PromotedType, ROWS, COLS> result;
- for (size_type i = 0; i < ROWS; ++i)
- for (size_type j = 0; j < COLS; ++j)
- result[i][j] = matrixA[i][j] + matrixB[i][j];
- return result;
- }
-
- //! vector space subtraction -- two-argument version
- template <class OtherScalar>
- friend auto operator-(const ThisType& matrixA, const Dune::FieldMatrix<OtherScalar, ROWS, COLS>& matrixB)
- {
- FieldMatrix<typename PromotionTraits<K, OtherScalar>::PromotedType, ROWS, COLS> result;
- for (size_type i = 0; i < ROWS; ++i)
- for (size_type j = 0; j < COLS; ++j)
- result[i][j] = matrixA[i][j] - matrixB[i][j];
- return result;
- }
-
- //! vector space multiplication with scalar
- template <class Scalar, std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
- friend auto operator*(const ThisType& matrix, Scalar scalar)
- {
- FieldMatrix<typename PromotionTraits<K, Scalar>::PromotedType, ROWS, COLS> result;
- for (size_type i = 0; i < ROWS; ++i)
- for (size_type j = 0; j < COLS; ++j)
- result[i][j] = matrix[i][j] * scalar;
- return result;
- }
-
- //! vector space multiplication with scalar
- template <class Scalar, std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
- friend auto operator*(Scalar scalar, const ThisType& matrix)
- {
- FieldMatrix<typename PromotionTraits<K, Scalar>::PromotedType, ROWS, COLS> result;
- for (size_type i = 0; i < ROWS; ++i)
- for (size_type j = 0; j < COLS; ++j)
- result[i][j] = scalar * matrix[i][j];
- return result;
- }
-
- //! vector space division by scalar
- template <class Scalar, std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
- friend auto operator/(const ThisType& matrix, Scalar scalar)
- {
- FieldMatrix<typename PromotionTraits<K, Scalar>::PromotedType, ROWS, COLS> result;
- for (size_type i = 0; i < ROWS; ++i)
- for (size_type j = 0; j < COLS; ++j)
- result[i][j] = matrix[i][j] / scalar;
- return result;
- }
-
- /** \brief Matrix-matrix multiplication
- */
- template <class OtherScalar, int otherCols>
- friend auto operator*(const ThisType& matrixA, const Dune::FieldMatrix<OtherScalar, COLS, otherCols>& matrixB)
- {
- FieldMatrix<typename PromotionTraits<K, OtherScalar>::PromotedType, ROWS, otherCols> result;
-
- for (size_type i = 0; i < matrixA.mat_rows(); ++i)
- for (size_type j = 0; j < matrixB.mat_cols(); ++j) {
- result[i][j] = 0;
- for (size_type k = 0; k < matrixA.mat_cols(); ++k)
- result[i][j] += matrixA[i][k] * matrixB[k][j];
- }
-
- return result;
- }
-
- template <class OtherScalar, int otherRows>
- friend auto operator*(const Dune::FieldMatrix<OtherScalar, otherRows, ROWS>& matrixA, const ThisType& matrixB)
- {
- FieldMatrix<typename PromotionTraits<K, OtherScalar>::PromotedType, otherRows, COLS> result;
- for (size_type i = 0; i < matrixA.mat_rows(); ++i)
- for (size_type j = 0; j < matrixB.mat_cols(); ++j) {
- result[i][j] = 0;
- for (size_type k = 0; k < matrixA.mat_cols(); ++k)
- result[i][j] += matrixA[i][k] * matrixB[k][j];
- }
- return result;
- }
-
- template <class OtherScalar, int otherCols>
- friend auto operator*(const ThisType& matrixA, const FieldMatrix<OtherScalar, COLS, otherCols>& matrixB)
- {
- FieldMatrix<typename PromotionTraits<K, OtherScalar>::PromotedType, ROWS, otherCols> result;
- for (size_type i = 0; i < matrixA.mat_rows(); ++i)
- for (size_type j = 0; j < matrixB.mat_cols(); ++j) {
- result[i][j] = 0;
- for (size_type k = 0; k < matrixA.mat_cols(); ++k)
- result[i][j] += matrixA[i][k] * matrixB[k][j];
- }
- return result;
- }
-}; // class FieldMatrix<...>
-
-
-/**
- * \todo We need to implement all operators from the base which return the base, to rather return ourselfes!
- */
-template <class K>
-class FieldMatrix<K, 1, 1> : public Dune::FieldMatrix<K, 1, 1>
-{
- static const int ROWS = 1;
- static const int COLS = 1;
- using BaseType = Dune::FieldMatrix<K, ROWS, COLS>;
- using ThisType = FieldMatrix;
-
-public:
- using typename BaseType::size_type;
-
- FieldMatrix(const K& kk = suitable_default<K>::value())
- : BaseType()
- {
- (*this)[0][0] = kk;
- }
-
- FieldMatrix(std::initializer_list<std::initializer_list<K>> list_of_rows)
- : BaseType()
- {
-#ifndef NDEBUG
- if (list_of_rows.size() != ROWS)
- DUNE_THROW(Exceptions::wrong_input_given,
- "You are trying to construct a FieldMatrix< ..., 1, 1 > (of "
- << "static size) from a list modeling a matrix with " << list_of_rows.size() << " rows!");
-#endif // NDEBUG
- for (auto row : list_of_rows) {
-#ifndef NDEBUG
- if (row.size() != COLS)
- DUNE_THROW(Exceptions::wrong_input_given,
- "You are trying to construct a FieldMatrix< ..., 1, 1 > (of "
- << "static size) from a list with a row of length " << row.size() << "!");
-#endif // NDEBUG
- for (auto entry : row)
- (*this)[0][0] = entry;
- }
- } // FieldMatrix(std::initializer_list<std::initializer_list<K>> list_of_rows)
-
- FieldMatrix(const size_t DXTC_DEBUG_ONLY(rr),
- const size_t DXTC_DEBUG_ONLY(cc),
- const K& kk = suitable_default<K>::value())
- : BaseType()
- {
-#ifndef NDEBUG
- if (rr != ROWS || cc != COLS)
- DUNE_THROW(Exceptions::wrong_input_given,
- "You are trying to construct a FieldMatrix< ..., " << ROWS << ", " << COLS << " > (of "
- << "static size) with " << rr << " rows and " << cc
- << " columns!");
-#endif // NDEBUG
- (*this)[0][0] = kk;
- } // ... FieldMatrix(...)
-
- FieldMatrix(const BaseType& other)
- : BaseType(other)
- {}
-
- FieldMatrix(const Dune::XT::Common::FieldVector<K, 1>& other)
- : BaseType()
- {
- (*this)[0][0] = other[0];
- }
-
- FieldMatrix(const Dune::FieldVector<K, 1>& other)
- : BaseType()
- {
- (*this)[0][0] = other[0];
- }
-
- Dune::XT::Common::FieldMatrix<K, COLS, ROWS> transpose() const
- {
- return *this;
- }
-
- using BaseType::operator=;
-
- ThisType& operator=(const FieldVector<K, 1>& other)
- {
- BaseType::operator=(other[0]);
- return *this;
- }
-
- //! vector space addition -- two-argument version
- template <class OtherScalar>
- friend auto operator+(const ThisType& matrixA, const Dune::FieldMatrix<OtherScalar, 1, 1>& matrixB)
- {
- return FieldMatrix<typename PromotionTraits<K, OtherScalar>::PromotedType, 1, 1>{matrixA[0][0] + matrixB[0][0]};
- }
-
- //! Binary addition when treating FieldMatrix<K,1,1> like K
- template <class Scalar, std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
- friend auto operator+(const ThisType& matrix, const Scalar& scalar)
- {
- return FieldMatrix<typename PromotionTraits<K, Scalar>::PromotedType, 1, 1>{matrix[0][0] + scalar};
- }
-
- //! Binary addition when treating FieldMatrix<K,1,1> like K
- template <class Scalar, std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
- friend auto operator+(const Scalar& scalar, const ThisType& matrix)
- {
- return FieldMatrix<typename PromotionTraits<Scalar, K>::PromotedType, 1, 1>{scalar + matrix[0][0]};
- }
-
- //! vector space subtraction -- two-argument version
- template <class OtherScalar>
- friend auto operator-(const ThisType& matrixA, const Dune::FieldMatrix<OtherScalar, 1, 1>& matrixB)
- {
- return FieldMatrix<typename PromotionTraits<K, OtherScalar>::PromotedType, 1, 1>{matrixA[0][0] - matrixB[0][0]};
- }
-
- //! Binary subtraction when treating FieldMatrix<K,1,1> like K
- template <class Scalar, std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
- friend auto operator-(const ThisType& matrix, const Scalar& scalar)
- {
- return FieldMatrix<typename PromotionTraits<K, Scalar>::PromotedType, 1, 1>{matrix[0][0] - scalar};
- }
-
- //! Binary subtraction when treating FieldMatrix<K,1,1> like K
- template <class Scalar, std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
- friend auto operator-(const Scalar& scalar, const ThisType& matrix)
- {
- return FieldMatrix<typename PromotionTraits<Scalar, K>::PromotedType, 1, 1>{scalar - matrix[0][0]};
- }
-
- //! vector space multiplication with scalar
- template <class Scalar, std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
- friend auto operator*(const ThisType& matrix, Scalar scalar)
- {
- return FieldMatrix<typename PromotionTraits<K, Scalar>::PromotedType, 1, 1>{matrix[0][0] * scalar};
- }
-
- //! vector space multiplication with scalar
- template <class Scalar, std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
- friend auto operator*(Scalar scalar, const ThisType& matrix)
- {
- return FieldMatrix<typename PromotionTraits<K, Scalar>::PromotedType, 1, 1>{scalar * matrix[0][0]};
- }
-
- //! vector space division by scalar
- template <class Scalar, std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
- friend auto operator/(const ThisType& matrix, Scalar scalar)
- {
- return FieldMatrix<typename PromotionTraits<K, Scalar>::PromotedType, 1, 1>{matrix[0][0] / scalar};
- }
-
- /** \brief Matrix-matrix multiplication
- */
- template <class OtherScalar, int otherCols>
- friend auto operator*(const ThisType& matrixA, const Dune::FieldMatrix<OtherScalar, 1, otherCols>& matrixB)
- {
- FieldMatrix<typename PromotionTraits<K, OtherScalar>::PromotedType, 1, otherCols> result;
- for (size_type j = 0; j < matrixB.mat_cols(); ++j)
- result[0][j] = matrixA[0][0] * matrixB[0][j];
- return result;
- }
-
- template <class OtherScalar, int otherRows>
- friend auto operator*(const Dune::FieldMatrix<OtherScalar, otherRows, 1>& matrixA, const ThisType& matrixB)
- {
- FieldMatrix<typename PromotionTraits<K, OtherScalar>::PromotedType, otherRows, 1> result;
- for (size_type i = 0; i < matrixA.mat_rows(); ++i)
- result[i][0] = matrixA[i][0] * matrixB[0][0];
- return result;
- }
-
- template <class OtherScalar, int otherCols>
- friend auto operator*(const ThisType& matrixA, const FieldMatrix<OtherScalar, 1, otherCols>& matrixB)
- {
- FieldMatrix<typename PromotionTraits<K, OtherScalar>::PromotedType, 1, otherCols> result;
- for (size_type j = 0; j < matrixB.mat_cols(); ++j)
- result[0][j] = matrixA[0][0] * matrixB[0][j];
- return result;
- }
-}; // class FieldMatrix
-
-
-template <class K, size_t num_blocks, size_t block_rows, size_t block_cols = block_rows>
-class BlockedFieldMatrix
-{
- using ThisType = BlockedFieldMatrix;
-
-public:
- static constexpr size_t num_rows = num_blocks * block_rows;
- static constexpr size_t num_cols = num_blocks * block_cols;
- using MatrixType = Dune::FieldMatrix<K, num_rows, num_cols>;
- using BlockType = FieldMatrix<K, block_rows, block_cols>;
-
- BlockedFieldMatrix(const K& val = K(0.))
- : backend_(BlockType(val))
- {}
-
- BlockedFieldMatrix(const size_t DXTC_DEBUG_ONLY(rows), const size_t DXTC_DEBUG_ONLY(cols), const K& val = K(0.))
- : backend_(BlockType(val))
- {
- assert(rows == num_rows && cols == num_cols && "Requested shape has to match static shape!");
- }
-
- template <class OtherMatrixType>
- BlockedFieldMatrix(const OtherMatrixType& other,
- std::enable_if_t<is_matrix<OtherMatrixType>::value, int> /*dummy*/ = 0)
- {
- *this = other;
- }
-
- template <class OtherMatrixType>
- std::enable_if_t<is_matrix<OtherMatrixType>::value, ThisType>& operator=(const OtherMatrixType& other)
- {
- for (size_t jj = 0; jj < num_blocks; ++jj)
- for (size_t ll = 0; ll < block_rows; ++ll)
- for (size_t mm = 0; mm < block_cols; ++mm)
- backend_[jj][ll][mm] =
- MatrixAbstraction<OtherMatrixType>::get_entry(other, jj * block_rows + ll, jj * block_cols + mm);
- return *this;
- }
-
-
- BlockedFieldMatrix(const MatrixType& other)
- {
- *this = other;
- }
-
- BlockedFieldMatrix(const BlockType& block)
- : backend_(block)
- {}
-
- bool operator==(const ThisType& other) const
- {
- for (size_t jj = 0; jj < num_blocks; ++jj)
- if (block(jj) != other.block(jj))
- return false;
- return true;
- }
-
- K get_entry(const size_t ii, const size_t jj) const
- {
- assert(ii < num_rows && jj < num_cols);
- if (!in_pattern(ii, jj))
- return K(0.);
- return backend_[ii / block_rows][ii % block_rows][jj % block_cols];
- }
-
- K& get_entry(const size_t jj, const size_t ll, const size_t mm)
- {
- assert(is_valid_entry(jj, ll, mm));
- return backend_[jj][ll][mm];
- }
-
- const K& get_entry(const size_t jj, const size_t ll, const size_t mm) const
- {
- assert(is_valid_entry(jj, ll, mm));
- return backend_[jj][ll][mm];
- }
-
- void set_entry(const size_t ii, const size_t jj, const K& val)
- {
- assert(ii < num_rows && jj < num_cols);
- if (in_pattern(ii, jj))
- backend_[ii / block_rows][ii % block_rows][jj % block_cols] = val;
- else if (XT::Common::FloatCmp::ne(val, K(0)))
- DUNE_THROW(Dune::MathError, "Tried to modify a value that is not in the pattern!");
- }
-
- void set_entry(const size_t jj, const size_t ll, const size_t mm, const K& val)
- {
- assert(is_valid_entry(jj, ll, mm));
- backend_[jj][ll][mm] = val;
- }
-
- void add_to_entry(const size_t ii, const size_t jj, const K& val)
- {
- assert(ii < num_rows && jj < num_cols);
- if (!in_pattern(ii, jj) && XT::Common::FloatCmp::ne(val, K(0)))
- DUNE_THROW(Dune::MathError, "Tried to modify a value that is not in the pattern!");
- backend_[ii / block_rows][ii % block_rows][jj % block_cols] += val;
- }
-
- void add_to_entry(const size_t jj, const size_t ll, const size_t mm, const K& val)
- {
- assert(is_valid_entry(jj, ll, mm));
- backend_[jj][ll][mm] += val;
- }
-
- BlockType& block(const size_t jj)
- {
- assert(jj < num_blocks);
- return backend_[jj];
- }
-
- const BlockType& block(const size_t jj) const
- {
- assert(jj < num_blocks);
- return backend_[jj];
- }
-
- void mv(const Dune::FieldVector<K, num_cols>& x, Dune::FieldVector<K, num_rows>& ret) const
- {
- std::fill(ret.begin(), ret.end(), 0.);
- for (size_t jj = 0; jj < num_blocks; ++jj) {
- const auto row_offset = block_rows * jj;
- const auto col_offset = block_cols * jj;
- for (size_t ll = 0; ll < block_rows; ++ll)
- for (size_t mm = 0; mm < block_cols; ++mm)
- ret[row_offset + ll] += backend_[jj][ll][mm] * x[col_offset + mm];
- } // jj
- } // void mv(...)
-
- void mv(const BlockedFieldVector<K, num_blocks, block_cols>& x,
- BlockedFieldVector<K, num_blocks, block_rows>& ret) const
- {
- for (size_t jj = 0; jj < num_blocks; ++jj)
- backend_[jj].mv(x.block(jj), ret.block(jj));
- } // void mv(...)
-
- void mtv(const Dune::FieldVector<K, num_rows>& x, Dune::FieldVector<K, num_cols>& ret) const
- {
- std::fill(ret.begin(), ret.end(), 0.);
- for (size_t jj = 0; jj < num_blocks; ++jj) {
- const auto row_offset = block_rows * jj;
- const auto col_offset = block_cols * jj;
- for (size_t mm = 0; mm < block_cols; ++mm)
- for (size_t ll = 0; ll < block_rows; ++ll)
- ret[col_offset + mm] += backend_[jj][ll][mm] * x[row_offset + ll];
- } // jj
- } // void mtv(...)
-
- void mtv(const BlockedFieldVector<K, num_blocks, block_rows>& x,
- BlockedFieldVector<K, num_blocks, block_cols>& ret) const
- {
- for (size_t jj = 0; jj < num_blocks; ++jj)
- backend_[jj].mtv(x.block(jj), ret.block(jj));
- } // void mv(...)
-
- template <size_t br, size_t bc>
- ThisType& rightmultiply(const BlockedFieldMatrix<K, num_blocks, br, bc>& other)
- {
- assert((this != &other) && "Multiplying a matrix by itself gives wrong results, please copy before!");
- static_assert(br == bc, "Cannot rightmultiply with non-square matrix");
- static_assert(br == block_cols, "Size mismatch");
- for (size_t jj = 0; jj < num_blocks; ++jj)
- backend_[jj].rightmultiply(other.backend_[jj]);
- return *this;
- }
-
- BlockedFieldMatrix<K, num_blocks, block_cols, block_rows> transpose()
- {
- BlockedFieldMatrix<K, num_blocks, block_cols, block_rows> ret;
- for (size_t jj = 0; jj < num_blocks; ++jj)
- ret.block(jj) = block(jj).transpose();
- return ret;
- }
-
- ThisType operator*(const ThisType& other) const
- {
- ThisType ret(*this);
- ret.rightmultiply(other);
- return ret;
- }
-
- ThisType& operator*=(const K& val)
- {
- for (size_t jj = 0; jj < num_blocks; ++jj)
- block(jj) *= val;
- return *this;
- }
-
- ThisType& operator+=(const ThisType& other)
- {
- for (size_t jj = 0; jj < num_blocks; ++jj)
- block(jj) += other.block(jj);
- return *this;
- }
-
- ThisType operator+(const ThisType& other) const
- {
- ThisType ret(*this);
- ret += other;
- return ret;
- }
-
- ThisType& operator-=(const ThisType& other)
- {
- for (size_t jj = 0; jj < num_blocks; ++jj)
- block(jj) -= other.block(jj);
- return *this;
- }
-
- ThisType operator-(const ThisType& other) const
- {
- ThisType ret(*this);
- ret -= other;
- return ret;
- }
-
- static bool in_pattern(const size_t ii, const size_t jj)
- {
- return (ii / block_rows == jj / block_cols);
- }
-
- static bool is_valid_entry(const size_t jj, const size_t ll, const size_t mm)
- {
- return (jj < num_blocks && ll < block_rows && mm < block_cols);
- }
-
- DynamicMatrix<K> convert_to_dynamic_matrix() const
- {
- DynamicMatrix<K> ret(num_rows, num_cols, 0.);
- for (size_t jj = 0; jj < num_blocks; ++jj) {
- const size_t row_offset = jj * block_rows;
- const size_t col_offset = jj * block_cols;
- for (size_t rr = 0; rr < block_rows; ++rr)
- for (size_t cc = 0; cc < block_cols; ++cc)
- ret[row_offset + rr][col_offset + cc] = block(jj)[rr][cc];
- } // jj
- return ret;
- }
-
- template <class CharType, class CharTraits>
- friend std::basic_ostream<CharType, CharTraits>& operator<<(std::basic_ostream<CharType, CharTraits>& out,
- const ThisType& mat)
- {
- return output_matrix(out, mat);
- } // ... operator<<(...)
-
-private:
- FieldVector<BlockType, num_blocks> backend_;
-};
-
-
-template <class K, int N, int M>
-struct MatrixAbstraction<Dune::XT::Common::FieldMatrix<K, N, M>>
-{
- typedef Dune::XT::Common::FieldMatrix<K, N, M> MatrixType;
- typedef typename Dune::FieldTraits<K>::field_type ScalarType;
- typedef typename Dune::FieldTraits<K>::real_type RealType;
- typedef ScalarType S;
- typedef RealType R;
- template <size_t rows = N, size_t cols = M, class FieldType = K>
- using MatrixTypeTemplate = Dune::XT::Common::FieldMatrix<FieldType, rows, cols>;
-
- static const bool is_matrix = true;
-
- static const bool has_static_size = true;
-
- static const size_t static_rows = N;
-
- static const size_t static_cols = M;
-
- static const constexpr StorageLayout storage_layout = StorageLayout::dense_row_major;
-
- static constexpr bool has_ostream = true;
-
- template <class SparsityPatternType = FullPattern>
- static inline MatrixType create(const size_t rows,
- const size_t cols,
- const ScalarType& val = suitable_default<ScalarType>::value(),
- const SparsityPatternType& /*pattern*/ = SparsityPatternType())
- {
- return MatrixType(rows, cols, val);
- }
-
- template <class SparsityPatternType = FullPattern>
- static inline std::unique_ptr<MatrixType> make_unique(const size_t rows,
- const size_t cols,
- const ScalarType& val = suitable_default<ScalarType>::value(),
- const SparsityPatternType& /*pattern*/ = SparsityPatternType())
- {
- return std::make_unique<MatrixType>(rows, cols, val);
- }
-
- static constexpr size_t rows(const MatrixType& /*mat*/)
- {
- return N;
- }
-
- static constexpr size_t cols(const MatrixType& /*mat*/)
- {
- return M;
- }
-
- static inline void set_entry(MatrixType& mat, const size_t row, const size_t col, const ScalarType& val)
- {
- mat[row][col] = val;
- }
-
- static inline ScalarType get_entry(const MatrixType& mat, const size_t row, const size_t col)
- {
- return mat[row][col];
- }
-
- static inline void add_to_entry(MatrixType& mat, const size_t row, const size_t col, const ScalarType& val)
- {
- mat[row][col] += val;
- }
-
-
- static inline ScalarType* data(MatrixType& mat)
- {
- return &(mat[0][0]);
- }
-
- static inline const ScalarType* data(const MatrixType& mat)
- {
- return &(mat[0][0]);
- }
-};
-
-template <class K, size_t num_blocks, size_t block_rows, size_t block_cols>
-struct MatrixAbstraction<Dune::XT::Common::BlockedFieldMatrix<K, num_blocks, block_rows, block_cols>>
-{
- typedef Dune::XT::Common::BlockedFieldMatrix<K, num_blocks, block_rows, block_cols> MatrixType;
- typedef typename Dune::FieldTraits<K>::field_type ScalarType;
- typedef typename Dune::FieldTraits<K>::real_type RealType;
- typedef ScalarType S;
- typedef RealType R;
-
- static const bool is_matrix = true;
-
- static const bool has_static_size = true;
-
- static const size_t static_rows = MatrixType::num_rows;
-
- static const size_t static_cols = MatrixType::num_cols;
-
- template <size_t rows = static_rows, size_t cols = static_cols, class FieldType = K>
- using MatrixTypeTemplate = Dune::XT::Common::BlockedFieldMatrix<FieldType, rows / block_rows, block_rows, block_cols>;
-
- static const constexpr StorageLayout storage_layout = StorageLayout::other;
-
- static constexpr bool has_ostream = true;
-
- template <class SparsityPatternType = FullPattern>
- static inline MatrixType create(const size_t rows,
- const size_t cols,
- const ScalarType& val = suitable_default<ScalarType>::value(),
- const SparsityPatternType& /*pattern*/ = SparsityPatternType())
- {
- return MatrixType(rows, cols, val);
- }
-
- template <class SparsityPatternType = FullPattern>
- static inline std::unique_ptr<MatrixType> make_unique(const size_t rows,
- const size_t cols,
- const ScalarType& val = suitable_default<ScalarType>::value(),
- const SparsityPatternType& /*pattern*/ = SparsityPatternType())
- {
- return std::make_unique<MatrixType>(rows, cols, val);
- }
-
- static constexpr size_t rows(const MatrixType& /*mat*/)
- {
- return static_rows;
- }
-
- static constexpr size_t cols(const MatrixType& /*mat*/)
- {
- return static_cols;
- }
-
- static inline void set_entry(MatrixType& mat, const size_t row, const size_t col, const ScalarType& val)
- {
- mat.set_entry(row, col, val);
- }
-
- static inline ScalarType get_entry(const MatrixType& mat, const size_t row, const size_t col)
- {
- return mat.get_entry(row, col);
- }
-
- static inline void add_to_entry(MatrixType& mat, const size_t row, const size_t col, const ScalarType& val)
- {
- mat.add_to_entry(row, col, val);
- }
-
- static inline ScalarType* data(MatrixType& /*mat*/)
- {
- DUNE_THROW(InvalidStateException, "Do not call me if storage layout is not dense!");
- return nullptr;
- }
-
- static inline const ScalarType* data(const MatrixType& /*mat*/)
- {
- DUNE_THROW(InvalidStateException, "Do not call me if storage_layout is not dense!");
- return nullptr;
- }
-};
-
-
-template <class M>
-typename std::enable_if<is_matrix<M>::value && MatrixAbstraction<M>::has_static_size,
- std::unique_ptr<FieldMatrix<typename MatrixAbstraction<M>::S,
- MatrixAbstraction<M>::static_rows,
- MatrixAbstraction<M>::static_cols>>>::type
-make_field_container_ptr(const M& mat)
-{
- static const size_t rows = MatrixAbstraction<M>::static_rows;
- static const size_t cols = MatrixAbstraction<M>::static_cols;
- auto ret = std::make_unique<FieldMatrix<typename MatrixAbstraction<M>::S, rows, cols>>;
- for (size_t ii = 0; ii < rows; ++ii)
- for (size_t jj = 0; jj < cols; ++jj)
- (*ret)[ii][jj] = get_matrix_entry(mat, ii, jj);
- return std::move(ret);
-}
-
-
-template <class M>
-typename std::enable_if<is_matrix<M>::value && MatrixAbstraction<M>::has_static_size,
- FieldMatrix<typename MatrixAbstraction<M>::S,
- MatrixAbstraction<M>::static_rows,
- MatrixAbstraction<M>::static_cols>>::type
-make_field_container(const M& mat)
-{
- static const size_t rows = MatrixAbstraction<M>::static_rows;
- static const size_t cols = MatrixAbstraction<M>::static_cols;
- FieldMatrix<typename MatrixAbstraction<M>::S, rows, cols> ret;
- for (size_t ii = 0; ii < rows; ++ii)
- for (size_t jj = 0; jj < cols; ++jj)
- ret[ii][jj] = get_matrix_entry(mat, ii, jj);
- return ret;
-}
-
-
-template <class K, int ROWS, int COLS>
-FieldMatrix<K, ROWS, COLS> make_field_container(Dune::FieldMatrix<K, ROWS, COLS>&& vec)
-{
- return std::move(vec);
-}
-
-
-template <class K, int ROWS, int COLS>
-FieldMatrix<K, ROWS, COLS> make_field_container(FieldMatrix<K, ROWS, COLS>&& vec)
-{
- return std::move(vec);
-}
-
-
-template <class K, int ROWS, int COLS>
-typename std::enable_if<is_arithmetic<K>::value && !is_complex<K>::value, FieldMatrix<K, ROWS, COLS>>::type
-real(const Dune::FieldMatrix<K, ROWS, COLS>& real_mat)
-{
- return real_mat;
-}
-
-template <class K, int ROWS, int COLS>
-typename std::enable_if<is_arithmetic<K>::value && !is_complex<K>::value, FieldMatrix<K, ROWS, COLS>>::type
-real(const FieldMatrix<K, ROWS, COLS>& real_mat)
-{
- return real_mat;
-}
-
-template <class K, int ROWS, int COLS>
-typename std::enable_if<is_arithmetic<K>::value && !is_complex<K>::value, FieldMatrix<K, ROWS, COLS>>::type
-real(Dune::FieldMatrix<K, ROWS, COLS>&& real_mat)
-{
- return std::move(real_mat);
-}
-
-template <class K, int ROWS, int COLS>
-typename std::enable_if<is_arithmetic<K>::value && !is_complex<K>::value, FieldMatrix<K, ROWS, COLS>>::type
-real(FieldMatrix<K, ROWS, COLS>&& real_mat)
-{
- return std::move(real_mat);
-}
-
-template <class K, int ROWS, int COLS>
-typename std::enable_if<is_complex<K>::value, FieldMatrix<real_t<K>, ROWS, COLS>>::type
-real(const Dune::FieldMatrix<K, ROWS, COLS>& complex_mat)
-{
- FieldMatrix<real_t<K>, ROWS, COLS> real_mat;
- for (size_t ii = 0; ii < ROWS; ++ii)
- for (size_t jj = 0; jj < COLS; ++jj)
- real_mat[ii][jj] = complex_mat[ii][jj].real();
- return real_mat;
-}
-
-template <class K, int ROWS, int COLS>
-typename std::enable_if<is_complex<K>::value, FieldMatrix<real_t<K>, ROWS, COLS>>::type
-real(const FieldMatrix<K, ROWS, COLS>& complex_mat)
-{
- FieldMatrix<real_t<K>, ROWS, COLS> real_mat;
- for (size_t ii = 0; ii < ROWS; ++ii)
- for (size_t jj = 0; jj < COLS; ++jj)
- real_mat[ii][jj] = complex_mat[ii][jj].real();
- return real_mat;
-}
-
-
-template <class K, int ROWS, int COLS>
-typename std::enable_if<is_arithmetic<K>::value && !is_complex<K>::value, FieldMatrix<K, ROWS, COLS>>::type
-imag(const Dune::FieldMatrix<K, ROWS, COLS>& /*real_mat*/)
-{
- return FieldMatrix<K, ROWS, COLS>(0);
-}
-
-template <class K, int ROWS, int COLS>
-typename std::enable_if<is_arithmetic<K>::value && !is_complex<K>::value, FieldMatrix<K, ROWS, COLS>>::type
-imag(const FieldMatrix<K, ROWS, COLS>& /*real_mat*/)
-{
- return FieldMatrix<K, ROWS, COLS>(0);
-}
-
-template <class K, int ROWS, int COLS>
-typename std::enable_if<is_complex<K>::value, FieldMatrix<real_t<K>, ROWS, COLS>>::type
-imag(const Dune::FieldMatrix<K, ROWS, COLS>& complex_mat)
-{
- FieldMatrix<real_t<K>, ROWS, COLS> real_mat;
- for (size_t ii = 0; ii < ROWS; ++ii)
- for (size_t jj = 0; jj < COLS; ++jj)
- real_mat[ii][jj] = complex_mat[ii][jj].imag();
- return real_mat;
-}
-
-template <class K, int ROWS, int COLS>
-typename std::enable_if<is_complex<K>::value, FieldMatrix<real_t<K>, ROWS, COLS>>::type
-imag(const FieldMatrix<K, ROWS, COLS>& complex_mat)
-{
- FieldMatrix<real_t<K>, ROWS, COLS> real_mat;
- for (size_t ii = 0; ii < ROWS; ++ii)
- for (size_t jj = 0; jj < COLS; ++jj)
- real_mat[ii][jj] = complex_mat[ii][jj].imag();
- return real_mat;
-}
-
-
-} // namespace Common
-} // namespace XT
-
-template <class K, int L_ROWS, int L_COLS, int R_COLS>
-void rightmultiply(Dune::FieldMatrix<K, L_ROWS, R_COLS>& ret,
- const Dune::FieldMatrix<K, L_ROWS, L_COLS>& left,
- const Dune::FieldMatrix<K, L_COLS, R_COLS>& right)
-{
- for (size_t ii = 0; ii < L_ROWS; ++ii) {
- for (size_t jj = 0; jj < R_COLS; ++jj) {
- ret[ii][jj] = 0.;
- for (size_t kk = 0; kk < L_COLS; ++kk)
- ret[ii][jj] += left[ii][kk] * right[kk][jj];
- }
- }
-}
-
-// versions that do not allocate matrices on the stack (for large matrices)
-template <class K, int L_ROWS, int L_COLS, int R_COLS>
-std::unique_ptr<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>
-operator*(const std::unique_ptr<Dune::FieldMatrix<K, L_ROWS, L_COLS>>& left,
- const Dune::FieldMatrix<K, L_COLS, R_COLS>& right)
-{
- auto ret = std::make_unique<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>();
- rightmultiply(*ret, *left, right);
- return ret;
-}
-
-template <class K, int L_ROWS, int L_COLS, int R_COLS>
-std::unique_ptr<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>
-operator*(const Dune::FieldMatrix<K, L_ROWS, L_COLS>& left,
- const std::unique_ptr<Dune::FieldMatrix<K, L_COLS, R_COLS>>& right)
-{
- auto ret = std::make_unique<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>();
- rightmultiply(*ret, left, *right);
- return ret;
-}
-
-template <class K, int L_ROWS, int L_COLS, int R_COLS>
-std::unique_ptr<Dune::XT::Common::FieldMatrix<K, L_ROWS, R_COLS>>
-operator*(const std::unique_ptr<Dune::FieldMatrix<K, L_ROWS, L_COLS>>& left,
- const std::unique_ptr<Dune::FieldMatrix<K, L_COLS, R_COLS>>& right)
-{
- return left * *right;
-}
-
-
-} // namespace Dune
+#if DUNE_VERSION_GTE(DUNE_COMMON, 2, 7)
+# include "fmatrix-2.7.hh"
+#else
+# include "fmatrix-2.6.hh"
+#endif
#endif // DUNE_XT_COMMON_FMATRIX_HH
diff --git a/dune/xt/common/fvector-2.6.hh b/dune/xt/common/fvector-2.6.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4f34056417a16658b6a36506beb720818008b20a
--- /dev/null
+++ b/dune/xt/common/fvector-2.6.hh
@@ -0,0 +1,856 @@
+// This file is part of the dune-xt project:
+// https://github.com/dune-community/dune-xt
+// Copyright 2009-2018 dune-xt developers and contributors. All rights reserved.
+// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+// or GPL-2.0+ (http://opensource.org/licenses/gpl-license)
+// with "runtime exception" (http://www.dune-project.org/license.html)
+// Authors:
+// Felix Schindler (2014 - 2017)
+// René Fritze (2014 - 2019)
+// Tim Keil (2018)
+// Tobias Leibner (2014, 2016 - 2019)
+
+#ifndef DUNE_XT_COMMON_FVECTOR_26_HH
+#define DUNE_XT_COMMON_FVECTOR_26_HH
+
+#include <functional>
+#include <initializer_list>
+#include <numeric>
+#include <type_traits>
+#include <vector>
+
+#include <boost/functional/hash.hpp>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/promotiontraits.hh>
+
+#include <dune/xt/common/debug.hh>
+#include <dune/xt/common/densevector.hh>
+#include <dune/xt/common/exceptions.hh>
+#include <dune/xt/common/float_cmp.hh>
+#include <dune/xt/common/vector.hh>
+#include <dune/xt/common/type_traits.hh>
+
+namespace Dune {
+namespace XT {
+namespace Common {
+
+
+/**
+ * \todo We need to implement all operators from the base which return the base, to rather return ourselfes!
+ */
+template <class K, int SIZE>
+class FieldVector : public Dune::FieldVector<K, SIZE>
+{
+ static_assert(SIZE >= 0, "Really?");
+
+ typedef Dune::FieldVector<K, SIZE> BaseType;
+ typedef FieldVector<K, SIZE> ThisType;
+
+public:
+ FieldVector(const K& kk = suitable_default<K>::value())
+ : BaseType(kk)
+ {}
+
+ FieldVector(const ThisType& other) = default;
+
+ FieldVector(const BaseType& other)
+ : BaseType(other)
+ {}
+
+ /* FieldMatrix< K, 1, 1 > is convertible to K, which in turn is convertible to FieldVector< K, 1 >. Without the
+ * following constructor, this leads to an "ambiguous constructor" error (candidates are copy constructor and
+ * constructor taking a K) */
+ template <class Type = K>
+ FieldVector(const typename std::enable_if<SIZE == 1 && std::is_same<K, Type>::value,
+ typename Dune::FieldMatrix<K, 1, 1>>::type& mat)
+ : BaseType(mat[0][0])
+ {}
+
+ FieldVector(const std::vector<K>& vec)
+ : BaseType()
+ {
+#ifndef NDEBUG
+ if (vec.size() != SIZE)
+ DUNE_THROW(Exceptions::wrong_input_given,
+ "You are trying to construct a FieldVector< ..., " << SIZE << " > (of "
+ << "static size) from a vector of size "
+ << vec.size() << "!");
+#endif // NDEBUG
+ for (size_t ii = 0; ii < SIZE; ++ii)
+ this->operator[](ii) = vec[ii];
+ } // FieldVector(...)
+
+ FieldVector(std::initializer_list<K> list)
+ : BaseType()
+ {
+#ifndef NDEBUG
+ if (list.size() != SIZE)
+ DUNE_THROW(Exceptions::wrong_input_given,
+ "You are trying to construct a FieldVector< ..., " << SIZE << " > (of "
+ << "static size) from a list of size "
+ << list.size() << "!");
+#endif // NDEBUG
+ size_t ii = 0;
+ for (auto element : list)
+ this->operator[](ii++) = element;
+ } // FieldVector(...)
+
+ template <class C>
+ FieldVector(
+ const DenseVector<C>& x,
+ typename std::enable_if<IsFieldVectorSizeCorrect<C, SIZE>::value
+ && std::is_convertible<typename DenseVector<C>::value_type, K>::value>::type* /*dummy*/
+ = nullptr)
+ : BaseType(x)
+ {}
+
+ operator std::vector<K>() const
+ {
+ std::vector<K> ret(SIZE);
+ for (size_t ii = 0; ii < SIZE; ++ii)
+ ret[ii] = this->operator[](ii);
+ return ret;
+ }
+
+ operator Dune::FieldMatrix<K, 1, SIZE>() const
+ {
+ Dune::FieldMatrix<K, 1, SIZE> ret;
+ ret[0] = *this;
+ return ret;
+ }
+
+ template <int S = SIZE>
+ operator typename std::enable_if<(S == SIZE) && (SIZE != 1), Dune::FieldMatrix<K, S, 1>>::type() const
+ {
+ Dune::FieldMatrix<K, SIZE, 1> ret;
+ for (size_t ii = 0; ii < SIZE; ++ii)
+ ret[ii][0] = this->operator[](ii);
+ return ret;
+ }
+
+ operator std::array<K, SIZE>() const
+ {
+ std::array<K, SIZE> ret;
+ for (size_t ii = 0; ii < SIZE; ++ii)
+ ret[ii] = this->operator[](ii);
+ return ret;
+ }
+
+ template <class K_Other>
+ typename multiplication_promotion<K, K_Other>::type operator*(const Dune::FieldVector<K_Other, SIZE>& other) const
+ {
+ typename multiplication_promotion<K, K_Other>::type result(0.);
+ for (size_t ii = 0; ii < SIZE; ++ii)
+ result += (*this)[ii] * other[ii];
+ return result;
+ }
+
+ template <int R>
+ typename std::enable_if<SIZE != 1 && R == SIZE, K>::type operator*(const Dune::FieldMatrix<K, R, 1>& mat) const
+ {
+ K ret = 0;
+ for (size_t ii = 0; ii < SIZE; ++ii)
+ ret += (*this)[ii] * mat[ii][0];
+ return ret;
+ }
+
+ //! This op is not redundant
+ ThisType operator*(const K& scalar) const
+ {
+ ThisType ret(*this);
+ ret *= scalar;
+ return ret;
+ }
+
+ ThisType operator/(const K& scalar) const
+ {
+ ThisType ret(*this);
+ ret /= scalar;
+ return ret;
+ }
+}; // class FieldVector
+
+
+template <class K, size_t block_num, size_t size_block>
+class BlockedFieldVector
+{
+ using ThisType = BlockedFieldVector;
+
+public:
+ static constexpr size_t num_blocks = block_num;
+ static constexpr size_t block_size = size_block;
+ static constexpr size_t static_size = num_blocks * block_size;
+ using VectorType = Dune::FieldVector<K, static_size>;
+ using XtVectorType = FieldVector<K, static_size>;
+ using BlockType = FieldVector<K, block_size>;
+
+ BlockedFieldVector(const K& val = K(0.))
+ : backend_(BlockType(val))
+ {}
+
+ template <class OtherVectorType>
+ BlockedFieldVector(const OtherVectorType& other,
+ std::enable_if_t<is_vector<OtherVectorType>::value, int> /*dummy*/ = 0)
+ {
+ *this = other;
+ }
+
+ BlockedFieldVector(const BlockType& block)
+ : backend_(block)
+ {}
+
+ BlockedFieldVector(const std::initializer_list<K>& init_list)
+ {
+ assert(init_list.size() == static_size || init_list.size() == num_blocks);
+ if (init_list.size() == static_size) {
+ *this = VectorType(init_list);
+ } else {
+ size_t jj = 0;
+ for (const auto& block_init : init_list)
+ block(jj++) = block_init;
+ }
+ }
+
+ template <class OtherVectorType>
+ std::enable_if_t<is_vector<OtherVectorType>::value, ThisType>& operator=(const OtherVectorType& other)
+ {
+ assert(other.size() == static_size);
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ for (size_t ii = 0; ii < block_size; ++ii)
+ backend_[jj][ii] = other[jj * block_size + ii];
+ return *this;
+ }
+
+ size_t size() const
+ {
+ return static_size;
+ }
+
+ K& operator[](const size_t ii)
+ {
+ assert(ii < static_size);
+ return backend_[ii / block_size][ii % block_size];
+ }
+
+ const K& operator[](const size_t ii) const
+ {
+ assert(ii < static_size);
+ return backend_[ii / block_size][ii % block_size];
+ }
+
+ K& get_entry(const size_t jj, const size_t ii)
+ {
+ assert(jj < num_blocks);
+ assert(ii < block_size);
+ return backend_[jj][ii];
+ }
+
+ const K& get_entry(const size_t jj, const size_t ii) const
+ {
+ assert(jj < num_blocks);
+ assert(ii < block_size);
+ return backend_[jj][ii];
+ }
+
+ BlockType& block(const size_t jj)
+ {
+ assert(jj < num_blocks);
+ return backend_[jj];
+ }
+
+ const BlockType& block(const size_t jj) const
+ {
+ assert(jj < num_blocks);
+ return backend_[jj];
+ }
+
+ K one_norm() const
+ {
+ K ret(0);
+ for (const auto& block_vec : backend_)
+ for (const auto& entry : block_vec)
+ ret += std::abs(entry);
+ return ret;
+ }
+
+ K two_norm() const
+ {
+ return std::sqrt(two_norm2());
+ }
+
+ K two_norm2() const
+ {
+ K ret(0);
+ for (const auto& block_vec : backend_)
+ for (const auto& entry : block_vec)
+ ret += std::pow(entry, 2);
+ return ret;
+ }
+
+ operator XtVectorType() const
+ {
+ XtVectorType ret(0.);
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ for (size_t ii = 0; ii < block_size; ++ii)
+ ret[jj * block_size + ii] = backend_[jj][ii];
+ return ret;
+ }
+
+ K* data()
+ {
+ return &(backend_[0][0]);
+ }
+
+ const K* data() const
+ {
+ return &(backend_[0][0]);
+ }
+
+ K* begin()
+ {
+ return data();
+ }
+
+ const K* begin() const
+ {
+ return data();
+ }
+
+ K* end()
+ {
+ return &(backend_[num_blocks - 1][block_size]);
+ }
+
+ const K* end() const
+ {
+ return &(backend_[num_blocks - 1][block_size]);
+ }
+
+ ThisType& operator*=(const K& val)
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ backend_[jj] *= val;
+ return *this;
+ }
+
+ ThisType operator*(const K& val) const
+ {
+ auto ret = *this;
+ ret *= val;
+ return ret;
+ }
+
+ K operator*(const ThisType& other) const
+ {
+ return std::inner_product(begin(), end(), other.begin(), 0.);
+ }
+
+ ThisType& operator+=(const ThisType& other)
+ {
+ backend_ += other.backend_;
+ return *this;
+ }
+
+ ThisType operator+(const ThisType& other) const
+ {
+ auto ret = *this;
+ ret += other;
+ return ret;
+ }
+
+ ThisType& operator-=(const ThisType& other)
+ {
+ backend_ -= other.backend_;
+ return *this;
+ }
+
+ ThisType operator-(const ThisType& other) const
+ {
+ auto ret = *this;
+ ret -= other;
+ return ret;
+ }
+
+private:
+ FieldVector<BlockType, num_blocks> backend_;
+};
+
+
+//! this allows to set the init value of the FieldVector at compile time
+template <class K, int SIZE, K value>
+class ValueInitFieldVector : public Dune::XT::Common::FieldVector<K, SIZE>
+{
+ typedef Dune::XT::Common::FieldVector<K, SIZE> BaseType;
+
+public:
+ ValueInitFieldVector()
+ : BaseType(value)
+ {}
+}; // class ValueInitFieldVector
+
+
+//! struct to be used as comparison function e.g. in a std::map<FieldVector<...>, ..., FieldVectorLess>
+struct FieldVectorLess
+{
+ template <class FieldType, int dimDomain>
+ bool operator()(const Dune::FieldVector<FieldType, dimDomain>& a,
+ const Dune::FieldVector<FieldType, dimDomain>& b) const
+ {
+ for (size_t dd = 0; dd < dimDomain; ++dd) {
+ if (a[dd] < b[dd])
+ return true;
+ else if (a[dd] > b[dd])
+ return false;
+ }
+ return false;
+ }
+};
+
+struct FieldVectorFloatLess
+{
+ template <class FieldType, int dimDomain>
+ bool operator()(const Dune::FieldVector<FieldType, dimDomain>& a,
+ const Dune::FieldVector<FieldType, dimDomain>& b) const
+ {
+ for (size_t dd = 0; dd < dimDomain; ++dd) {
+ if (XT::Common::FloatCmp::lt(a[dd], b[dd]))
+ return true;
+ else if (XT::Common::FloatCmp::gt(a[dd], b[dd]))
+ return false;
+ }
+ return false;
+ }
+};
+
+
+//! Specialization of VectorAbstraction for Dune::XT::Common::FieldVector
+template <class K, int SIZE>
+struct VectorAbstraction<Dune::XT::Common::FieldVector<K, SIZE>>
+ : public internal::VectorAbstractionBase<Dune::XT::Common::FieldVector<K, SIZE>, K>
+ , public internal::HasSubscriptOperatorForVectorAbstraction<Dune::XT::Common::FieldVector<K, SIZE>,
+ typename Dune::FieldTraits<K>::field_type>
+{
+ static const bool has_static_size = true;
+ static const size_t static_size = SIZE;
+ static const bool is_contiguous = true;
+
+ template <size_t SZ = SIZE, class Field = K>
+ using VectorTypeTemplate = Dune::XT::Common::FieldVector<Field, SZ>;
+
+ template <size_t SZ = SIZE>
+ static inline VectorTypeTemplate<SZ> create(const size_t sz, const K& val = suitable_default<K>::value())
+ {
+ if (sz != SZ)
+ DUNE_THROW(Exceptions::wrong_input_given,
+ "You are trying to construct a FieldVector< ..., " << SZ << " > (of "
+ << "static size) with " << sz << " elements!");
+ return VectorTypeTemplate<SZ>(val);
+ }
+};
+
+
+//! Specialization of VectorAbstraction for Dune::XT::Common::BlockedFieldVector
+template <class K, size_t num_blocks, size_t block_size>
+struct VectorAbstraction<Dune::XT::Common::BlockedFieldVector<K, num_blocks, block_size>>
+ : public internal::VectorAbstractionBase<Dune::XT::Common::BlockedFieldVector<K, num_blocks, block_size>, K>
+ , public internal::HasSubscriptOperatorForVectorAbstraction<
+ Dune::XT::Common::BlockedFieldVector<K, num_blocks, block_size>,
+ typename Dune::FieldTraits<K>::field_type>
+{
+ using VectorType = Dune::XT::Common::BlockedFieldVector<K, num_blocks, block_size>;
+ static constexpr bool has_static_size = true;
+ static constexpr size_t static_size = VectorType::static_size;
+ static constexpr bool is_contiguous = true;
+
+ template <size_t SZ = static_size>
+ static inline VectorType create(const size_t sz, const K& val = suitable_default<K>::value())
+ {
+ static_assert(SZ == static_size, "Creation of Vector with different size not implemented!");
+ if (sz != SZ)
+ DUNE_THROW(Exceptions::wrong_input_given,
+ "You are trying to construct a FieldVector< ..., " << SZ << " > (of "
+ << "static size) with " << sz << " elements!");
+ return VectorType(val);
+ }
+};
+
+
+template <class V>
+typename std::enable_if<
+ is_vector<V>::value && VectorAbstraction<V>::has_static_size,
+ std::unique_ptr<FieldVector<typename VectorAbstraction<V>::S, VectorAbstraction<V>::static_size>>>::type
+make_field_container_ptr(const V& vec)
+{
+ auto ret = std::make_unique<FieldVector<typename VectorAbstraction<V>::S, VectorAbstraction<V>::static_size>>;
+ for (size_t ii = 0; ii < ret->size(); ++ii)
+ (*ret)[ii] = vec[ii];
+ return std::move(ret);
+}
+
+
+template <class V>
+typename std::enable_if<is_vector<V>::value && VectorAbstraction<V>::has_static_size,
+ FieldVector<typename VectorAbstraction<V>::S, VectorAbstraction<V>::static_size>>::type
+make_field_container(const V& vec)
+{
+ FieldVector<typename VectorAbstraction<V>::S, VectorAbstraction<V>::static_size> ret;
+ for (size_t ii = 0; ii < ret.size(); ++ii)
+ ret[ii] = vec[ii];
+ return ret;
+}
+
+
+template <class K, int SIZE>
+FieldVector<K, SIZE> make_field_container(Dune::FieldVector<K, SIZE>&& vec)
+{
+ return std::move(vec);
+}
+
+
+template <class K, int SIZE>
+FieldVector<K, SIZE> make_field_container(FieldVector<K, SIZE>&& vec)
+{
+ return std::move(vec);
+}
+
+template <class K>
+FieldVector<K, 3> cross_product(const FieldVector<K, 3>& u, const FieldVector<K, 3>& v)
+{
+ FieldVector<K, 3> ret;
+ ret[0] = u[1] * v[2] - u[2] * v[1];
+ ret[1] = u[2] * v[0] - u[0] * v[2];
+ ret[2] = u[0] * v[1] - u[1] * v[0];
+ return ret;
+}
+
+
+namespace internal {
+
+
+/**
+ * Most likely, you do not want to use this class directly, but \sa hstack instead.
+ */
+template <class K>
+struct hstack_decay
+{
+ static_assert(is_arithmetic<K>::value, "");
+ using type = K;
+};
+
+template <class K, int SIZE>
+struct hstack_decay<Dune::FieldVector<K, SIZE>>
+{
+ using type = Dune::XT::Common::FieldVector<K, SIZE>;
+};
+
+template <class K, int SIZE>
+struct hstack_decay<const Dune::FieldVector<K, SIZE>&> : public hstack_decay<Dune::FieldVector<K, SIZE>>
+{};
+
+template <class K, int SIZE>
+struct hstack_decay<Dune::XT::Common::FieldVector<K, SIZE>> : public hstack_decay<Dune::FieldVector<K, SIZE>>
+{};
+
+template <class K, int SIZE>
+struct hstack_decay<const Dune::XT::Common::FieldVector<K, SIZE>&> : public hstack_decay<Dune::FieldVector<K, SIZE>>
+{};
+
+
+/**
+ * Most likely, you do not want to use this, but \sa hstack instead.
+ */
+template <class T>
+using hstack_decay_t = typename hstack_decay<T>::type;
+
+
+/**
+ * Most likely, you do not want to use this class directly, but \sa hstack instead.
+ */
+template <class... Vectors>
+struct hstack_helper;
+
+template <class L, class R, class... Vectors>
+struct hstack_helper<L, R, Vectors...>
+{
+ using type =
+ typename hstack_helper<typename hstack_helper<hstack_decay_t<L>, hstack_decay_t<R>>::type, Vectors...>::type;
+};
+
+template <class K, int SIZE>
+struct hstack_helper<Dune::FieldVector<K, SIZE>>
+{
+ using type = hstack_decay_t<Dune::FieldVector<K, SIZE>>;
+};
+
+template <class K, int SIZE>
+struct hstack_helper<Dune::XT::Common::FieldVector<K, SIZE>> : public hstack_helper<Dune::FieldVector<K, SIZE>>
+{};
+
+template <class KL, class KR, int r>
+struct hstack_helper<KL, Dune::FieldVector<KR, r>>
+{
+ using type = hstack_decay_t<Dune::FieldVector<typename PromotionTraits<KL, KR>::PromotedType, 1 + r>>;
+};
+
+template <class KL, class KR, int r>
+struct hstack_helper<KL, Dune::XT::Common::FieldVector<KR, r>> : public hstack_helper<KL, Dune::FieldVector<KR, r>>
+{};
+
+template <class KL, int l, class KR>
+struct hstack_helper<Dune::FieldVector<KL, l>, KR>
+{
+ using type = hstack_decay_t<Dune::FieldVector<typename PromotionTraits<KL, KR>::PromotedType, l + 1>>;
+};
+
+template <class KL, int l, class KR>
+struct hstack_helper<Dune::XT::Common::FieldVector<KL, l>, KR> : public hstack_helper<Dune::FieldVector<KL, l>, KR>
+{};
+
+template <class KL, int l, class KR, int r>
+struct hstack_helper<Dune::FieldVector<KL, l>, Dune::FieldVector<KR, r>>
+{
+ using type = hstack_decay_t<Dune::FieldVector<typename PromotionTraits<KL, KR>::PromotedType, l + r>>;
+};
+
+template <class KL, int l, class KR, int r>
+struct hstack_helper<Dune::XT::Common::FieldVector<KL, l>, Dune::FieldVector<KR, r>>
+ : public hstack_helper<Dune::FieldVector<KL, l>, Dune::FieldVector<KR, r>>
+{};
+
+template <class KL, int l, class KR, int r>
+struct hstack_helper<Dune::FieldVector<KL, l>, Dune::XT::Common::FieldVector<KR, r>>
+ : public hstack_helper<Dune::FieldVector<KL, l>, Dune::FieldVector<KR, r>>
+{};
+
+template <class KL, int l, class KR, int r>
+struct hstack_helper<Dune::XT::Common::FieldVector<KL, l>, Dune::XT::Common::FieldVector<KR, r>>
+ : public hstack_helper<Dune::FieldVector<KL, l>, Dune::FieldVector<KR, r>>
+{};
+
+// not sure why this is required, would've hoped the decay above would take care of it
+template <class KL, int l, class KR, int r>
+struct hstack_helper<Dune::FieldVector<KL, l>, const Dune::FieldVector<KR, r>&>
+ : public hstack_helper<Dune::FieldVector<KL, l>, Dune::FieldVector<KR, r>>
+{};
+
+// not sure why this is required, would've hoped the decay above would take care of it
+template <class KL, int l, class KR, int r>
+struct hstack_helper<Dune::FieldVector<KL, l>, const Dune::XT::Common::FieldVector<KR, r>&>
+ : public hstack_helper<Dune::FieldVector<KL, l>, Dune::FieldVector<KR, r>>
+{};
+
+
+} // namespace internal
+
+
+/// \brief Specialization of \sa hstack to stop the recursion.
+template <class K, int SIZE>
+Dune::XT::Common::FieldVector<K, SIZE> hstack(Dune::FieldVector<K, SIZE>&& vec)
+{
+ return std::move(vec);
+}
+
+/// \brief Specialization of \sa hstack to stop the recursion.
+template <class K, int SIZE>
+Dune::XT::Common::FieldVector<K, SIZE> hstack(Dune::XT::Common::FieldVector<K, SIZE>&& vec)
+{
+ return std::move(vec);
+}
+
+/// \brief Specialization of \sa hstack to stop the recursion.
+template <class K, int SIZE>
+Dune::XT::Common::FieldVector<K, SIZE> hstack(const Dune::FieldVector<K, SIZE>& vec)
+{
+ return vec;
+}
+
+/// \brief Specialization of \sa hstack to stop the recursion.
+template <class K, int SIZE>
+Dune::XT::Common::FieldVector<K, SIZE> hstack(const Dune::XT::Common::FieldVector<K, SIZE>& vec)
+{
+ return vec;
+}
+
+
+/// \brief Specialization of \sa hstack for two arguments.
+template <class KL, class KR, int r>
+typename std::enable_if<
+ is_arithmetic<KL>::value && !is_field_vector<KL>::value,
+ typename internal::hstack_helper<internal::hstack_decay_t<KL>,
+ internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type>::type
+hstack(const KL& left, const Dune::FieldVector<KR, r>& right)
+{
+ typename internal::hstack_helper<internal::hstack_decay_t<KL>,
+ internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type ret;
+ ret[0] = left;
+ for (size_t ii = 0; ii < r; ++ii)
+ ret[ii + 1] = right[ii];
+ return ret;
+} // ... hstack<1, r>(...)
+
+/// \brief Specialization of \sa hstack for two arguments.
+template <class KL, class KR, int r>
+typename std::enable_if<
+ is_arithmetic<KL>::value && !is_field_vector<KL>::value,
+ typename internal::hstack_helper<internal::hstack_decay_t<KL>,
+ internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type>::type
+hstack(const KL& left, const Dune::XT::Common::FieldVector<KR, r>& right)
+{
+ return hstack(left, static_cast<const Dune::FieldVector<KR, r>&>(right));
+}
+
+/// \brief Specialization of \sa hstack for two arguments.
+template <class KL, int l, class KR>
+typename std::enable_if<is_arithmetic<KR>::value && !is_field_vector<KR>::value,
+ typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
+ internal::hstack_decay_t<KR>>::type>::type
+hstack(const Dune::FieldVector<KL, l>& left, const KR& right)
+{
+ typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
+ internal::hstack_decay_t<KR>>::type ret;
+ for (size_t ii = 0; ii < l; ++ii)
+ ret[ii] = left[ii];
+ ret[l] = right;
+ return ret;
+} // ... hstack<l, 1>(...)
+
+/// \brief Specialization of \sa hstack for two arguments.
+template <class KL, int l, class KR>
+typename std::enable_if<is_arithmetic<KR>::value && !is_field_vector<KR>::value,
+ typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
+ internal::hstack_decay_t<KR>>::type>::type
+hstack(const Dune::XT::Common::FieldVector<KL, l>& left, const KR& right)
+{
+ return hstack(static_cast<const Dune::FieldVector<KL, l>&>(left), right);
+}
+
+/// \brief Specialization of \sa hstack for two arguments.
+template <class KL, int l, class KR, int r>
+typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
+ internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type
+hstack(const Dune::FieldVector<KL, l>& left, const Dune::FieldVector<KR, r>& right)
+{
+ typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
+ internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type ret;
+ size_t ii = 0;
+ for (size_t jj = 0; jj < l; ++jj, ++ii)
+ ret[ii] = left[jj];
+ for (size_t jj = 0; jj < r; ++jj, ++ii)
+ ret[ii] = right[jj];
+ return ret;
+} // ... hstack<l, r>(...)
+
+/// \brief Specialization of \sa hstack for two arguments.
+template <class KL, int l, class KR, int r>
+typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
+ internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type
+hstack(const Dune::XT::Common::FieldVector<KL, l>& left, const Dune::FieldVector<KR, r>& right)
+{
+ return hstack(static_cast<const Dune::FieldVector<KL, l>&>(left), right);
+}
+
+/// \brief Specialization of \sa hstack for two arguments.
+template <class KL, int l, class KR, int r>
+typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
+ internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type
+hstack(const Dune::FieldVector<KL, l>& left, const Dune::XT::Common::FieldVector<KR, r>& right)
+{
+ return hstack(left, static_cast<const Dune::FieldVector<KR, r>&>(right));
+}
+
+/// \brief Specialization of \sa hstack for two arguments.
+template <class KL, int l, class KR, int r>
+typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
+ internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type
+hstack(const Dune::XT::Common::FieldVector<KL, l>& left, const Dune::XT::Common::FieldVector<KR, r>& right)
+{
+ return hstack(static_cast<const Dune::FieldVector<KL, l>&>(left),
+ static_cast<const Dune::FieldVector<KR, r>&>(right));
+}
+
+
+/**
+ * \brief Stacks an arbitrary number of numbers and FieldVectors horizontally into one appropriate FieldVector.
+ * \note At least one argument must be a FieldVector!
+ */
+template <class L, class R, class... Vectors>
+typename std::enable_if<
+ is_field_vector<L>::value || is_field_vector<R>::value,
+ typename internal::hstack_helper<internal::hstack_decay_t<L>, internal::hstack_decay_t<R>, Vectors...>::type>::type
+hstack(const L& left, const R& right, Vectors&&... vectors)
+{
+ return hstack(hstack(left, right), std::forward<Vectors>(vectors)...);
+}
+
+
+} // namespace Common
+} // namespace XT
+
+
+template <class K, int SIZE>
+typename std::enable_if<SIZE != 1, K>::type operator*(const Dune::FieldVector<K, SIZE>& vec,
+ const Dune::FieldMatrix<K, SIZE, 1>& mat)
+{
+ K ret = 0;
+ for (size_t ii = 0; ii < SIZE; ++ii)
+ ret += vec[ii] * mat[ii][0];
+ return ret;
+}
+
+template <class K, int SIZE>
+struct FieldTraits<XT::Common::FieldVector<K, SIZE>>
+{
+ typedef typename FieldTraits<K>::field_type field_type;
+ typedef typename FieldTraits<K>::real_type real_type;
+};
+
+
+} // namespace Dune
+namespace std {
+
+
+/**
+ * \sa http://en.cppreference.com/w/cpp/utility/hash
+ * \sa http://www.boost.org/doc/libs/1_55_0/doc/html/hash/combine.html
+ */
+template <class K, int SIZE>
+struct hash<Dune::FieldVector<K, SIZE>>
+{
+ typedef Dune::FieldVector<K, SIZE> argument_type;
+ typedef std::size_t result_type;
+
+ result_type operator()(argument_type const& s) const noexcept
+ {
+ std::size_t seed = 0;
+ for (size_t ii = 0; ii < s.size(); ++ii)
+ boost::hash_combine(seed, s[ii]);
+ return seed;
+ }
+}; // struct hash<Dune::FieldVector<...>>
+
+
+/**
+ * \sa http://en.cppreference.com/w/cpp/utility/hash
+ * \sa http://www.boost.org/doc/libs/1_55_0/doc/html/hash/combine.html
+ */
+template <class K, int SIZE>
+struct hash<Dune::XT::Common::FieldVector<K, SIZE>>
+{
+ typedef Dune::XT::Common::FieldVector<K, SIZE> argument_type;
+ typedef std::size_t result_type;
+
+ result_type operator()(argument_type const& s) const noexcept
+ {
+ std::size_t seed = 0;
+ for (size_t ii = 0; ii < s.size(); ++ii)
+ boost::hash_combine(seed, s[ii]);
+ return seed;
+ }
+}; // struct hash<Dune::XT::Common::FieldVector<...>>
+
+
+} // namespace std
+
+#endif // DUNE_XT_COMMON_FVECTOR_26_HH
diff --git a/dune/xt/common/fvector-2.7.hh b/dune/xt/common/fvector-2.7.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4bcb0807f80d95bd8b9f4e17b8328091994dab57
--- /dev/null
+++ b/dune/xt/common/fvector-2.7.hh
@@ -0,0 +1,856 @@
+// This file is part of the dune-xt project:
+// https://github.com/dune-community/dune-xt
+// Copyright 2009-2018 dune-xt developers and contributors. All rights reserved.
+// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+// or GPL-2.0+ (http://opensource.org/licenses/gpl-license)
+// with "runtime exception" (http://www.dune-project.org/license.html)
+// Authors:
+// Felix Schindler (2014 - 2017)
+// René Fritze (2014 - 2019)
+// Tim Keil (2018)
+// Tobias Leibner (2014, 2016 - 2019)
+
+#ifndef DUNE_XT_COMMON_FVECTOR_27_HH
+#define DUNE_XT_COMMON_FVECTOR_27_HH
+
+#include <functional>
+#include <initializer_list>
+#include <numeric>
+#include <type_traits>
+#include <vector>
+
+#include <boost/functional/hash.hpp>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/promotiontraits.hh>
+
+#include <dune/xt/common/debug.hh>
+#include <dune/xt/common/densevector.hh>
+#include <dune/xt/common/exceptions.hh>
+#include <dune/xt/common/float_cmp.hh>
+#include <dune/xt/common/vector.hh>
+#include <dune/xt/common/type_traits.hh>
+
+namespace Dune {
+namespace XT {
+namespace Common {
+
+
+/**
+ * \todo We need to implement all operators from the base which return the base, to rather return ourselfes!
+ */
+template <class K, int SIZE>
+class FieldVector : public Dune::FieldVector<K, SIZE>
+{
+ static_assert(SIZE >= 0, "Really?");
+
+ typedef Dune::FieldVector<K, SIZE> BaseType;
+ typedef FieldVector<K, SIZE> ThisType;
+
+public:
+ FieldVector(const K& kk = suitable_default<K>::value())
+ : BaseType(kk)
+ {}
+
+ FieldVector(const ThisType& other) = default;
+
+ FieldVector(const BaseType& other)
+ : BaseType(other)
+ {}
+
+ /* FieldMatrix< K, 1, 1 > is convertible to K, which in turn is convertible to FieldVector< K, 1 >. Without the
+ * following constructor, this leads to an "ambiguous constructor" error (candidates are copy constructor and
+ * constructor taking a K) */
+ template <class Type = K>
+ FieldVector(const typename std::enable_if<SIZE == 1 && std::is_same<K, Type>::value,
+ typename Dune::FieldMatrix<K, 1, 1>>::type& mat)
+ : BaseType(mat[0][0])
+ {}
+
+ FieldVector(const std::vector<K>& vec)
+ : BaseType()
+ {
+#ifndef NDEBUG
+ if (vec.size() != SIZE)
+ DUNE_THROW(Exceptions::wrong_input_given,
+ "You are trying to construct a FieldVector< ..., " << SIZE << " > (of "
+ << "static size) from a vector of size "
+ << vec.size() << "!");
+#endif // NDEBUG
+ for (size_t ii = 0; ii < SIZE; ++ii)
+ this->operator[](ii) = vec[ii];
+ } // FieldVector(...)
+
+ FieldVector(std::initializer_list<K> list)
+ : BaseType()
+ {
+#ifndef NDEBUG
+ if (list.size() != SIZE)
+ DUNE_THROW(Exceptions::wrong_input_given,
+ "You are trying to construct a FieldVector< ..., " << SIZE << " > (of "
+ << "static size) from a list of size "
+ << list.size() << "!");
+#endif // NDEBUG
+ size_t ii = 0;
+ for (auto element : list)
+ this->operator[](ii++) = element;
+ } // FieldVector(...)
+
+ template <class C>
+ FieldVector(
+ const DenseVector<C>& x,
+ typename std::enable_if<IsFieldVectorSizeCorrect<C, SIZE>::value
+ && std::is_convertible<typename DenseVector<C>::value_type, K>::value>::type* /*dummy*/
+ = nullptr)
+ : BaseType(x)
+ {}
+
+ operator std::vector<K>() const
+ {
+ std::vector<K> ret(SIZE);
+ for (size_t ii = 0; ii < SIZE; ++ii)
+ ret[ii] = this->operator[](ii);
+ return ret;
+ }
+
+ operator Dune::FieldMatrix<K, 1, SIZE>() const
+ {
+ Dune::FieldMatrix<K, 1, SIZE> ret;
+ ret[0] = *this;
+ return ret;
+ }
+
+ template <int S = SIZE>
+ operator typename std::enable_if<(S == SIZE) && (SIZE != 1), Dune::FieldMatrix<K, S, 1>>::type() const
+ {
+ Dune::FieldMatrix<K, SIZE, 1> ret;
+ for (size_t ii = 0; ii < SIZE; ++ii)
+ ret[ii][0] = this->operator[](ii);
+ return ret;
+ }
+
+ operator std::array<K, SIZE>() const
+ {
+ std::array<K, SIZE> ret;
+ for (size_t ii = 0; ii < SIZE; ++ii)
+ ret[ii] = this->operator[](ii);
+ return ret;
+ }
+
+ template <class K_Other>
+ typename multiplication_promotion<K, K_Other>::type operator*(const Dune::FieldVector<K_Other, SIZE>& other) const
+ {
+ typename multiplication_promotion<K, K_Other>::type result(0.);
+ for (size_t ii = 0; ii < SIZE; ++ii)
+ result += (*this)[ii] * other[ii];
+ return result;
+ }
+
+ template <int R>
+ typename std::enable_if<R == SIZE, K>::type operator*(const Dune::FieldMatrix<K, R, 1>& mat) const
+ {
+ K ret = 0;
+ for (size_t ii = 0; ii < SIZE; ++ii)
+ ret += (*this)[ii] * mat[ii][0];
+ return ret;
+ }
+
+ //! This op is not redundant
+ ThisType operator*(const K& scalar) const
+ {
+ ThisType ret(*this);
+ ret *= scalar;
+ return ret;
+ }
+
+ ThisType operator/(const K& scalar) const
+ {
+ ThisType ret(*this);
+ ret /= scalar;
+ return ret;
+ }
+}; // class FieldVector
+
+
+template <class K, size_t block_num, size_t size_block>
+class BlockedFieldVector
+{
+ using ThisType = BlockedFieldVector;
+
+public:
+ static constexpr size_t num_blocks = block_num;
+ static constexpr size_t block_size = size_block;
+ static constexpr size_t static_size = num_blocks * block_size;
+ using VectorType = Dune::FieldVector<K, static_size>;
+ using XtVectorType = FieldVector<K, static_size>;
+ using BlockType = FieldVector<K, block_size>;
+
+ BlockedFieldVector(const K& val = K(0.))
+ : backend_(BlockType(val))
+ {}
+
+ template <class OtherVectorType>
+ BlockedFieldVector(const OtherVectorType& other,
+ std::enable_if_t<is_vector<OtherVectorType>::value, int> /*dummy*/ = 0)
+ {
+ *this = other;
+ }
+
+ BlockedFieldVector(const BlockType& block)
+ : backend_(block)
+ {}
+
+ BlockedFieldVector(const std::initializer_list<K>& init_list)
+ {
+ assert(init_list.size() == static_size || init_list.size() == num_blocks);
+ if (init_list.size() == static_size) {
+ *this = VectorType(init_list);
+ } else {
+ size_t jj = 0;
+ for (const auto& block_init : init_list)
+ block(jj++) = block_init;
+ }
+ }
+
+ template <class OtherVectorType>
+ std::enable_if_t<is_vector<OtherVectorType>::value, ThisType>& operator=(const OtherVectorType& other)
+ {
+ assert(other.size() == static_size);
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ for (size_t ii = 0; ii < block_size; ++ii)
+ backend_[jj][ii] = other[jj * block_size + ii];
+ return *this;
+ }
+
+ size_t size() const
+ {
+ return static_size;
+ }
+
+ K& operator[](const size_t ii)
+ {
+ assert(ii < static_size);
+ return backend_[ii / block_size][ii % block_size];
+ }
+
+ const K& operator[](const size_t ii) const
+ {
+ assert(ii < static_size);
+ return backend_[ii / block_size][ii % block_size];
+ }
+
+ K& get_entry(const size_t jj, const size_t ii)
+ {
+ assert(jj < num_blocks);
+ assert(ii < block_size);
+ return backend_[jj][ii];
+ }
+
+ const K& get_entry(const size_t jj, const size_t ii) const
+ {
+ assert(jj < num_blocks);
+ assert(ii < block_size);
+ return backend_[jj][ii];
+ }
+
+ BlockType& block(const size_t jj)
+ {
+ assert(jj < num_blocks);
+ return backend_[jj];
+ }
+
+ const BlockType& block(const size_t jj) const
+ {
+ assert(jj < num_blocks);
+ return backend_[jj];
+ }
+
+ K one_norm() const
+ {
+ K ret(0);
+ for (const auto& block_vec : backend_)
+ for (const auto& entry : block_vec)
+ ret += std::abs(entry);
+ return ret;
+ }
+
+ K two_norm() const
+ {
+ return std::sqrt(two_norm2());
+ }
+
+ K two_norm2() const
+ {
+ K ret(0);
+ for (const auto& block_vec : backend_)
+ for (const auto& entry : block_vec)
+ ret += std::pow(entry, 2);
+ return ret;
+ }
+
+ operator XtVectorType() const
+ {
+ XtVectorType ret(0.);
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ for (size_t ii = 0; ii < block_size; ++ii)
+ ret[jj * block_size + ii] = backend_[jj][ii];
+ return ret;
+ }
+
+ K* data()
+ {
+ return &(backend_[0][0]);
+ }
+
+ const K* data() const
+ {
+ return &(backend_[0][0]);
+ }
+
+ K* begin()
+ {
+ return data();
+ }
+
+ const K* begin() const
+ {
+ return data();
+ }
+
+ K* end()
+ {
+ return &(backend_[num_blocks - 1][block_size]);
+ }
+
+ const K* end() const
+ {
+ return &(backend_[num_blocks - 1][block_size]);
+ }
+
+ ThisType& operator*=(const K& val)
+ {
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ backend_[jj] *= val;
+ return *this;
+ }
+
+ ThisType operator*(const K& val) const
+ {
+ auto ret = *this;
+ ret *= val;
+ return ret;
+ }
+
+ K operator*(const ThisType& other) const
+ {
+ return std::inner_product(begin(), end(), other.begin(), 0.);
+ }
+
+ ThisType& operator+=(const ThisType& other)
+ {
+ backend_ += other.backend_;
+ return *this;
+ }
+
+ ThisType operator+(const ThisType& other) const
+ {
+ auto ret = *this;
+ ret += other;
+ return ret;
+ }
+
+ ThisType& operator-=(const ThisType& other)
+ {
+ backend_ -= other.backend_;
+ return *this;
+ }
+
+ ThisType operator-(const ThisType& other) const
+ {
+ auto ret = *this;
+ ret -= other;
+ return ret;
+ }
+
+private:
+ FieldVector<BlockType, num_blocks> backend_;
+};
+
+
+//! this allows to set the init value of the FieldVector at compile time
+template <class K, int SIZE, K value>
+class ValueInitFieldVector : public Dune::XT::Common::FieldVector<K, SIZE>
+{
+ typedef Dune::XT::Common::FieldVector<K, SIZE> BaseType;
+
+public:
+ ValueInitFieldVector()
+ : BaseType(value)
+ {}
+}; // class ValueInitFieldVector
+
+
+//! struct to be used as comparison function e.g. in a std::map<FieldVector<...>, ..., FieldVectorLess>
+struct FieldVectorLess
+{
+ template <class FieldType, int dimDomain>
+ bool operator()(const Dune::FieldVector<FieldType, dimDomain>& a,
+ const Dune::FieldVector<FieldType, dimDomain>& b) const
+ {
+ for (size_t dd = 0; dd < dimDomain; ++dd) {
+ if (a[dd] < b[dd])
+ return true;
+ else if (a[dd] > b[dd])
+ return false;
+ }
+ return false;
+ }
+};
+
+struct FieldVectorFloatLess
+{
+ template <class FieldType, int dimDomain>
+ bool operator()(const Dune::FieldVector<FieldType, dimDomain>& a,
+ const Dune::FieldVector<FieldType, dimDomain>& b) const
+ {
+ for (size_t dd = 0; dd < dimDomain; ++dd) {
+ if (XT::Common::FloatCmp::lt(a[dd], b[dd]))
+ return true;
+ else if (XT::Common::FloatCmp::gt(a[dd], b[dd]))
+ return false;
+ }
+ return false;
+ }
+};
+
+
+//! Specialization of VectorAbstraction for Dune::XT::Common::FieldVector
+template <class K, int SIZE>
+struct VectorAbstraction<Dune::XT::Common::FieldVector<K, SIZE>>
+ : public internal::VectorAbstractionBase<Dune::XT::Common::FieldVector<K, SIZE>, K>
+ , public internal::HasSubscriptOperatorForVectorAbstraction<Dune::XT::Common::FieldVector<K, SIZE>,
+ typename Dune::FieldTraits<K>::field_type>
+{
+ static const bool has_static_size = true;
+ static const size_t static_size = SIZE;
+ static const bool is_contiguous = true;
+
+ template <size_t SZ = SIZE, class Field = K>
+ using VectorTypeTemplate = Dune::XT::Common::FieldVector<Field, SZ>;
+
+ template <size_t SZ = SIZE>
+ static inline VectorTypeTemplate<SZ> create(const size_t sz, const K& val = suitable_default<K>::value())
+ {
+ if (sz != SZ)
+ DUNE_THROW(Exceptions::wrong_input_given,
+ "You are trying to construct a FieldVector< ..., " << SZ << " > (of "
+ << "static size) with " << sz << " elements!");
+ return VectorTypeTemplate<SZ>(val);
+ }
+};
+
+
+//! Specialization of VectorAbstraction for Dune::XT::Common::BlockedFieldVector
+template <class K, size_t num_blocks, size_t block_size>
+struct VectorAbstraction<Dune::XT::Common::BlockedFieldVector<K, num_blocks, block_size>>
+ : public internal::VectorAbstractionBase<Dune::XT::Common::BlockedFieldVector<K, num_blocks, block_size>, K>
+ , public internal::HasSubscriptOperatorForVectorAbstraction<
+ Dune::XT::Common::BlockedFieldVector<K, num_blocks, block_size>,
+ typename Dune::FieldTraits<K>::field_type>
+{
+ using VectorType = Dune::XT::Common::BlockedFieldVector<K, num_blocks, block_size>;
+ static constexpr bool has_static_size = true;
+ static constexpr size_t static_size = VectorType::static_size;
+ static constexpr bool is_contiguous = true;
+
+ template <size_t SZ = static_size>
+ static inline VectorType create(const size_t sz, const K& val = suitable_default<K>::value())
+ {
+ static_assert(SZ == static_size, "Creation of Vector with different size not implemented!");
+ if (sz != SZ)
+ DUNE_THROW(Exceptions::wrong_input_given,
+ "You are trying to construct a FieldVector< ..., " << SZ << " > (of "
+ << "static size) with " << sz << " elements!");
+ return VectorType(val);
+ }
+};
+
+
+template <class V>
+typename std::enable_if<
+ is_vector<V>::value && VectorAbstraction<V>::has_static_size,
+ std::unique_ptr<FieldVector<typename VectorAbstraction<V>::S, VectorAbstraction<V>::static_size>>>::type
+make_field_container_ptr(const V& vec)
+{
+ auto ret = std::make_unique<FieldVector<typename VectorAbstraction<V>::S, VectorAbstraction<V>::static_size>>;
+ for (size_t ii = 0; ii < ret->size(); ++ii)
+ (*ret)[ii] = vec[ii];
+ return std::move(ret);
+}
+
+
+template <class V>
+typename std::enable_if<is_vector<V>::value && VectorAbstraction<V>::has_static_size,
+ FieldVector<typename VectorAbstraction<V>::S, VectorAbstraction<V>::static_size>>::type
+make_field_container(const V& vec)
+{
+ FieldVector<typename VectorAbstraction<V>::S, VectorAbstraction<V>::static_size> ret;
+ for (size_t ii = 0; ii < ret.size(); ++ii)
+ ret[ii] = vec[ii];
+ return ret;
+}
+
+
+template <class K, int SIZE>
+FieldVector<K, SIZE> make_field_container(Dune::FieldVector<K, SIZE>&& vec)
+{
+ return std::move(vec);
+}
+
+
+template <class K, int SIZE>
+FieldVector<K, SIZE> make_field_container(FieldVector<K, SIZE>&& vec)
+{
+ return std::move(vec);
+}
+
+template <class K>
+FieldVector<K, 3> cross_product(const FieldVector<K, 3>& u, const FieldVector<K, 3>& v)
+{
+ FieldVector<K, 3> ret;
+ ret[0] = u[1] * v[2] - u[2] * v[1];
+ ret[1] = u[2] * v[0] - u[0] * v[2];
+ ret[2] = u[0] * v[1] - u[1] * v[0];
+ return ret;
+}
+
+
+namespace internal {
+
+
+/**
+ * Most likely, you do not want to use this class directly, but \sa hstack instead.
+ */
+template <class K>
+struct hstack_decay
+{
+ static_assert(is_arithmetic<K>::value, "");
+ using type = K;
+};
+
+template <class K, int SIZE>
+struct hstack_decay<Dune::FieldVector<K, SIZE>>
+{
+ using type = Dune::XT::Common::FieldVector<K, SIZE>;
+};
+
+template <class K, int SIZE>
+struct hstack_decay<const Dune::FieldVector<K, SIZE>&> : public hstack_decay<Dune::FieldVector<K, SIZE>>
+{};
+
+template <class K, int SIZE>
+struct hstack_decay<Dune::XT::Common::FieldVector<K, SIZE>> : public hstack_decay<Dune::FieldVector<K, SIZE>>
+{};
+
+template <class K, int SIZE>
+struct hstack_decay<const Dune::XT::Common::FieldVector<K, SIZE>&> : public hstack_decay<Dune::FieldVector<K, SIZE>>
+{};
+
+
+/**
+ * Most likely, you do not want to use this, but \sa hstack instead.
+ */
+template <class T>
+using hstack_decay_t = typename hstack_decay<T>::type;
+
+
+/**
+ * Most likely, you do not want to use this class directly, but \sa hstack instead.
+ */
+template <class... Vectors>
+struct hstack_helper;
+
+template <class L, class R, class... Vectors>
+struct hstack_helper<L, R, Vectors...>
+{
+ using type =
+ typename hstack_helper<typename hstack_helper<hstack_decay_t<L>, hstack_decay_t<R>>::type, Vectors...>::type;
+};
+
+template <class K, int SIZE>
+struct hstack_helper<Dune::FieldVector<K, SIZE>>
+{
+ using type = hstack_decay_t<Dune::FieldVector<K, SIZE>>;
+};
+
+template <class K, int SIZE>
+struct hstack_helper<Dune::XT::Common::FieldVector<K, SIZE>> : public hstack_helper<Dune::FieldVector<K, SIZE>>
+{};
+
+template <class KL, class KR, int r>
+struct hstack_helper<KL, Dune::FieldVector<KR, r>>
+{
+ using type = hstack_decay_t<Dune::FieldVector<typename PromotionTraits<KL, KR>::PromotedType, 1 + r>>;
+};
+
+template <class KL, class KR, int r>
+struct hstack_helper<KL, Dune::XT::Common::FieldVector<KR, r>> : public hstack_helper<KL, Dune::FieldVector<KR, r>>
+{};
+
+template <class KL, int l, class KR>
+struct hstack_helper<Dune::FieldVector<KL, l>, KR>
+{
+ using type = hstack_decay_t<Dune::FieldVector<typename PromotionTraits<KL, KR>::PromotedType, l + 1>>;
+};
+
+template <class KL, int l, class KR>
+struct hstack_helper<Dune::XT::Common::FieldVector<KL, l>, KR> : public hstack_helper<Dune::FieldVector<KL, l>, KR>
+{};
+
+template <class KL, int l, class KR, int r>
+struct hstack_helper<Dune::FieldVector<KL, l>, Dune::FieldVector<KR, r>>
+{
+ using type = hstack_decay_t<Dune::FieldVector<typename PromotionTraits<KL, KR>::PromotedType, l + r>>;
+};
+
+template <class KL, int l, class KR, int r>
+struct hstack_helper<Dune::XT::Common::FieldVector<KL, l>, Dune::FieldVector<KR, r>>
+ : public hstack_helper<Dune::FieldVector<KL, l>, Dune::FieldVector<KR, r>>
+{};
+
+template <class KL, int l, class KR, int r>
+struct hstack_helper<Dune::FieldVector<KL, l>, Dune::XT::Common::FieldVector<KR, r>>
+ : public hstack_helper<Dune::FieldVector<KL, l>, Dune::FieldVector<KR, r>>
+{};
+
+template <class KL, int l, class KR, int r>
+struct hstack_helper<Dune::XT::Common::FieldVector<KL, l>, Dune::XT::Common::FieldVector<KR, r>>
+ : public hstack_helper<Dune::FieldVector<KL, l>, Dune::FieldVector<KR, r>>
+{};
+
+// not sure why this is required, would've hoped the decay above would take care of it
+template <class KL, int l, class KR, int r>
+struct hstack_helper<Dune::FieldVector<KL, l>, const Dune::FieldVector<KR, r>&>
+ : public hstack_helper<Dune::FieldVector<KL, l>, Dune::FieldVector<KR, r>>
+{};
+
+// not sure why this is required, would've hoped the decay above would take care of it
+template <class KL, int l, class KR, int r>
+struct hstack_helper<Dune::FieldVector<KL, l>, const Dune::XT::Common::FieldVector<KR, r>&>
+ : public hstack_helper<Dune::FieldVector<KL, l>, Dune::FieldVector<KR, r>>
+{};
+
+
+} // namespace internal
+
+
+/// \brief Specialization of \sa hstack to stop the recursion.
+template <class K, int SIZE>
+Dune::XT::Common::FieldVector<K, SIZE> hstack(Dune::FieldVector<K, SIZE>&& vec)
+{
+ return std::move(vec);
+}
+
+/// \brief Specialization of \sa hstack to stop the recursion.
+template <class K, int SIZE>
+Dune::XT::Common::FieldVector<K, SIZE> hstack(Dune::XT::Common::FieldVector<K, SIZE>&& vec)
+{
+ return std::move(vec);
+}
+
+/// \brief Specialization of \sa hstack to stop the recursion.
+template <class K, int SIZE>
+Dune::XT::Common::FieldVector<K, SIZE> hstack(const Dune::FieldVector<K, SIZE>& vec)
+{
+ return vec;
+}
+
+/// \brief Specialization of \sa hstack to stop the recursion.
+template <class K, int SIZE>
+Dune::XT::Common::FieldVector<K, SIZE> hstack(const Dune::XT::Common::FieldVector<K, SIZE>& vec)
+{
+ return vec;
+}
+
+
+/// \brief Specialization of \sa hstack for two arguments.
+template <class KL, class KR, int r>
+typename std::enable_if<
+ is_arithmetic<KL>::value && !is_field_vector<KL>::value,
+ typename internal::hstack_helper<internal::hstack_decay_t<KL>,
+ internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type>::type
+hstack(const KL& left, const Dune::FieldVector<KR, r>& right)
+{
+ typename internal::hstack_helper<internal::hstack_decay_t<KL>,
+ internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type ret;
+ ret[0] = left;
+ for (size_t ii = 0; ii < r; ++ii)
+ ret[ii + 1] = right[ii];
+ return ret;
+} // ... hstack<1, r>(...)
+
+/// \brief Specialization of \sa hstack for two arguments.
+template <class KL, class KR, int r>
+typename std::enable_if<
+ is_arithmetic<KL>::value && !is_field_vector<KL>::value,
+ typename internal::hstack_helper<internal::hstack_decay_t<KL>,
+ internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type>::type
+hstack(const KL& left, const Dune::XT::Common::FieldVector<KR, r>& right)
+{
+ return hstack(left, static_cast<const Dune::FieldVector<KR, r>&>(right));
+}
+
+/// \brief Specialization of \sa hstack for two arguments.
+template <class KL, int l, class KR>
+typename std::enable_if<is_arithmetic<KR>::value && !is_field_vector<KR>::value,
+ typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
+ internal::hstack_decay_t<KR>>::type>::type
+hstack(const Dune::FieldVector<KL, l>& left, const KR& right)
+{
+ typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
+ internal::hstack_decay_t<KR>>::type ret;
+ for (size_t ii = 0; ii < l; ++ii)
+ ret[ii] = left[ii];
+ ret[l] = right;
+ return ret;
+} // ... hstack<l, 1>(...)
+
+/// \brief Specialization of \sa hstack for two arguments.
+template <class KL, int l, class KR>
+typename std::enable_if<is_arithmetic<KR>::value && !is_field_vector<KR>::value,
+ typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
+ internal::hstack_decay_t<KR>>::type>::type
+hstack(const Dune::XT::Common::FieldVector<KL, l>& left, const KR& right)
+{
+ return hstack(static_cast<const Dune::FieldVector<KL, l>&>(left), right);
+}
+
+/// \brief Specialization of \sa hstack for two arguments.
+template <class KL, int l, class KR, int r>
+typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
+ internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type
+hstack(const Dune::FieldVector<KL, l>& left, const Dune::FieldVector<KR, r>& right)
+{
+ typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
+ internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type ret;
+ size_t ii = 0;
+ for (size_t jj = 0; jj < l; ++jj, ++ii)
+ ret[ii] = left[jj];
+ for (size_t jj = 0; jj < r; ++jj, ++ii)
+ ret[ii] = right[jj];
+ return ret;
+} // ... hstack<l, r>(...)
+
+/// \brief Specialization of \sa hstack for two arguments.
+template <class KL, int l, class KR, int r>
+typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
+ internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type
+hstack(const Dune::XT::Common::FieldVector<KL, l>& left, const Dune::FieldVector<KR, r>& right)
+{
+ return hstack(static_cast<const Dune::FieldVector<KL, l>&>(left), right);
+}
+
+/// \brief Specialization of \sa hstack for two arguments.
+template <class KL, int l, class KR, int r>
+typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
+ internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type
+hstack(const Dune::FieldVector<KL, l>& left, const Dune::XT::Common::FieldVector<KR, r>& right)
+{
+ return hstack(left, static_cast<const Dune::FieldVector<KR, r>&>(right));
+}
+
+/// \brief Specialization of \sa hstack for two arguments.
+template <class KL, int l, class KR, int r>
+typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
+ internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type
+hstack(const Dune::XT::Common::FieldVector<KL, l>& left, const Dune::XT::Common::FieldVector<KR, r>& right)
+{
+ return hstack(static_cast<const Dune::FieldVector<KL, l>&>(left),
+ static_cast<const Dune::FieldVector<KR, r>&>(right));
+}
+
+
+/**
+ * \brief Stacks an arbitrary number of numbers and FieldVectors horizontally into one appropriate FieldVector.
+ * \note At least one argument must be a FieldVector!
+ */
+template <class L, class R, class... Vectors>
+typename std::enable_if<
+ is_field_vector<L>::value || is_field_vector<R>::value,
+ typename internal::hstack_helper<internal::hstack_decay_t<L>, internal::hstack_decay_t<R>, Vectors...>::type>::type
+hstack(const L& left, const R& right, Vectors&&... vectors)
+{
+ return hstack(hstack(left, right), std::forward<Vectors>(vectors)...);
+}
+
+
+} // namespace Common
+} // namespace XT
+
+
+template <class K, int SIZE>
+typename std::enable_if<SIZE != 1, K>::type operator*(const Dune::FieldVector<K, SIZE>& vec,
+ const Dune::FieldMatrix<K, SIZE, 1>& mat)
+{
+ K ret = 0;
+ for (size_t ii = 0; ii < SIZE; ++ii)
+ ret += vec[ii] * mat[ii][0];
+ return ret;
+}
+
+template <class K, int SIZE>
+struct FieldTraits<XT::Common::FieldVector<K, SIZE>>
+{
+ typedef typename FieldTraits<K>::field_type field_type;
+ typedef typename FieldTraits<K>::real_type real_type;
+};
+
+
+} // namespace Dune
+namespace std {
+
+
+/**
+ * \sa http://en.cppreference.com/w/cpp/utility/hash
+ * \sa http://www.boost.org/doc/libs/1_55_0/doc/html/hash/combine.html
+ */
+template <class K, int SIZE>
+struct hash<Dune::FieldVector<K, SIZE>>
+{
+ typedef Dune::FieldVector<K, SIZE> argument_type;
+ typedef std::size_t result_type;
+
+ result_type operator()(argument_type const& s) const noexcept
+ {
+ std::size_t seed = 0;
+ for (size_t ii = 0; ii < s.size(); ++ii)
+ boost::hash_combine(seed, s[ii]);
+ return seed;
+ }
+}; // struct hash<Dune::FieldVector<...>>
+
+
+/**
+ * \sa http://en.cppreference.com/w/cpp/utility/hash
+ * \sa http://www.boost.org/doc/libs/1_55_0/doc/html/hash/combine.html
+ */
+template <class K, int SIZE>
+struct hash<Dune::XT::Common::FieldVector<K, SIZE>>
+{
+ typedef Dune::XT::Common::FieldVector<K, SIZE> argument_type;
+ typedef std::size_t result_type;
+
+ result_type operator()(argument_type const& s) const noexcept
+ {
+ std::size_t seed = 0;
+ for (size_t ii = 0; ii < s.size(); ++ii)
+ boost::hash_combine(seed, s[ii]);
+ return seed;
+ }
+}; // struct hash<Dune::XT::Common::FieldVector<...>>
+
+
+} // namespace std
+
+#endif // DUNE_XT_COMMON_FVECTOR_27_HH
diff --git a/dune/xt/common/fvector.hh b/dune/xt/common/fvector.hh
index a33cda87b482761aa0bb971595134d804e3ae70c..f4ab76dcf8fc148e699657d1ac075cc23be00de7 100644
--- a/dune/xt/common/fvector.hh
+++ b/dune/xt/common/fvector.hh
@@ -13,844 +13,12 @@
#ifndef DUNE_XT_COMMON_FVECTOR_HH
#define DUNE_XT_COMMON_FVECTOR_HH
-#include <functional>
-#include <initializer_list>
-#include <numeric>
-#include <type_traits>
-#include <vector>
+#include <dune/common/version.hh>
-#include <boost/functional/hash.hpp>
-
-#include <dune/common/fvector.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/common/promotiontraits.hh>
-
-#include <dune/xt/common/debug.hh>
-#include <dune/xt/common/densevector.hh>
-#include <dune/xt/common/exceptions.hh>
-#include <dune/xt/common/float_cmp.hh>
-#include <dune/xt/common/vector.hh>
-#include <dune/xt/common/type_traits.hh>
-
-namespace Dune {
-namespace XT {
-namespace Common {
-
-
-/**
- * \todo We need to implement all operators from the base which return the base, to rather return ourselfes!
- */
-template <class K, int SIZE>
-class FieldVector : public Dune::FieldVector<K, SIZE>
-{
- static_assert(SIZE >= 0, "Really?");
-
- typedef Dune::FieldVector<K, SIZE> BaseType;
- typedef FieldVector<K, SIZE> ThisType;
-
-public:
- FieldVector(const K& kk = suitable_default<K>::value())
- : BaseType(kk)
- {}
-
- FieldVector(const ThisType& other) = default;
-
- FieldVector(const BaseType& other)
- : BaseType(other)
- {}
-
- /* FieldMatrix< K, 1, 1 > is convertible to K, which in turn is convertible to FieldVector< K, 1 >. Without the
- * following constructor, this leads to an "ambiguous constructor" error (candidates are copy constructor and
- * constructor taking a K) */
- template <class Type = K>
- FieldVector(const typename std::enable_if<SIZE == 1 && std::is_same<K, Type>::value,
- typename Dune::FieldMatrix<K, 1, 1>>::type& mat)
- : BaseType(mat[0][0])
- {}
-
- FieldVector(const std::vector<K>& vec)
- : BaseType()
- {
-#ifndef NDEBUG
- if (vec.size() != SIZE)
- DUNE_THROW(Exceptions::wrong_input_given,
- "You are trying to construct a FieldVector< ..., " << SIZE << " > (of "
- << "static size) from a vector of size "
- << vec.size() << "!");
-#endif // NDEBUG
- for (size_t ii = 0; ii < SIZE; ++ii)
- this->operator[](ii) = vec[ii];
- } // FieldVector(...)
-
- FieldVector(std::initializer_list<K> list)
- : BaseType()
- {
-#ifndef NDEBUG
- if (list.size() != SIZE)
- DUNE_THROW(Exceptions::wrong_input_given,
- "You are trying to construct a FieldVector< ..., " << SIZE << " > (of "
- << "static size) from a list of size "
- << list.size() << "!");
-#endif // NDEBUG
- size_t ii = 0;
- for (auto element : list)
- this->operator[](ii++) = element;
- } // FieldVector(...)
-
- template <class C>
- FieldVector(
- const DenseVector<C>& x,
- typename std::enable_if<IsFieldVectorSizeCorrect<C, SIZE>::value
- && std::is_convertible<typename DenseVector<C>::value_type, K>::value>::type* /*dummy*/
- = nullptr)
- : BaseType(x)
- {}
-
- operator std::vector<K>() const
- {
- std::vector<K> ret(SIZE);
- for (size_t ii = 0; ii < SIZE; ++ii)
- ret[ii] = this->operator[](ii);
- return ret;
- }
-
- operator Dune::FieldMatrix<K, 1, SIZE>() const
- {
- Dune::FieldMatrix<K, 1, SIZE> ret;
- ret[0] = *this;
- return ret;
- }
-
- template <int S = SIZE>
- operator typename std::enable_if<(S == SIZE) && (SIZE != 1), Dune::FieldMatrix<K, S, 1>>::type() const
- {
- Dune::FieldMatrix<K, SIZE, 1> ret;
- for (size_t ii = 0; ii < SIZE; ++ii)
- ret[ii][0] = this->operator[](ii);
- return ret;
- }
-
- operator std::array<K, SIZE>() const
- {
- std::array<K, SIZE> ret;
- for (size_t ii = 0; ii < SIZE; ++ii)
- ret[ii] = this->operator[](ii);
- return ret;
- }
-
- template <class K_Other>
- typename multiplication_promotion<K, K_Other>::type operator*(const Dune::FieldVector<K_Other, SIZE>& other) const
- {
- typename multiplication_promotion<K, K_Other>::type result(0.);
- for (size_t ii = 0; ii < SIZE; ++ii)
- result += (*this)[ii] * other[ii];
- return result;
- }
-
- template <int R>
- typename std::enable_if<R == SIZE, K>::type operator*(const Dune::FieldMatrix<K, R, 1>& mat) const
- {
- K ret = 0;
- for (size_t ii = 0; ii < SIZE; ++ii)
- ret += (*this)[ii] * mat[ii][0];
- return ret;
- }
-
- //! This op is not redundant
- ThisType operator*(const K& scalar) const
- {
- ThisType ret(*this);
- ret *= scalar;
- return ret;
- }
-
- ThisType operator/(const K& scalar) const
- {
- ThisType ret(*this);
- ret /= scalar;
- return ret;
- }
-}; // class FieldVector
-
-
-template <class K, size_t block_num, size_t size_block>
-class BlockedFieldVector
-{
- using ThisType = BlockedFieldVector;
-
-public:
- static constexpr size_t num_blocks = block_num;
- static constexpr size_t block_size = size_block;
- static constexpr size_t static_size = num_blocks * block_size;
- using VectorType = Dune::FieldVector<K, static_size>;
- using XtVectorType = FieldVector<K, static_size>;
- using BlockType = FieldVector<K, block_size>;
-
- BlockedFieldVector(const K& val = K(0.))
- : backend_(BlockType(val))
- {}
-
- template <class OtherVectorType>
- BlockedFieldVector(const OtherVectorType& other,
- std::enable_if_t<is_vector<OtherVectorType>::value, int> /*dummy*/ = 0)
- {
- *this = other;
- }
-
- BlockedFieldVector(const BlockType& block)
- : backend_(block)
- {}
-
- BlockedFieldVector(const std::initializer_list<K>& init_list)
- {
- assert(init_list.size() == static_size || init_list.size() == num_blocks);
- if (init_list.size() == static_size) {
- *this = VectorType(init_list);
- } else {
- size_t jj = 0;
- for (const auto& block_init : init_list)
- block(jj++) = block_init;
- }
- }
-
- template <class OtherVectorType>
- std::enable_if_t<is_vector<OtherVectorType>::value, ThisType>& operator=(const OtherVectorType& other)
- {
- assert(other.size() == static_size);
- for (size_t jj = 0; jj < num_blocks; ++jj)
- for (size_t ii = 0; ii < block_size; ++ii)
- backend_[jj][ii] = other[jj * block_size + ii];
- return *this;
- }
-
- size_t size() const
- {
- return static_size;
- }
-
- K& operator[](const size_t ii)
- {
- assert(ii < static_size);
- return backend_[ii / block_size][ii % block_size];
- }
-
- const K& operator[](const size_t ii) const
- {
- assert(ii < static_size);
- return backend_[ii / block_size][ii % block_size];
- }
-
- K& get_entry(const size_t jj, const size_t ii)
- {
- assert(jj < num_blocks);
- assert(ii < block_size);
- return backend_[jj][ii];
- }
-
- const K& get_entry(const size_t jj, const size_t ii) const
- {
- assert(jj < num_blocks);
- assert(ii < block_size);
- return backend_[jj][ii];
- }
-
- BlockType& block(const size_t jj)
- {
- assert(jj < num_blocks);
- return backend_[jj];
- }
-
- const BlockType& block(const size_t jj) const
- {
- assert(jj < num_blocks);
- return backend_[jj];
- }
-
- K one_norm() const
- {
- K ret(0);
- for (const auto& block_vec : backend_)
- for (const auto& entry : block_vec)
- ret += std::abs(entry);
- return ret;
- }
-
- K two_norm() const
- {
- return std::sqrt(two_norm2());
- }
-
- K two_norm2() const
- {
- K ret(0);
- for (const auto& block_vec : backend_)
- for (const auto& entry : block_vec)
- ret += std::pow(entry, 2);
- return ret;
- }
-
- operator XtVectorType() const
- {
- XtVectorType ret(0.);
- for (size_t jj = 0; jj < num_blocks; ++jj)
- for (size_t ii = 0; ii < block_size; ++ii)
- ret[jj * block_size + ii] = backend_[jj][ii];
- return ret;
- }
-
- K* data()
- {
- return &(backend_[0][0]);
- }
-
- const K* data() const
- {
- return &(backend_[0][0]);
- }
-
- K* begin()
- {
- return data();
- }
-
- const K* begin() const
- {
- return data();
- }
-
- K* end()
- {
- return &(backend_[num_blocks - 1][block_size]);
- }
-
- const K* end() const
- {
- return &(backend_[num_blocks - 1][block_size]);
- }
-
- ThisType& operator*=(const K& val)
- {
- for (size_t jj = 0; jj < num_blocks; ++jj)
- backend_[jj] *= val;
- return *this;
- }
-
- ThisType operator*(const K& val) const
- {
- auto ret = *this;
- ret *= val;
- return ret;
- }
-
- K operator*(const ThisType& other) const
- {
- return std::inner_product(begin(), end(), other.begin(), 0.);
- }
-
- ThisType& operator+=(const ThisType& other)
- {
- backend_ += other.backend_;
- return *this;
- }
-
- ThisType operator+(const ThisType& other) const
- {
- auto ret = *this;
- ret += other;
- return ret;
- }
-
- ThisType& operator-=(const ThisType& other)
- {
- backend_ -= other.backend_;
- return *this;
- }
-
- ThisType operator-(const ThisType& other) const
- {
- auto ret = *this;
- ret -= other;
- return ret;
- }
-
-private:
- FieldVector<BlockType, num_blocks> backend_;
-};
-
-
-//! this allows to set the init value of the FieldVector at compile time
-template <class K, int SIZE, K value>
-class ValueInitFieldVector : public Dune::XT::Common::FieldVector<K, SIZE>
-{
- typedef Dune::XT::Common::FieldVector<K, SIZE> BaseType;
-
-public:
- ValueInitFieldVector()
- : BaseType(value)
- {}
-}; // class ValueInitFieldVector
-
-
-//! struct to be used as comparison function e.g. in a std::map<FieldVector<...>, ..., FieldVectorLess>
-struct FieldVectorLess
-{
- template <class FieldType, int dimDomain>
- bool operator()(const Dune::FieldVector<FieldType, dimDomain>& a,
- const Dune::FieldVector<FieldType, dimDomain>& b) const
- {
- for (size_t dd = 0; dd < dimDomain; ++dd) {
- if (a[dd] < b[dd])
- return true;
- else if (a[dd] > b[dd])
- return false;
- }
- return false;
- }
-};
-
-struct FieldVectorFloatLess
-{
- template <class FieldType, int dimDomain>
- bool operator()(const Dune::FieldVector<FieldType, dimDomain>& a,
- const Dune::FieldVector<FieldType, dimDomain>& b) const
- {
- for (size_t dd = 0; dd < dimDomain; ++dd) {
- if (XT::Common::FloatCmp::lt(a[dd], b[dd]))
- return true;
- else if (XT::Common::FloatCmp::gt(a[dd], b[dd]))
- return false;
- }
- return false;
- }
-};
-
-
-//! Specialization of VectorAbstraction for Dune::XT::Common::FieldVector
-template <class K, int SIZE>
-struct VectorAbstraction<Dune::XT::Common::FieldVector<K, SIZE>>
- : public internal::VectorAbstractionBase<Dune::XT::Common::FieldVector<K, SIZE>, K>
- , public internal::HasSubscriptOperatorForVectorAbstraction<Dune::XT::Common::FieldVector<K, SIZE>,
- typename Dune::FieldTraits<K>::field_type>
-{
- static const bool has_static_size = true;
- static const size_t static_size = SIZE;
- static const bool is_contiguous = true;
-
- template <size_t SZ = SIZE, class Field = K>
- using VectorTypeTemplate = Dune::XT::Common::FieldVector<Field, SZ>;
-
- template <size_t SZ = SIZE>
- static inline VectorTypeTemplate<SZ> create(const size_t sz, const K& val = suitable_default<K>::value())
- {
- if (sz != SZ)
- DUNE_THROW(Exceptions::wrong_input_given,
- "You are trying to construct a FieldVector< ..., " << SZ << " > (of "
- << "static size) with " << sz << " elements!");
- return VectorTypeTemplate<SZ>(val);
- }
-};
-
-
-//! Specialization of VectorAbstraction for Dune::XT::Common::BlockedFieldVector
-template <class K, size_t num_blocks, size_t block_size>
-struct VectorAbstraction<Dune::XT::Common::BlockedFieldVector<K, num_blocks, block_size>>
- : public internal::VectorAbstractionBase<Dune::XT::Common::BlockedFieldVector<K, num_blocks, block_size>, K>
- , public internal::HasSubscriptOperatorForVectorAbstraction<
- Dune::XT::Common::BlockedFieldVector<K, num_blocks, block_size>,
- typename Dune::FieldTraits<K>::field_type>
-{
- using VectorType = Dune::XT::Common::BlockedFieldVector<K, num_blocks, block_size>;
- static constexpr bool has_static_size = true;
- static constexpr size_t static_size = VectorType::static_size;
- static constexpr bool is_contiguous = true;
-
- template <size_t SZ = static_size>
- static inline VectorType create(const size_t sz, const K& val = suitable_default<K>::value())
- {
- static_assert(SZ == static_size, "Creation of Vector with different size not implemented!");
- if (sz != SZ)
- DUNE_THROW(Exceptions::wrong_input_given,
- "You are trying to construct a FieldVector< ..., " << SZ << " > (of "
- << "static size) with " << sz << " elements!");
- return VectorType(val);
- }
-};
-
-
-template <class V>
-typename std::enable_if<
- is_vector<V>::value && VectorAbstraction<V>::has_static_size,
- std::unique_ptr<FieldVector<typename VectorAbstraction<V>::S, VectorAbstraction<V>::static_size>>>::type
-make_field_container_ptr(const V& vec)
-{
- auto ret = std::make_unique<FieldVector<typename VectorAbstraction<V>::S, VectorAbstraction<V>::static_size>>;
- for (size_t ii = 0; ii < ret->size(); ++ii)
- (*ret)[ii] = vec[ii];
- return std::move(ret);
-}
-
-
-template <class V>
-typename std::enable_if<is_vector<V>::value && VectorAbstraction<V>::has_static_size,
- FieldVector<typename VectorAbstraction<V>::S, VectorAbstraction<V>::static_size>>::type
-make_field_container(const V& vec)
-{
- FieldVector<typename VectorAbstraction<V>::S, VectorAbstraction<V>::static_size> ret;
- for (size_t ii = 0; ii < ret.size(); ++ii)
- ret[ii] = vec[ii];
- return ret;
-}
-
-
-template <class K, int SIZE>
-FieldVector<K, SIZE> make_field_container(Dune::FieldVector<K, SIZE>&& vec)
-{
- return std::move(vec);
-}
-
-
-template <class K, int SIZE>
-FieldVector<K, SIZE> make_field_container(FieldVector<K, SIZE>&& vec)
-{
- return std::move(vec);
-}
-
-template <class K>
-FieldVector<K, 3> cross_product(const FieldVector<K, 3>& u, const FieldVector<K, 3>& v)
-{
- FieldVector<K, 3> ret;
- ret[0] = u[1] * v[2] - u[2] * v[1];
- ret[1] = u[2] * v[0] - u[0] * v[2];
- ret[2] = u[0] * v[1] - u[1] * v[0];
- return ret;
-}
-
-
-namespace internal {
-
-
-/**
- * Most likely, you do not want to use this class directly, but \sa hstack instead.
- */
-template <class K>
-struct hstack_decay
-{
- static_assert(is_arithmetic<K>::value, "");
- using type = K;
-};
-
-template <class K, int SIZE>
-struct hstack_decay<Dune::FieldVector<K, SIZE>>
-{
- using type = Dune::XT::Common::FieldVector<K, SIZE>;
-};
-
-template <class K, int SIZE>
-struct hstack_decay<const Dune::FieldVector<K, SIZE>&> : public hstack_decay<Dune::FieldVector<K, SIZE>>
-{};
-
-template <class K, int SIZE>
-struct hstack_decay<Dune::XT::Common::FieldVector<K, SIZE>> : public hstack_decay<Dune::FieldVector<K, SIZE>>
-{};
-
-template <class K, int SIZE>
-struct hstack_decay<const Dune::XT::Common::FieldVector<K, SIZE>&> : public hstack_decay<Dune::FieldVector<K, SIZE>>
-{};
-
-
-/**
- * Most likely, you do not want to use this, but \sa hstack instead.
- */
-template <class T>
-using hstack_decay_t = typename hstack_decay<T>::type;
-
-
-/**
- * Most likely, you do not want to use this class directly, but \sa hstack instead.
- */
-template <class... Vectors>
-struct hstack_helper;
-
-template <class L, class R, class... Vectors>
-struct hstack_helper<L, R, Vectors...>
-{
- using type =
- typename hstack_helper<typename hstack_helper<hstack_decay_t<L>, hstack_decay_t<R>>::type, Vectors...>::type;
-};
-
-template <class K, int SIZE>
-struct hstack_helper<Dune::FieldVector<K, SIZE>>
-{
- using type = hstack_decay_t<Dune::FieldVector<K, SIZE>>;
-};
-
-template <class K, int SIZE>
-struct hstack_helper<Dune::XT::Common::FieldVector<K, SIZE>> : public hstack_helper<Dune::FieldVector<K, SIZE>>
-{};
-
-template <class KL, class KR, int r>
-struct hstack_helper<KL, Dune::FieldVector<KR, r>>
-{
- using type = hstack_decay_t<Dune::FieldVector<typename PromotionTraits<KL, KR>::PromotedType, 1 + r>>;
-};
-
-template <class KL, class KR, int r>
-struct hstack_helper<KL, Dune::XT::Common::FieldVector<KR, r>> : public hstack_helper<KL, Dune::FieldVector<KR, r>>
-{};
-
-template <class KL, int l, class KR>
-struct hstack_helper<Dune::FieldVector<KL, l>, KR>
-{
- using type = hstack_decay_t<Dune::FieldVector<typename PromotionTraits<KL, KR>::PromotedType, l + 1>>;
-};
-
-template <class KL, int l, class KR>
-struct hstack_helper<Dune::XT::Common::FieldVector<KL, l>, KR> : public hstack_helper<Dune::FieldVector<KL, l>, KR>
-{};
-
-template <class KL, int l, class KR, int r>
-struct hstack_helper<Dune::FieldVector<KL, l>, Dune::FieldVector<KR, r>>
-{
- using type = hstack_decay_t<Dune::FieldVector<typename PromotionTraits<KL, KR>::PromotedType, l + r>>;
-};
-
-template <class KL, int l, class KR, int r>
-struct hstack_helper<Dune::XT::Common::FieldVector<KL, l>, Dune::FieldVector<KR, r>>
- : public hstack_helper<Dune::FieldVector<KL, l>, Dune::FieldVector<KR, r>>
-{};
-
-template <class KL, int l, class KR, int r>
-struct hstack_helper<Dune::FieldVector<KL, l>, Dune::XT::Common::FieldVector<KR, r>>
- : public hstack_helper<Dune::FieldVector<KL, l>, Dune::FieldVector<KR, r>>
-{};
-
-template <class KL, int l, class KR, int r>
-struct hstack_helper<Dune::XT::Common::FieldVector<KL, l>, Dune::XT::Common::FieldVector<KR, r>>
- : public hstack_helper<Dune::FieldVector<KL, l>, Dune::FieldVector<KR, r>>
-{};
-
-// not sure why this is required, would've hoped the decay above would take care of it
-template <class KL, int l, class KR, int r>
-struct hstack_helper<Dune::FieldVector<KL, l>, const Dune::FieldVector<KR, r>&>
- : public hstack_helper<Dune::FieldVector<KL, l>, Dune::FieldVector<KR, r>>
-{};
-
-// not sure why this is required, would've hoped the decay above would take care of it
-template <class KL, int l, class KR, int r>
-struct hstack_helper<Dune::FieldVector<KL, l>, const Dune::XT::Common::FieldVector<KR, r>&>
- : public hstack_helper<Dune::FieldVector<KL, l>, Dune::FieldVector<KR, r>>
-{};
-
-
-} // namespace internal
-
-
-/// \brief Specialization of \sa hstack to stop the recursion.
-template <class K, int SIZE>
-Dune::XT::Common::FieldVector<K, SIZE> hstack(Dune::FieldVector<K, SIZE>&& vec)
-{
- return std::move(vec);
-}
-
-/// \brief Specialization of \sa hstack to stop the recursion.
-template <class K, int SIZE>
-Dune::XT::Common::FieldVector<K, SIZE> hstack(Dune::XT::Common::FieldVector<K, SIZE>&& vec)
-{
- return std::move(vec);
-}
-
-/// \brief Specialization of \sa hstack to stop the recursion.
-template <class K, int SIZE>
-Dune::XT::Common::FieldVector<K, SIZE> hstack(const Dune::FieldVector<K, SIZE>& vec)
-{
- return vec;
-}
-
-/// \brief Specialization of \sa hstack to stop the recursion.
-template <class K, int SIZE>
-Dune::XT::Common::FieldVector<K, SIZE> hstack(const Dune::XT::Common::FieldVector<K, SIZE>& vec)
-{
- return vec;
-}
-
-
-/// \brief Specialization of \sa hstack for two arguments.
-template <class KL, class KR, int r>
-typename std::enable_if<
- is_arithmetic<KL>::value && !is_field_vector<KL>::value,
- typename internal::hstack_helper<internal::hstack_decay_t<KL>,
- internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type>::type
-hstack(const KL& left, const Dune::FieldVector<KR, r>& right)
-{
- typename internal::hstack_helper<internal::hstack_decay_t<KL>,
- internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type ret;
- ret[0] = left;
- for (size_t ii = 0; ii < r; ++ii)
- ret[ii + 1] = right[ii];
- return ret;
-} // ... hstack<1, r>(...)
-
-/// \brief Specialization of \sa hstack for two arguments.
-template <class KL, class KR, int r>
-typename std::enable_if<
- is_arithmetic<KL>::value && !is_field_vector<KL>::value,
- typename internal::hstack_helper<internal::hstack_decay_t<KL>,
- internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type>::type
-hstack(const KL& left, const Dune::XT::Common::FieldVector<KR, r>& right)
-{
- return hstack(left, static_cast<const Dune::FieldVector<KR, r>&>(right));
-}
-
-/// \brief Specialization of \sa hstack for two arguments.
-template <class KL, int l, class KR>
-typename std::enable_if<is_arithmetic<KR>::value && !is_field_vector<KR>::value,
- typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
- internal::hstack_decay_t<KR>>::type>::type
-hstack(const Dune::FieldVector<KL, l>& left, const KR& right)
-{
- typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
- internal::hstack_decay_t<KR>>::type ret;
- for (size_t ii = 0; ii < l; ++ii)
- ret[ii] = left[ii];
- ret[l] = right;
- return ret;
-} // ... hstack<l, 1>(...)
-
-/// \brief Specialization of \sa hstack for two arguments.
-template <class KL, int l, class KR>
-typename std::enable_if<is_arithmetic<KR>::value && !is_field_vector<KR>::value,
- typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
- internal::hstack_decay_t<KR>>::type>::type
-hstack(const Dune::XT::Common::FieldVector<KL, l>& left, const KR& right)
-{
- return hstack(static_cast<const Dune::FieldVector<KL, l>&>(left), right);
-}
-
-/// \brief Specialization of \sa hstack for two arguments.
-template <class KL, int l, class KR, int r>
-typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
- internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type
-hstack(const Dune::FieldVector<KL, l>& left, const Dune::FieldVector<KR, r>& right)
-{
- typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
- internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type ret;
- size_t ii = 0;
- for (size_t jj = 0; jj < l; ++jj, ++ii)
- ret[ii] = left[jj];
- for (size_t jj = 0; jj < r; ++jj, ++ii)
- ret[ii] = right[jj];
- return ret;
-} // ... hstack<l, r>(...)
-
-/// \brief Specialization of \sa hstack for two arguments.
-template <class KL, int l, class KR, int r>
-typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
- internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type
-hstack(const Dune::XT::Common::FieldVector<KL, l>& left, const Dune::FieldVector<KR, r>& right)
-{
- return hstack(static_cast<const Dune::FieldVector<KL, l>&>(left), right);
-}
-
-/// \brief Specialization of \sa hstack for two arguments.
-template <class KL, int l, class KR, int r>
-typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
- internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type
-hstack(const Dune::FieldVector<KL, l>& left, const Dune::XT::Common::FieldVector<KR, r>& right)
-{
- return hstack(left, static_cast<const Dune::FieldVector<KR, r>&>(right));
-}
-
-/// \brief Specialization of \sa hstack for two arguments.
-template <class KL, int l, class KR, int r>
-typename internal::hstack_helper<internal::hstack_decay_t<Dune::FieldVector<KL, l>>,
- internal::hstack_decay_t<Dune::FieldVector<KR, r>>>::type
-hstack(const Dune::XT::Common::FieldVector<KL, l>& left, const Dune::XT::Common::FieldVector<KR, r>& right)
-{
- return hstack(static_cast<const Dune::FieldVector<KL, l>&>(left),
- static_cast<const Dune::FieldVector<KR, r>&>(right));
-}
-
-
-/**
- * \brief Stacks an arbitrary number of numbers and FieldVectors horizontally into one appropriate FieldVector.
- * \note At least one argument must be a FieldVector!
- */
-template <class L, class R, class... Vectors>
-typename std::enable_if<
- is_field_vector<L>::value || is_field_vector<R>::value,
- typename internal::hstack_helper<internal::hstack_decay_t<L>, internal::hstack_decay_t<R>, Vectors...>::type>::type
-hstack(const L& left, const R& right, Vectors&&... vectors)
-{
- return hstack(hstack(left, right), std::forward<Vectors>(vectors)...);
-}
-
-
-} // namespace Common
-} // namespace XT
-
-
-template <class K, int SIZE>
-typename std::enable_if<SIZE != 1, K>::type operator*(const Dune::FieldVector<K, SIZE>& vec,
- const Dune::FieldMatrix<K, SIZE, 1>& mat)
-{
- K ret = 0;
- for (size_t ii = 0; ii < SIZE; ++ii)
- ret += vec[ii] * mat[ii][0];
- return ret;
-}
-
-template <class K, int SIZE>
-struct FieldTraits<XT::Common::FieldVector<K, SIZE>>
-{
- typedef typename FieldTraits<K>::field_type field_type;
- typedef typename FieldTraits<K>::real_type real_type;
-};
-
-
-} // namespace Dune
-namespace std {
-
-
-/**
- * \sa http://en.cppreference.com/w/cpp/utility/hash
- * \sa http://www.boost.org/doc/libs/1_55_0/doc/html/hash/combine.html
- */
-template <class K, int SIZE>
-struct hash<Dune::FieldVector<K, SIZE>>
-{
- typedef Dune::FieldVector<K, SIZE> argument_type;
- typedef std::size_t result_type;
-
- result_type operator()(argument_type const& s) const noexcept
- {
- std::size_t seed = 0;
- for (size_t ii = 0; ii < s.size(); ++ii)
- boost::hash_combine(seed, s[ii]);
- return seed;
- }
-}; // struct hash<Dune::FieldVector<...>>
-
-
-/**
- * \sa http://en.cppreference.com/w/cpp/utility/hash
- * \sa http://www.boost.org/doc/libs/1_55_0/doc/html/hash/combine.html
- */
-template <class K, int SIZE>
-struct hash<Dune::XT::Common::FieldVector<K, SIZE>>
-{
- typedef Dune::XT::Common::FieldVector<K, SIZE> argument_type;
- typedef std::size_t result_type;
-
- result_type operator()(argument_type const& s) const noexcept
- {
- std::size_t seed = 0;
- for (size_t ii = 0; ii < s.size(); ++ii)
- boost::hash_combine(seed, s[ii]);
- return seed;
- }
-}; // struct hash<Dune::XT::Common::FieldVector<...>>
-
-
-} // namespace std
+#if DUNE_VERSION_GTE(DUNE_COMMON, 2, 7)
+# include "fvector-2.7.hh"
+#else
+# include "fvector-2.6.hh"
+#endif
#endif // DUNE_XT_COMMON_FVECTOR_HH