added implementation for quadmath, concurrentcache and strtonumber

cmake/modules/AddQuadMathFlags.cmake

+# Defines the functions to use QuadMath
+#
+# .. cmake_function:: add_dune_quadmath_flags
+#
+#    .. cmake_param:: targets
+#       :positional:
+#       :single:
+#       :required:
+#
+#       A list of targets to use QuadMath with.
+#
+
+
+function(add_dune_quadmath_flags _targets)
+  if(QUADMATH_FOUND)
+    foreach(_target ${_targets})
+      target_link_libraries(${_target} "quadmath")
+      set_property(TARGET ${_target}
+        APPEND_STRING
+        PROPERTY COMPILE_FLAGS "-DENABLE_QUADMATH=1 -D_GLIBCXX_USE_FLOAT128=1 ")
+      if(${CMAKE_CXX_COMPILER_ID} STREQUAL GNU)
+        set_property(TARGET ${_target}
+          APPEND_STRING
+          PROPERTY COMPILE_FLAGS "-fext-numeric-literals ")
+      endif()
+    endforeach(_target ${_targets})
+  endif(QUADMATH_FOUND)
+endfunction(add_dune_quadmath_flags)

cmake/modules/CMakeLists.txt

 set(modules "DuneCommonExtensionsMacros.cmake")
 
+find_package(QuadMath)
+include(AddQuadMathFlags)
+
 install(FILES ${modules} DESTINATION ${DUNE_INSTALL_MODULEDIR})

cmake/modules/FindQuadMath.cmake

