Add getScalar to extract the scale factor of some scaling.

dune/matrix-vector/CMakeLists.txt

-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

dune/matrix-vector/getscalar.hh

+#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