From dc9153fb2ebb179731bd37051c80e77c78dc1e4c Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Mon, 25 Jul 2016 11:47:01 +0200
Subject: [PATCH] Incorporate addProduct from Arithmetic
---
dune/matrix-vector/CMakeLists.txt | 3 +
dune/matrix-vector/axpy.hh | 462 +++++++++++++++++++++++++++++
dune/matrix-vector/matrixtraits.hh | 55 ++++
dune/matrix-vector/scalartraits.hh | 44 +++
4 files changed, 564 insertions(+)
create mode 100644 dune/matrix-vector/axpy.hh
create mode 100644 dune/matrix-vector/matrixtraits.hh
create mode 100644 dune/matrix-vector/scalartraits.hh
diff --git a/dune/matrix-vector/CMakeLists.txt b/dune/matrix-vector/CMakeLists.txt
index 7a7517e..3ffe97f 100644
--- a/dune/matrix-vector/CMakeLists.txt
+++ b/dune/matrix-vector/CMakeLists.txt
@@ -1,7 +1,10 @@
#install headers
install(FILES
+ axpy.hh
componentwisematrixmap.hh
genericvectortools.hh
+ matrixtraits.hh
+ scalartraits.hh
singlenonzerocolumnmatrix.hh
singlenonzerorowmatrix.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/matrix-vector)
\ No newline at end of file
diff --git a/dune/matrix-vector/axpy.hh b/dune/matrix-vector/axpy.hh
new file mode 100644
index 0000000..792dc03
--- /dev/null
+++ b/dune/matrix-vector/axpy.hh
@@ -0,0 +1,462 @@
+#ifndef DUNE_MATRIX_VECTOR_AXPY_HH
+#define DUNE_MATRIX_VECTOR_AXPY_HH
+
+#include <type_traits>
+
+#include <dune/common/diagonalmatrix.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/scaledidmatrix.hh>
+
+#include "matrixtraits.hh"
+#include "scalartraits.hh"
+
+namespace Dune {
+namespace MatrixVector {
+ /** \brief Internal helper class for product operations
+ *
+ */
+ template <class A, class B, class C, bool AisScalar, bool BisScalar,
+ bool CisScalar>
+ struct ProductHelper {
+ template <
+ class ADummy = A,
+ std::enable_if_t<!MatrixTraits<ADummy>::isMatrix, int> SFINAE_Dummy = 0>
+ static void addProduct(A& a, const B& b, const C& c) {
+ b.umv(c, a);
+ }
+
+ template <
+ class ADummy = A,
+ std::enable_if_t<MatrixTraits<ADummy>::isMatrix, int> SFINAE_Dummy = 0>
+ static void addProduct(A& a, const B& b, const C& c) {
+ typename B::ConstRowIterator bi = b.begin();
+ typename B::ConstRowIterator bEnd = b.end();
+ for (; bi != bEnd; ++bi) {
+ typename B::ConstColIterator bik = bi->begin();
+ typename B::ConstColIterator biEnd = bi->end();
+ for (; bik != biEnd; ++bik) {
+ typename C::ConstColIterator ckj = c[bik.index()].begin();
+ typename C::ConstColIterator ckEnd = c[bik.index()].end();
+ for (; ckj != ckEnd; ++ckj)
+ a[bi.index()][ckj.index()] += (*bik) * (*ckj);
+ }
+ }
+ }
+ };
+
+ // Internal helper class for scaled product operations (i.e., b is always a
+ // scalar)
+ template <class A, class Scalar, class B, class C, bool AisScalar,
+ bool BisScalar, bool CisScalar>
+ struct ScaledProductHelper {
+ template <
+ class ADummy = A,
+ std::enable_if_t<!MatrixTraits<ADummy>::isMatrix, int> SFINAE_Dummy = 0>
+ static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c) {
+ b.usmv(scalar, c, a);
+ }
+
+ template <
+ class ADummy = A,
+ std::enable_if_t<MatrixTraits<ADummy>::isMatrix, int> SFINAE_Dummy = 0>
+ static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c) {
+ typename B::ConstRowIterator bi = b.begin();
+ typename B::ConstRowIterator bEnd = b.end();
+ for (; bi != bEnd; ++bi) {
+ typename B::ConstColIterator bik = b[bi.index()].begin();
+ typename B::ConstColIterator biEnd = b[bi.index()].end();
+ for (; bik != biEnd; ++bik) {
+ typename C::ConstColIterator ckj = c[bik.index()].begin();
+ typename C::ConstColIterator ckEnd = c[bik.index()].end();
+ for (; ckj != ckEnd; ++ckj)
+ a[bi.index()][ckj.index()] += scalar * (*bik) * (*ckj);
+ }
+ }
+ }
+ };
+
+ template <class T, int n, int m, int k>
+ struct ProductHelper<Dune::FieldMatrix<T, n, k>, Dune::FieldMatrix<T, n, m>,
+ Dune::FieldMatrix<T, m, k>, false, false, false> {
+ typedef Dune::FieldMatrix<T, n, k> A;
+ typedef Dune::FieldMatrix<T, n, m> B;
+ typedef Dune::FieldMatrix<T, m, k> C;
+ static void addProduct(A& a, const B& b, const C& c) {
+ for (size_t row = 0; row < n; ++row) {
+ for (size_t col = 0; col < k; ++col) {
+ for (size_t i = 0; i < m; ++i)
+ a[row][col] += b[row][i] * c[i][col];
+ }
+ }
+ }
+ };
+
+ template <class T, int n, int m, int k>
+ struct ScaledProductHelper<Dune::FieldMatrix<T, n, k>, T,
+ Dune::FieldMatrix<T, n, m>,
+ Dune::FieldMatrix<T, m, k>, false, false, false> {
+ typedef Dune::FieldMatrix<T, n, k> A;
+ typedef Dune::FieldMatrix<T, n, m> C;
+ typedef Dune::FieldMatrix<T, m, k> D;
+ static void addProduct(A& a, const T& b, const C& c, const D& d) {
+ for (size_t row = 0; row < n; ++row) {
+ for (size_t col = 0; col < k; ++col) {
+ for (size_t i = 0; i < m; ++i)
+ a[row][col] += b * c[row][i] * d[i][col];
+ }
+ }
+ }
+ };
+
+ template <class T, int n>
+ struct ProductHelper<Dune::DiagonalMatrix<T, n>, Dune::DiagonalMatrix<T, n>,
+ Dune::DiagonalMatrix<T, n>, false, false, false> {
+ typedef Dune::DiagonalMatrix<T, n> A;
+ typedef A B;
+ typedef B C;
+ static void addProduct(A& a, const B& b, const C& c) {
+ for (size_t i = 0; i < n; ++i)
+ a.diagonal(i) += b.diagonal(i) * c.diagonal(i);
+ }
+ };
+
+ template <class T, int n>
+ struct ScaledProductHelper<Dune::DiagonalMatrix<T, n>, T,
+ Dune::DiagonalMatrix<T, n>,
+ Dune::DiagonalMatrix<T, n>, false, false, false> {
+ typedef Dune::DiagonalMatrix<T, n> A;
+ typedef A C;
+ typedef C D;
+ static void addProduct(A& a, const T& b, const C& c, const D& d) {
+ for (size_t i = 0; i < n; ++i)
+ a.diagonal(i) += b * c.diagonal(i) * d.diagonal(i);
+ }
+ };
+
+ template <class T, int n>
+ struct ProductHelper<Dune::FieldMatrix<T, n, n>, Dune::DiagonalMatrix<T, n>,
+ Dune::DiagonalMatrix<T, n>, false, false, false> {
+ typedef Dune::FieldMatrix<T, n, n> A;
+ typedef Dune::DiagonalMatrix<T, n> B;
+ typedef B C;
+ static void addProduct(A& a, const B& b, const C& c) {
+ for (size_t i = 0; i < n; ++i)
+ a[i][i] += b.diagonal(i) * c.diagonal(i);
+ }
+ };
+
+ template <class T, int n>
+ struct ScaledProductHelper<Dune::FieldMatrix<T, n, n>, T,
+ Dune::DiagonalMatrix<T, n>,
+ Dune::DiagonalMatrix<T, n>, false, false, false> {
+ typedef Dune::FieldMatrix<T, n, n> A;
+ typedef Dune::DiagonalMatrix<T, n> C;
+ typedef C D;
+ static void addProduct(A& a, const T& b, const C& c, const D& d) {
+ for (size_t i = 0; i < n; ++i)
+ a[i][i] += b * c.diagonal(i) * d.diagonal(i);
+ }
+ };
+
+ template <class T, int n, int m>
+ struct ProductHelper<Dune::FieldMatrix<T, n, m>, Dune::DiagonalMatrix<T, n>,
+ Dune::FieldMatrix<T, n, m>, false, false, false> {
+ typedef Dune::FieldMatrix<T, n, m> A;
+ typedef Dune::DiagonalMatrix<T, n> B;
+ typedef A C;
+ static void addProduct(A& a, const B& b, const C& c) {
+ for (size_t i = 0; i < n; ++i)
+ a[i].axpy(b.diagonal(i), c[i]);
+ }
+ };
+
+ template <class T, int n, int m>
+ struct ScaledProductHelper<Dune::FieldMatrix<T, n, m>, T,
+ Dune::DiagonalMatrix<T, n>,
+ Dune::FieldMatrix<T, n, m>, false, false, false> {
+ typedef Dune::FieldMatrix<T, n, m> A;
+ typedef Dune::DiagonalMatrix<T, n> C;
+ typedef A D;
+ static void addProduct(A& a, const T& b, const C& c, const D& d) {
+ for (size_t i = 0; i < n; ++i)
+ a[i].axpy(b * c.diagonal(i), d[i]);
+ }
+ };
+
+ template <class T, int n, int m>
+ struct ProductHelper<Dune::FieldMatrix<T, n, m>, Dune::FieldMatrix<T, n, m>,
+ Dune::DiagonalMatrix<T, m>, false, false, false> {
+ typedef Dune::FieldMatrix<T, n, m> A;
+ typedef Dune::DiagonalMatrix<T, m> C;
+ typedef A B;
+
+ static void addProduct(A& a, const B& b, const C& c) {
+ for (size_t i = 0; i < n; ++i)
+ for (size_t j = 0; j < m; j++)
+ a[i][j] += c.diagonal(j) * b[i][j];
+ }
+ };
+
+ template <class T, int n, int m>
+ struct ScaledProductHelper<Dune::FieldMatrix<T, n, m>, T,
+ Dune::FieldMatrix<T, n, m>,
+ Dune::DiagonalMatrix<T, m>, false, false, false> {
+ typedef Dune::FieldMatrix<T, n, m> A;
+ typedef Dune::DiagonalMatrix<T, m> D;
+ typedef A C;
+
+ static void addProduct(A& a, const T& b, const C& c, const D& d) {
+ for (size_t i = 0; i < n; ++i)
+ for (size_t j = 0; j < m; j++)
+ a[i][j] += b * d.diagonal(j) * c[i][j];
+ }
+ };
+
+ template <class T, int n>
+ struct ProductHelper<Dune::ScaledIdentityMatrix<T, n>,
+ Dune::ScaledIdentityMatrix<T, n>,
+ Dune::ScaledIdentityMatrix<T, n>, false, false, false> {
+ typedef Dune::ScaledIdentityMatrix<T, n> A;
+ typedef A B;
+ typedef B C;
+
+ static void addProduct(A& a, const B& b, const C& c) {
+ a.scalar() += b.scalar() * c.scalar();
+ }
+ };
+
+ template <class T, class Scalar, int n>
+ struct ScaledProductHelper<Dune::ScaledIdentityMatrix<T, n>, Scalar,
+ Dune::ScaledIdentityMatrix<T, n>,
+ Dune::ScaledIdentityMatrix<T, n>, false, false,
+ false> {
+ typedef Dune::ScaledIdentityMatrix<T, n> A;
+ typedef A C;
+ typedef C D;
+
+ static void addProduct(A& a, const Scalar& b, const C& c, const D& d) {
+ a.scalar() += b * c.scalar() * d.scalar();
+ }
+ };
+
+ /** \brief Specialization for b being a scalar type
+ */
+ template <class A, class ScalarB, class C, bool AisScalar, bool CisScalar>
+ struct ProductHelper<A, ScalarB, C, AisScalar, true, CisScalar> {
+ typedef ScalarB B;
+
+ template <
+ class ADummy = A,
+ std::enable_if_t<!MatrixTraits<ADummy>::isMatrix, int> SFINAE_Dummy = 0>
+ static void addProduct(A& a, const B& b, const C& c) {
+ a.axpy(b, c);
+ }
+
+ template <
+ class ADummy = A,
+ std::enable_if_t<MatrixTraits<ADummy>::isMatrix, int> SFINAE_Dummy = 0>
+ static void addProduct(A& a, const B& b, const C& c) {
+ typename C::ConstRowIterator ci = c.begin();
+ typename C::ConstRowIterator cEnd = c.end();
+ for (; ci != cEnd; ++ci) {
+ typename C::ConstColIterator cik = c[ci.index()].begin();
+ typename C::ConstColIterator ciEnd = c[ci.index()].end();
+ for (; cik != ciEnd; ++cik) {
+ a[ci.index()][cik.index()] += b * (*cik);
+ }
+ }
+ }
+ };
+
+ template <class A, class Scalar, class ScalarB, class C, bool AisScalar,
+ bool CisScalar>
+ struct ScaledProductHelper<A, Scalar, ScalarB, C, AisScalar, true,
+ CisScalar> {
+ typedef ScalarB B;
+
+ template <
+ class ADummy = A,
+ std::enable_if_t<!MatrixTraits<ADummy>::isMatrix, int> SFINAE_Dummy = 0>
+ static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c) {
+ a.axpy(scalar * b, c);
+ }
+
+ template <
+ class ADummy = A,
+ std::enable_if_t<MatrixTraits<ADummy>::isMatrix, int> SFINAE_Dummy = 0>
+ static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c) {
+ typename C::ConstRowIterator ci = c.begin();
+ typename C::ConstRowIterator cEnd = c.end();
+ for (; ci != cEnd; ++ci) {
+ typename C::ConstColIterator cik = c[ci.index()].begin();
+ typename C::ConstColIterator ciEnd = c[ci.index()].end();
+ for (; cik != ciEnd; ++cik) {
+ a[ci.index()][cik.index()] += scalar * b * (*cik);
+ }
+ }
+ }
+ };
+
+ template <class ABlock, class ScalarB, class CBlock>
+ struct ProductHelper<Dune::BCRSMatrix<ABlock>, ScalarB,
+ Dune::BCRSMatrix<CBlock>, false, true, false> {
+ typedef Dune::BCRSMatrix<ABlock> A;
+ typedef ScalarB B;
+ typedef Dune::BCRSMatrix<CBlock> C;
+
+ static void addProduct(A& a, const B& b, const C& c) {
+ // Workaround for the fact that Dune::BCRSMatrix assumes
+ // its blocks to have an axpy() implementation, yet
+ // Dune::DiagonalMatrix does not.
+
+ // a.axpy(b, c);
+
+ for (typename C::ConstIterator i = c.begin(); i != c.end(); ++i) {
+ const size_t iindex = i.index();
+ for (typename C::row_type::ConstIterator j = i->begin(); j != i->end();
+ ++j)
+ ProductHelper<typename A::block_type, B, typename C::block_type,
+ false, true, false>::addProduct(a[iindex][j.index()], b,
+ *j);
+ }
+ }
+ };
+
+ template <class ABlock, class Scalar, class ScalarB, class CBlock>
+ struct ScaledProductHelper<Dune::BCRSMatrix<ABlock>, Scalar, ScalarB,
+ Dune::BCRSMatrix<CBlock>, false, true, false> {
+ typedef Dune::BCRSMatrix<ABlock> A;
+ typedef ScalarB B;
+ typedef Dune::BCRSMatrix<CBlock> C;
+
+ static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c) {
+ // Workaround for the fact that Dune::BCRSMatrix assumes
+ // its blocks to have an axpy() implementation, yet
+ // Dune::DiagonalMatrix does not.
+
+ // a.axpy(scalar*b, c);
+
+ for (typename C::ConstIterator i = c.begin(); i != c.end(); ++i) {
+ const size_t iindex = i.index();
+ for (typename C::row_type::ConstIterator j = i->begin(); j != i->end();
+ ++j)
+ ProductHelper<typename A::block_type, B, typename C::block_type,
+ false, true, false>::addProduct(a[iindex][j.index()],
+ scalar * b, *j);
+ }
+ }
+ };
+
+ template <class T, int n, class ScalarB>
+ struct ProductHelper<Dune::ScaledIdentityMatrix<T, n>, ScalarB,
+ Dune::ScaledIdentityMatrix<T, n>, false, true, false> {
+ typedef Dune::ScaledIdentityMatrix<T, n> A;
+ typedef ScalarB B;
+ typedef Dune::ScaledIdentityMatrix<T, n> C;
+
+ static void addProduct(A& a, const B& b, const C& c) {
+ a.scalar() += b * c.scalar();
+ }
+ };
+
+ template <class T, int n, class Scalar, class ScalarB>
+ struct ScaledProductHelper<Dune::ScaledIdentityMatrix<T, n>, Scalar, ScalarB,
+ Dune::ScaledIdentityMatrix<T, n>, false, true,
+ false> {
+ typedef Dune::ScaledIdentityMatrix<T, n> A;
+ typedef ScalarB B;
+ typedef Dune::ScaledIdentityMatrix<T, n> C;
+
+ static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c) {
+ a.scalar() += scalar * b * c.scalar();
+ }
+ };
+
+ template <class A, class B, class C>
+ struct ProductHelper<A, B, C, true, true, true> {
+ static void addProduct(A& a, const B& b, const C& c) { a += b * c; }
+ };
+
+ template <class A, class B, class C, class D>
+ struct ScaledProductHelper<A, B, C, D, true, true, true> {
+ static void addProduct(A& a, const B& b, const C& c, const D& d) {
+ a += b * c * d;
+ }
+ };
+
+ /** \brief Add a product to some matrix or vector
+ *
+ * This function computes a+=b*c.
+ *
+ * This function should tolerate all meaningful
+ * combinations of scalars, vectors, and matrices.
+ *
+ * a,b,c could be matrices with appropriate
+ * dimensions. But b can also always be a scalar
+ * represented by a 1-dim vector or a 1 by 1 matrix.
+ */
+ template <class A, class B, class C>
+ void addProduct(A& a, const B& b, const C& c) {
+ ProductHelper<A, B, C, ScalarTraits<A>::isScalar, ScalarTraits<B>::isScalar,
+ ScalarTraits<C>::isScalar>::addProduct(a, b, c);
+ }
+
+ /** \brief Subtract a product from some matrix or vector
+ *
+ * This function computes a-=b*c.
+ *
+ * This function should tolerate all meaningful
+ * combinations of scalars, vectors, and matrices.
+ *
+ * a,b,c could be matrices with appropriate
+ * dimensions. But b can also always be a scalar
+ * represented by a 1-dim vector or a 1 by 1 matrix.
+ */
+ template <class A, class B, class C>
+ void subtractProduct(A& a, const B& b, const C& c) {
+ ScaledProductHelper<A, int, B, C, ScalarTraits<A>::isScalar,
+ ScalarTraits<B>::isScalar,
+ ScalarTraits<C>::isScalar>::addProduct(a, -1, b, c);
+ }
+
+ /** \brief Add a scaled product to some matrix or vector
+ *
+ * This function computes a+=b*c*d.
+ *
+ * This function should tolerate all meaningful
+ * combinations of scalars, vectors, and matrices.
+ *
+ * a,c,d could be matrices with appropriate dimensions. But b must
+ * (currently) and c can also always be a scalar represented by a
+ * 1-dim vector or a 1 by 1 matrix.
+ */
+ template <class A, class B, class C, class D>
+ typename std::enable_if<ScalarTraits<B>::isScalar, void>::type addProduct(
+ A& a, const B& b, const C& c, const D& d) {
+ ScaledProductHelper<A, B, C, D, ScalarTraits<A>::isScalar,
+ ScalarTraits<C>::isScalar,
+ ScalarTraits<D>::isScalar>::addProduct(a, b, c, d);
+ }
+
+ /** \brief Subtract a scaled product from some matrix or vector
+ *
+ * This function computes a-=b*c*d.
+ *
+ * This function should tolerate all meaningful
+ * combinations of scalars, vectors, and matrices.
+ *
+ * a,c,d could be matrices with appropriate dimensions. But b must
+ * (currently) and c can also always be a scalar represented by a
+ * 1-dim vector or a 1 by 1 matrix.
+ */
+ template <class A, class B, class C, class D>
+ typename std::enable_if<ScalarTraits<B>::isScalar, void>::type
+ subtractProduct(A& a, const B& b, const C& c, const D& d) {
+ ScaledProductHelper<A, B, C, D, ScalarTraits<A>::isScalar,
+ ScalarTraits<C>::isScalar,
+ ScalarTraits<D>::isScalar>::addProduct(a, -b, c, d);
+ }
+}
+}
+#endif
diff --git a/dune/matrix-vector/matrixtraits.hh b/dune/matrix-vector/matrixtraits.hh
new file mode 100644
index 0000000..da72f65
--- /dev/null
+++ b/dune/matrix-vector/matrixtraits.hh
@@ -0,0 +1,55 @@
+#ifndef DUNE_MATRIX_VECTOR_MATRIXTRAITS_HH
+#define DUNE_MATRIX_VECTOR_MATRIXTRAITS_HH
+
+#include <dune/common/diagonalmatrix.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/scaledidmatrix.hh>
+
+namespace Dune {
+ namespace MatrixVector {
+
+ /** \brief Class to identify matrix types and extract information
+ *
+ * Specialize this class for all types that can be used like a matrix.
+ */
+ template<class T>
+ struct MatrixTraits
+ {
+ enum { isMatrix=false };
+ enum { rows=-1};
+ enum { cols=-1};
+ };
+
+ template<class T, int n, int m>
+ struct MatrixTraits<Dune::FieldMatrix<T,n,m> >
+ {
+ enum { isMatrix=true };
+ enum { rows=n};
+ enum { cols=m};
+ };
+
+ template<class T, int n>
+ struct MatrixTraits<Dune::DiagonalMatrix<T,n> >
+ {
+ enum { isMatrix=true };
+ enum { rows=n};
+ enum { cols=n};
+ };
+
+ template<class T, int n>
+ struct MatrixTraits<Dune::ScaledIdentityMatrix<T,n> >
+ {
+ enum { isMatrix=true };
+ enum { rows=n};
+ enum { cols=n};
+ };
+
+ template<class T>
+ struct MatrixTraits<Dune::BCRSMatrix<T> >
+ {
+ enum { isMatrix=true };
+ };
+ }
+}
+#endif
diff --git a/dune/matrix-vector/scalartraits.hh b/dune/matrix-vector/scalartraits.hh
new file mode 100644
index 0000000..5fabfc1
--- /dev/null
+++ b/dune/matrix-vector/scalartraits.hh
@@ -0,0 +1,44 @@
+#ifndef DUNE_MATRIX_VECTOR_SCALARTRAITS_HH
+#define DUNE_MATRIX_VECTOR_SCALARTRAITS_HH
+
+#include <dune/common/diagonalmatrix.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/scaledidmatrix.hh>
+
+namespace Dune {
+namespace MatrixVector {
+ /** \brief Class to identify scalar types
+ *
+ * Specialize this class for all types that can be used
+ * like scalar quantities.
+ */
+ template <class T>
+ struct ScalarTraits {
+ enum {
+ isScalar = (std::is_scalar<T>::value and not std::is_pointer<T>::value)
+ };
+ };
+
+ template <class T>
+ struct ScalarTraits<Dune::FieldVector<T, 1>> {
+ enum { isScalar = true };
+ };
+
+ template <class T>
+ struct ScalarTraits<Dune::FieldMatrix<T, 1, 1>> {
+ enum { isScalar = true };
+ };
+
+ template <class T>
+ struct ScalarTraits<Dune::DiagonalMatrix<T, 1>> {
+ enum { isScalar = true };
+ };
+
+ template <class T>
+ struct ScalarTraits<Dune::ScaledIdentityMatrix<T, 1>> {
+ enum { isScalar = true };
+ };
+}
+}
+#endif
--
GitLab