+# .. cmake_module::
+#
+#    Find the GCC Quad-Precision library
+#
+#    Sets the following variables:
+#
+#    :code:`QUADMATH_FOUND`
+#       True if the Quad-Precision library was found.
+#
+#
+
+
+# search for the header quadmath.h
+include(CheckIncludeFile)
+check_include_file(quadmath.h QUADMATH_HEADER)
+
+include(CheckCSourceCompiles)
+include(CMakePushCheckState)
+
+cmake_push_check_state() # Save variables
+set(CMAKE_REQUIRED_LIBRARIES quadmath)
+check_c_source_compiles("
+#include <quadmath.h>
+
+int main ()
+{
+  __float128 r = 1.0q;
+  r = strtoflt128(\"1.2345678\", NULL);
+  return 0;
+}" QUADMATH_COMPILES)
+cmake_pop_check_state()
+
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(
+  "QuadMath"
+  DEFAULT_MSG
+  QUADMATH_HEADER
+  QUADMATH_COMPILES
+)
+
+# text for feature summary
+set_package_properties("QuadMath" PROPERTIES
+  DESCRIPTION "GCC Quad-Precision library")
+
+# set HAVE_QUADMATH for config.h
+set(HAVE_QUADMATH ${QUADMATH_FOUND})
+
+# -fext-numeric-literals is a GCC extension not available in other compilers like clang
+if(${CMAKE_CXX_COMPILER_ID} STREQUAL GNU)
+  set(_QUADMATH_EXT_NUMERIC_LITERALS "-fext-numeric-literals")
+endif()
+
+# register all QuadMath related flags
+if(HAVE_QUADMATH)
+  dune_register_package_flags(COMPILE_DEFINITIONS "ENABLE_QUADMATH=1" "_GLIBCXX_USE_FLOAT128=1"
+                              COMPILE_OPTIONS ${_QUADMATH_EXT_NUMERIC_LITERALS}
+                              LIBRARIES "quadmath")
+endif()

dune/common/CMakeLists.txt

 add_subdirectory(std)
 add_subdirectory(test)
+
+install(FILES
+  concurrentcache.hh
+  quadmath.hh
+  strtonumber.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/common)
\ No newline at end of file

dune/common/concurrentcache.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_COMMON_CONCURRENT_CACHE_HH
+#define DUNE_COMMON_CONCURRENT_CACHE_HH
+
+#include <mutex>
+#include <shared_mutex>
+#include <thread>
+#include <tuple>
+#include <unordered_map>
+
+#include <dune/common/hash.hh>
+#include <dune/common/std/type_traits.hh>
+
+namespace Dune
+{
+  /// Store cache in instance.
+  template <class Container>
+  struct ConsecutivePolicy;
+
+  /// Store cache thread local, requires no locking.
+  template <class Container>
+  struct ThreadLocalPolicy;
+
+  /// Stores cache global static, requires locking on write access.
+  template <class Container>
+  struct StaticLockedPolicy;
+
+
+  /// \brief The class template ConcurrentCache describes an associative static container that allows the
+  /// concurrent access to the stored data.
+  /**
+   * Cache data of arbitray type that needs initialization on the first access. The data is thereby
+   * initialized thread-wise or globally only once, and guarantees that you always get initialized data.
+   *
+   * \tparam Key        The type of key to access the data.
+   * \tparam Data       The type of the data stored in the cache. The behaviur is undefined if Data is not
+   *                    the same type as Container::mapped_type.
+   * \tparam Policy     A policy class template implementing the method `get_or_init()`. Three implementations
+   *                    are provided: \ref ConsecutivePolicy, \ref ThreadLocalPolicy and \ref StaticLockedPolicy.
+   *                    By default, if no policy class template is specified, the `ThreadLocalPolicy` is used.
+   *                    \see ConcurrentCachePolicy
+   * \tparam Container  The type of the underlying associative container to use to store the data. The
+   *                    container must satisfy the requirements of AssociativeContainer. The standard
+   *                    containers `std::map` and `std::unordered_map` satisfie this requirement. By default,
+   *                    if not container class is specified, the standard container `std::unordered_map<Key,Data>`
+   *                    is used. Note, an unordered_map requires the key to be hashable.
+   *
+   * The `Policy` class template is a template parametrizable with the container type, that provides a static `get_or_init()`
+   * method that is called with the key, and a functor for creation of new data elements.
+   **/
+  template <class Key,
+            class Data,
+            template <class> class Policy = ThreadLocalPolicy,
+            class Container = std::unordered_map<Key, Data>>
+  class ConcurrentCache;
+
+
+#ifdef DOXYGEN
+  /// \brief The class template ConcurrentCachePolicy describes a concrete policies for the use in \ref ConcurrentCache.
+  /**
+   * Provide a static cache and a `get_or_init()` static method that extracts the data from the cache if it exists or
+   * creates a new extry by using an initialization functor.
+   *
+   * Realizations of this template are \ref ConsecutivePolicy, \ref ThreadLocalPolicy and \ref StaticLockedPolicy.
+   *
+   * \tparam Container  The Type of the associative container key->data to store the cached data.
+   **/
+  template <class Container>
+  class ConcurrentCachePolicy;
+#endif
+
+  // implementation of the consecutive policy. Data is stored in instance variable.
+  template <class Container>
+  struct ConsecutivePolicy
+  {
+    using key_type = typename Container::key_type;
+    using data_type = typename Container::mapped_type;
+    using container_type = Container;
+
+    template <class F, class... Args>
+    data_type const& get_or_init(key_type const& key, F&& f, Args&&... args) const
+    {
+      return impl(std::is_default_constructible<data_type>{},
+        key, std::forward<F>(f), std::forward<Args>(args)...);
+    }
+
+  private:
+    // data_type is default_constructible
+    template <class F, class... Args>
+    data_type const& impl(std::true_type, key_type const& key, F&& f, Args&&... args) const
+    {
+      data_type empty;
+      auto it = cachedData_.emplace(key, std::move(empty));
+      if (it.second) {
+        data_type data = f(key, std::forward<Args>(args)...);
+        it.first->second = std::move(data);
+      }
+      return it.first->second;
+    }
+
+    // data_type is not default_constructible
+    template <class F, class... Args>
+    data_type const& impl(std::false_type, key_type const& key, F&& f, Args&&... args) const
+    {
+      auto it = cachedData_.find(key);
+      if (it != cachedData_.end())
+        return it->second;
+      else {
+        data_type data = f(key, std::forward<Args>(args)...);
+        auto it = cachedData_.emplace(key, std::move(data));
+        return it.first->second;
+      }
+    }
+
+    mutable container_type cachedData_;
+  };
+
+  // implementation of the ThreadLocal policy. Data is stored in thread_local variable.
+  template <class Container>
+  struct ThreadLocalPolicy
+  {
+    using key_type = typename Container::key_type;
+    using data_type = typename Container::mapped_type;
+    using container_type = Container;
+
+    template <class F, class... Args>
+    static data_type const& get_or_init(key_type const& key, F&& f, Args&&... args)
+    {
+      return impl(std::is_default_constructible<data_type>{},
+        key, std::forward<F>(f), std::forward<Args>(args)...);
+    }
+
+  private:
+    // data_type is default_constructible
+    template <class F, class... Args>
+    static data_type const& impl(std::true_type, key_type const& key, F&& f, Args&&... args)
+    {
+      // Container to store the cached values
+      thread_local container_type cached_data;
+
+      data_type empty;
+      auto it = cached_data.emplace(key, std::move(empty));
+      if (it.second) {
+        data_type data = f(key, std::forward<Args>(args)...);
+        it.first->second = std::move(data);
+      }
+      return it.first->second;
+    }
+
+    // data_type is not default_constructible
+    template <class F, class... Args>
+    static data_type const& impl(std::false_type, key_type const& key, F&& f, Args&&... args)
+    {
+      // Container to store the cached values
+      thread_local container_type cached_data;
+
+      auto it = cached_data.find(key);
+      if (it != cached_data.end())
+        return it->second;
+      else {
+        data_type data = f(key, std::forward<Args>(args)...);
+        auto it = cached_data.emplace(key, std::move(data));
+        return it.first->second;
+      }
+    }
+  };
+
+
+  // implementation of the Shared policy. Data is stored in static variable.
+  template <class Container>
+  struct StaticLockedPolicy
+  {
+    using key_type = typename Container::key_type;
+    using data_type = typename Container::mapped_type;
+    using container_type = Container;
+
+    template <class F, class... Args>
+    static data_type const& get_or_init(key_type const& key, F&& f, Args&&... args)
+    {
+      // Container to store the cached values
+      static container_type cached_data;
+
+      // mutex used to access the data in the container, necessary since
+      // access emplace is read-write.
+      using mutex_type = std::shared_timed_mutex;
+      static mutex_type access_mutex;
+
+      // first try to lock for read-only, if an element for key is found, return it,
+      // if not, obtain a unique_lock to insert a new element and initialize it.
+      std::shared_lock<mutex_type> read_lock(access_mutex);
+      auto it = cached_data.find(key);
+      if (it != cached_data.end())
+        return it->second;
+      else {
+        read_lock.unlock();
+        data_type data = f(key, std::forward<Args>(args)...);
+        std::unique_lock<mutex_type> write_lock(access_mutex);
+        auto new_it = cached_data.emplace(key, std::move(data));
+        return new_it.first->second;
+      }
+    }
+  };
+
+
+  template <class Key, class Data, template <class> class Policy, class Container>
+  class ConcurrentCache
+      : protected Policy<Container>
+  {
+    using key_type = Key;
+    using data_type = Data;
+
+  public:
+
+    /// \brief Return the data associated to the `key`.
+    /**
+     * Return the data associated to key. If no data is found, create a new entry in the container
+     * with a value obtained from the functor, by calling `f(key, args...)`.
+     *
+     * \param f        A functor of signature data_type(key_type, Args...)
+     * \param args...  Arguments passed additionally to the functor f
+     **/
+    template <class F, class... Args>
+    data_type const& get(key_type const& key, F&& f, Args&&... args) const
+    {
+      static_assert(Std::is_callable<F(key_type, Args...), data_type>::value,
+        "Functor F must have the signature data_type(key_type, Args...)");
+
+      return ConcurrentCache::get_or_init(key, std::forward<F>(f), std::forward<Args>(args)...);
+    }
+  };
+
+} // end namespace Dune
+
+#endif // DUNE_COMMON_CONCURRENT_CACHE_HH

dune/common/strtonumber.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_COMMON_STRTONUMBER_HH
+#define DUNE_COMMON_STRTONUMBER_HH
+
+#include <cerrno>
+#include <cstdlib>
+#include <limits>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/std/charconv.hh>
+#include <dune/common/std/type_traits.hh>
+
+namespace Dune
+{
+  /// An exception thrown whenever interpretation of the argument is not possible.
+  class InvalidArgument : public Exception {};
+
+  namespace Impl
+  {
+    // forward declaration of default parser
+    template <class T>
+    struct FromCharsParser;
+
+    template <class T>
+    using HasFromChars
+      = decltype(Std::from_chars(std::declval<const char*>(), std::declval<const char*>(), std::declval<T&>()));
+
+
+    // default behavior: call the string constructor
+    template<typename T, class = void>
+    struct StrToNumber
+    {
+      static T eval (const char* first, const char* /*last*/)
+      {
+        return T(first);;
+      }
+    };
+
+    // if possible, use std::from_chars to convert string to number
+    template<typename T>
+    struct StrToNumber<T, std::enable_if_t<Std::is_detected<HasFromChars, T>::value>>
+    {
+      static T eval (const char* first, const char* last)
+      {
+        return FromCharsParser<T>::eval(first, last);
+      }
+    };
+
+
+    /// \brief Implementation of a parser using std::from_chars functions
+    /**
+     * \tparam T  Target numeric type
+     **/
+    template<typename T>
+    struct FromCharsParser
+    {
+      static T eval (const char* first, const char* last)
+      {
+        T value;
+        auto result = Std::from_chars(first, last, value);
+
+        if (result.ec == std::errc::result_out_of_range) {
+          DUNE_THROW(RangeError, std::error_condition(result.ec).message());
+        }
+        else if (result.ec == std::errc::invalid_argument) {
+          DUNE_THROW(InvalidArgument,
+            "Conversion of '" << first << "' to number failed. Possible reason: invalid string or locale format.");
+        }
+        else if (result.ec == std::errc::not_supported) {
+          DUNE_THROW(NotImplemented,
+            "Parameter passed to from_chars not supported.");
+        }
+
+        return value;
+      }
+    };
+
+    /// \brief Implementation of a fallback parser using std::strtoXXX like functions
+    /**
+     * Allows leading and trailing whitespace characters.
+     * NOTE: This parser is not locale independent, but throws an error if whole
+     *       string can not be parsed completely.
+     *
+     * \tparam T  Target numeric type
+     **/
+    template<typename T>
+    struct StrToNumberParser
+    {
+      // The parser has the signature `T(const char*, char**)` and may set the errno
+      // in case of a range error.
+      template<typename Parser>
+      T eval (const char* first, const char* /*last*/, Parser parser)
+      {
+        char* end;
+        auto old_errno = errno;
+        errno = 0;
+        auto x = parser(first, &end); // maybe sets errno
+        std::swap(errno, old_errno);
+
+        if (old_errno == ERANGE) {
+          DUNE_THROW(RangeError, std::strerror(old_errno));
+        }
+
+        // test whether all non-space characters are consumed during conversion
+        bool all_consumed = (end != first);
+        while (all_consumed && (*end != '\0'))
+          all_consumed = std::isspace(*end++);
+
+        if (!all_consumed) {
+          DUNE_THROW(InvalidArgument,
+            "Conversion of '" << first << "' to number failed. Possible reason: invalid string or locale format");
+        }
+
+        return convertToRange<T>(x);
+      }
+
+    private:
+      // Check whether a numeric conversion is safe
+      template<typename U>
+      static T convertToRange (U const& u)
+      {
+        if (sizeof(U) > sizeof(T) && (u > U(std::numeric_limits<T>::max()) ||
+                                      u < U(std::numeric_limits<T>::min()))) {
+          DUNE_THROW(RangeError, "Numerical result out of range");
+        }
+        return T(u);
+      }
+    };
+
+  } // end namespace Impl
+
+
+  /// \brief Convert a character sequence to a number type `T`
+  /**
+   * The cast from character sequence to any numeric type is implemented in terms of the std library function
+   * `from_chars()` but parametrized with the concrete target type for ease of usage
+   * in generic code.
+   *
+   * Only types that provide a specialization of `Impl::StrToNumber` can call the special functions
+   * for conversion. All other types by default call the constructor with strings.
+   *
+   * The conversion is locale-independent if an overload of from_chars() is found and
+   * throws in the case of the fallback implementation an \ref InvalidArgument exception if not all
+   * characters are consumed during conversion, except leading and trailing whitespaces.
+   *
+   * In case the represented number is out of range of the number type T, a \ref RangeError exception
+   * is thrown.
+   *
+   * \tparam T   The target number type to convert the string to.
+   * \param str  A pointer to the null-terminated byte string to be interpreted.
+   *
+   * \throws InvalidArgument
+   * \throws RangeError
+   *
+   * Example of usage:
+   * \code{.cpp}
+   * double x = strTo<double>("1.2345");
+   * \endcode
+   **/
+  template<typename T>
+  T strTo (const char* str)
+  {
+    return Impl::StrToNumber<T>::eval(str, str + std::strlen(str));
+  }
+
+  /// Overload of \ref strTo for `std::string` arguments.
+  template<typename T>
+  T strTo (const std::string& str)
+  {
+    return Impl::StrToNumber<T>::eval(str.c_str(), str.c_str() + str.size());
+  }
+
+} // end namespace Dune
+
+#endif // DUNE_COMMON_STRTONUMBER_HH

dune/common/test/CMakeLists.txt

 dune_add_test(SOURCES charconvtest.cc
+              LINK_LIBRARIES dunecommon)
+
+dune_add_test(SOURCES concurrentcachetest.cc
+              LINK_LIBRARIES dunecommon)
+
+dune_add_test(SOURCES quadmathtest.cc
               LINK_LIBRARIES dunecommon
-              LABELS quick)
+              CMAKE_GUARD HAVE_QUADMATH)
+
+dune_add_test(SOURCES strtonumbertest.cc
+              LINK_LIBRARIES dunecommon)

dune/common/test/charconvtest.cc

@@ -167,7 +167,7 @@ int main()
         auto result3 = std::from_chars(bufPtr, result1.ptr, new_value2, base);
         test.check(bool(new_value == new_value2), "comparison with std::from_chars");
 #else
-        std::setlocale(LC_NUMERIC, "C");
+        setlocale(LC_NUMERIC, "C");
         T new_value2;
         if (std::is_unsigned<T>::value)
           new_value2 = std::strtoull(bufPtr, nullptr, base);

dune/common/test/concurrentcachetest.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <algorithm>
+#include <array>
+#include <iostream>
+#include <mutex>
+#include <random>
+#include <sstream>
+#include <thread>
+
+#include <dune/common/concurrentcache.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/hash.hh>
+
+namespace Dune
+{
+  struct quadrature_key
+  {
+    int id; // topologyId
+    int p;  // order
+    int qt; // quadrature type
+
+    struct hasher
+    {
+      std::size_t operator()(quadrature_key const& t) const
+      {
+        std::size_t seed = 0;
+        hash_combine(seed, t.id);
+        hash_combine(seed, t.p);
+        hash_combine(seed, t.qt);
+        return seed;
+      }
+    };
+
+    friend bool operator==(quadrature_key const& lhs, quadrature_key const& rhs)
+    {
+      return std::tie(lhs.id,lhs.p,lhs.qt) == std::tie(rhs.id,rhs.p,rhs.qt);
+    }
+  };
+
+} // end namespace Dune
+
+
+using quadrature_key = Dune::quadrature_key;
+using quadrature_data = std::vector<double>;
+
+quadrature_data init_data(quadrature_key const& key)
+{
+  quadrature_data data(100);
+  std::generate(data.begin(), data.end(), []{ return std::fmod(double(std::rand()), 10.0); });
+  std::stringstream ss;
+  ss << "init [" << key.id << "," << key.p << "," << key.qt << "]\n";
+  std::cout << ss.str();
+  return data;
+}
+
+template <template <class> class Policy, unsigned threads_count = 16>
+void test()
+{
+  auto cache = Dune::ConcurrentCache<quadrature_key, quadrature_data, Policy,
+    std::unordered_map<quadrature_key,quadrature_data,quadrature_key::hasher>>{};
+
+  std::random_device rd;
+
+  std::thread threads[threads_count];
+  for (auto& t: threads) {
+    t = std::thread([&rd,&cache]() {
+      std::mt19937 gen(rd());
+      std::uniform_int_distribution<int> uniform_dist(1, 3);
+      for (auto i = 0; i < 100; ++i) {
+        int id = uniform_dist(gen);
+        int p = uniform_dist(gen);
+        int qt = uniform_dist(gen);
+
+        auto const& data = cache.get(quadrature_key{id,p,qt}, init_data);
+
+        if (data.size() != 100) {
+          DUNE_THROW(Dune::Exception, "Data must be initialized to size 100");
+        }
+      }
+    });
+  }
+
+  for (auto& t: threads) {
+    t.join();
+  }
+}
+
+int main()
+{
+  test<Dune::ThreadLocalPolicy>();
+  test<Dune::StaticLockedPolicy>();
+  test<Dune::ConsecutivePolicy,1>();
+}

dune/common/test/quadmathtest.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <cassert>
+#include <iostream>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/float_cmp.hh>
+#include <dune/common/quadmath.hh>
+#include <dune/common/test/testsuite.hh>
+
+using namespace Dune;
+
+template <class T>
+struct Comparator
+{
+  Comparator(T tol)
+    : tol_(tol)
+  {}
+
+  bool operator()(T const& x, T const& y)
+  {
+    return Dune::FloatCmp::eq<T,FloatCmp::absolute>(x, y, tol_);
+  }
+
+private:
+  T tol_;
+};
+
+int main()
+{
+  // check vector and matrix type with Float128 field type
+  TestSuite test{};
+  Comparator<Float128> cmp{std::numeric_limits<Float128>::epsilon() * 8};
+  Comparator<Float128> weakcmp{cbrt(std::numeric_limits<Float128>::epsilon())};
+
+  // implicit conversion
+  Float128 x1 = 1;
+  Float128 x2 = 1.0f;
+  Float128 x3 = 1.0;
+  Float128 x4 = 1.0l;
+
+  int z1 = x1;
+  float z2 = x2;
+  double z3 = x3;
+  long double z4 = x4;
+
+  // field-vector
+  FieldVector<Float128,3> v{1,2,3}, x;
+  FieldMatrix<Float128,3,3> M{ {1,2,3}, {2,3,4}, {3,4,6} }, A;
+  FieldMatrix<Float128,3,3> M2{ {1,2,3}, {2,3,4}, {3,4,7} };
+
+  auto y1 = v.one_norm();
+  test.check(cmp(y1, 6.q), "vec.one_norm()");
+
+  auto y2 = v.two_norm();
+  test.check(cmp(y2, sqrtq(14.q)), "vec.two_norm()");
+
+  auto y3 = v.infinity_norm();
+  test.check(cmp(y3, 3.q), "vec.infinity_norm()");
+
+  M.mv(v, x);   // x = M*v
+  M.mtv(v, x);  // x = M^T*v
+  M.umv(v, x);  // x+= M*v
+  M.umtv(v, x); // x+= M^T*v
+  M.mmv(v, x);  // x-= M*v
+  M.mmtv(v, x); // x-= M^T*v
+
+  auto w1 = M.infinity_norm();
+  test.check(cmp(w1, 13.q), "mat.infinity_norm()");
+
+  auto w2 = M.determinant();
+  test.check(cmp(w2, -1.q), "mat.determinant()");
+
+  M.solve(v, x);  // x = M^(-1)*v
+
+  auto M3 = M.leftmultiplyany(M2);
+