From c0d7b6e03fc2abf420c8e2c4bce9ad3ddc1db0fc Mon Sep 17 00:00:00 2001 From: Max Kahnt <max.kahnt@fu-berlin.de> Date: Tue, 3 Oct 2017 15:53:07 +0200 Subject: [PATCH] Add getScalar to extract the scale factor of some scaling. --- dune/matrix-vector/CMakeLists.txt | 3 +- dune/matrix-vector/getscalar.hh | 83 +++++++++++++++++++++++++++++++ 2 files changed, 85 insertions(+), 1 deletion(-) create mode 100644 dune/matrix-vector/getscalar.hh diff --git a/dune/matrix-vector/CMakeLists.txt b/dune/matrix-vector/CMakeLists.txt index 9e4da1f..c1f8cb1 100644 --- a/dune/matrix-vector/CMakeLists.txt +++ b/dune/matrix-vector/CMakeLists.txt @@ -1,4 +1,4 @@ -add_subdirectory(test) +add_subdirectory(test) add_subdirectory(traits) #install headers @@ -12,6 +12,7 @@ install(FILES concepts.hh crossproduct.hh genericvectortools.hh + getscalar.hh ldlt.hh promote.hh singlenonzerocolumnmatrix.hh diff --git a/dune/matrix-vector/getscalar.hh b/dune/matrix-vector/getscalar.hh new file mode 100644 index 0000000..608042c --- /dev/null +++ b/dune/matrix-vector/getscalar.hh @@ -0,0 +1,83 @@ +#ifndef DUNE_MATRIX_VECTOR_GETSCALAR_HH +#define DUNE_MATRIX_VECTOR_GETSCALAR_HH + +#include <utility> + +#include <dune/common/fmatrix.hh> +#include <dune/istl/scaledidmatrix.hh> + +#include <dune/matrix-vector/traits/utilities.hh> + +namespace Dune { +namespace MatrixVector { + +template <class, class Enable = void> +struct Helper; + +/** + * \brief Extract the scalar of some scale provider. + * \param tol Is the acceptable difference for ambiguous extracted scalings. + **/ +template <class Scalar, class Provider> +static Scalar getScalar(Provider&& p, Scalar tol = 0.0) { + static_assert(isScalar<Scalar>(), "Type not know as scalar."); + + return Helper<std::decay_t<Provider>>::template get<Scalar>( + std::forward<Provider>(p), tol); +} + +//! \brief If the scaling is given by a scalar, that is what we return. +template <class Provider> +struct Helper<Provider, EnableScalar<Provider>> { + template <class Scalar, class Provider_> + static Scalar get(Provider_&& p, Scalar) { + return p; + } +}; + +//! \brief If the scaling is given by a matrix a*Id, extract a. +//! Fail if given matrix is not a multiple of Id. +template <class Provider> +struct Helper<Provider, EnableMatrix<Provider>> { + template <class Scalar, class K, int n> + static Scalar get(const ScaledIdentityMatrix<K, n>& m, Scalar) { + return m.scalar(); + } + + template <class Scalar, class K, int n> + static Scalar get(const FieldMatrix<K, n, n>& m, Scalar tol) { + // TODO make istl-with-checking like mechanism where one avoids all + // the following computations and just returns the first value. + + using namespace std; + Scalar result_min = std::numeric_limits<Scalar>::infinity(); + Scalar result_max = -std::numeric_limits<Scalar>::infinity(); + Scalar others_max = 0; + for (size_t i = 0; i < n; ++i) { + for (size_t j = 0; j < n; ++j) { + if (i == j) { + result_min = min(result_min, m[i][j]); + result_max = max(result_max, m[i][j]); + } else { + others_max = max(others_max, fabs(m[i][j])); + } + } + } + if(others_max > tol) + DUNE_THROW(Exception, "Offdiagonal is not (close enought to) zero."); + if(fabs(result_max - result_min) > tol) + DUNE_THROW(Exception, "Diagonal scaling is (too) ambiguous."); + return result_max; + } + + template <class Scalar, class M> + static Scalar get(const M&, Scalar) { + static_assert(AlwaysFalse<M>::value, "Not yet implemented."); + return 1.0; + } +}; + +} // end namespace MatrixVector +} // end namespace Dune + +#endif // DUNE_MATRIX_VECTOR_GETSCALAR_HH -- GitLab