Merge branch 'feature/float128' into 'master'

Feature/float128

See merge request !518

cmake/modules/DuneCommonMacros.cmake

@@ -28,6 +28,7 @@ set_package_properties("LAPACK" PROPERTIES
 
 find_package(GMP)
 include(AddGMPFlags)
+find_package(QuadMath)
 find_package(Inkscape)
 include(UseInkscape)
 include(FindMProtect)

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})
+
+# register all QuadMath related flags
+if(HAVE_QUADMATH)
+  dune_register_package_flags(COMPILE_DEFINITIONS "ENABLE_QUADMATH=1" "_GLIBCXX_USE_FLOAT128=1"
+                              COMPILE_OPTIONS "-fext-numeric-literals"
+                              LIBRARIES "quadmath")
+endif()

config.h.cmake

@@ -78,6 +78,10 @@
    to facilitate activating and deactivating GMP using compile flags. */
 #cmakedefine HAVE_GMP ENABLE_GMP
 
+/* Define if you have the GCC Quad-Precision library. The value should be ENABLE_QUADMATH
+   to facilitate activating and deactivating QuadMath using compile flags. */
+#cmakedefine HAVE_QUADMATH ENABLE_QUADMATH
+
 /* Define if you have the Vc library. The value should be ENABLE_VC
    to facilitate activating and deactivating Vc using compile flags. */
 #cmakedefine HAVE_VC ENABLE_VC

dune/common/CMakeLists.txt

@@ -92,6 +92,7 @@ install(FILES
         propertymap.hh
         promotiontraits.hh
         proxymemberaccess.hh
+        quadmath.hh
         rangeutilities.hh
         reservedvector.hh
         shared_ptr.hh

dune/common/float_cmp.cc

@@ -71,21 +71,30 @@ namespace Dune {
         static bool eq(const T &first,
                        const T &second,
                        typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T>::value())
-        { return std::abs(first - second) <= epsilon*std::max(std::abs(first), std::abs(second)); }
+        {
+          using std::abs;
+          return abs(first - second) <= epsilon*std::max(abs(first), abs(second));
+        }
       };
       template<class T>
       struct eq_t<T, relativeStrong> {
         static bool eq(const T &first,
                        const T &second,
                        typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T>::value())
-        { return std::abs(first - second) <= epsilon*std::min(std::abs(first), std::abs(second)); }
+        {
+          using std::abs;
+          return abs(first - second) <= epsilon*std::min(abs(first), abs(second));
+        }
       };
       template<class T>
       struct eq_t<T, absolute> {
         static bool eq(const T &first,
                        const T &second,
                        typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T>::value())
-        { return std::abs(first-second) <= epsilon; }
+        {
+          using std::abs;
+          return abs(first-second) <= epsilon;
+        }
       };
       template<class T, CmpStyle cstyle>
       struct eq_t_std_vec {

dune/common/quadmath.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_QUADMATH_HH
+#define DUNE_QUADMATH_HH
+
+#if HAVE_QUADMATH
+#include <quadmath.h>
+
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
+#include <cstdlib> // abs
+#include <istream>
+#include <ostream>
+#include <type_traits>
+#include <utility>
+
+#include <dune/common/typetraits.hh>
+
+namespace Dune
+{
+  namespace Impl
+  {
+    // forward declaration
+    class Float128;
+
+  } // end namespace Impl
+
+  using Impl::Float128;
+
+  // The purpose of this namespace is to move the `<cmath>` function overloads
+  // out of namespace `Dune`, see AlignedNumber in debugalign.hh.
+  namespace Impl
+  {
+    using float128_t = __float128;
+
+    /// Wrapper for quad-precision type __float128
+    class Float128
+    {
+      float128_t value_ = 0.0q;
+
+    public:
+      constexpr Float128() = default;
+      constexpr Float128(const float128_t& value) noexcept
+        : value_(value)
+      {}
+
+      // constructor from any floating-point or integer type
+      template <class T,
+        std::enable_if_t<std::is_arithmetic<T>::value, int> = 0>
+      constexpr Float128(const T& value) noexcept
+        : value_(value)
+      {}
+
+      // accessors
+      constexpr operator float128_t() const noexcept { return value_; }
+
+      constexpr float128_t const& value() const noexcept { return value_; }
+      constexpr float128_t&       value() noexcept       { return value_; }
+
+      // I/O
+      template<class CharT, class Traits>
+      friend std::basic_istream<CharT, Traits>&
+      operator>>(std::basic_istream<CharT, Traits>& in, Float128& x)
+      {
+        std::string buf;
+        buf.reserve(128);
+        in >> buf;
+        x.value() = strtoflt128(buf.c_str(), NULL);
+        return in;
+      }
+
+      template<class CharT, class Traits>
+      friend std::basic_ostream<CharT, Traits>&
+      operator<<(std::basic_ostream<CharT, Traits>& out, const Float128& x)
+      {
+        const std::size_t bufSize = 128;
+        CharT buf[128];
+
+        std::string format = "%." + std::to_string(out.precision()) + "Q" +
+                              ((out.flags() | std::ios_base::scientific) ? "e" : "f");
+        const int numChars = quadmath_snprintf(buf, bufSize, format.c_str(), x.value());
+        if (std::size_t(numChars) >= bufSize) {
+          DUNE_THROW(Dune::RangeError, "Failed to print Float128 value: buffer overflow");
+        }
+        out << buf;
+        return out;
+      }
+
+      // Increment, decrement
+      constexpr Float128& operator++() noexcept { ++value_; return *this; }
+      constexpr Float128& operator--() noexcept { --value_; return *this; }
+
+      constexpr Float128 operator++(int) noexcept { Float128 tmp{*this}; ++value_; return tmp; }
+      constexpr Float128 operator--(int) noexcept { Float128 tmp{*this}; --value_; return tmp; }
+
+      // unary operators
+      constexpr Float128 operator+() const noexcept { return Float128{+value_}; }
+      constexpr Float128 operator-() const noexcept { return Float128{-value_}; }
+
+      // assignment operators
+#define DUNE_ASSIGN_OP(OP)                                              \
+      constexpr Float128& operator OP(const Float128& u) noexcept       \
+      {                                                                 \
+        value_ OP float128_t(u);                                        \
+        return *this;                                                   \
+      }                                                                 \
+      static_assert(true, "Require semicolon to unconfuse editors")
+
+      DUNE_ASSIGN_OP(+=);
+      DUNE_ASSIGN_OP(-=);
+
+      DUNE_ASSIGN_OP(*=);
+      DUNE_ASSIGN_OP(/=);
+
+#undef DUNE_ASSIGN_OP
+
+    }; // end class Float128
+
+    // binary operators:
+    // For symmetry provide overloads with arithmetic types
+    // in the first or second argument.
+#define DUNE_BINARY_OP(OP)                                              \
+    constexpr Float128 operator OP(const Float128& t,                   \
+                                   const Float128& u) noexcept          \
+    {                                                                   \
+      return Float128{float128_t(t) OP float128_t(u)};                  \
+    }                                                                   \
+    template <class T,                                                  \
+      std::enable_if_t<std::is_arithmetic<T>::value, int> = 0>          \
+    constexpr Float128 operator OP(const T& t,                          \
+                                   const Float128& u) noexcept          \
+    {                                                                   \
+      return Float128{float128_t(t) OP float128_t(u)};                  \
+    }                                                                   \
+    template <class U,                                                  \
+      std::enable_if_t<std::is_arithmetic<U>::value, int> = 0>          \
+    constexpr Float128 operator OP(const Float128& t,                   \
+                                   const U& u) noexcept                 \
+    {                                                                   \
+      return Float128{float128_t(t) OP float128_t(u)};                  \
+    }                                                                   \
+    static_assert(true, "Require semicolon to unconfuse editors")
+
+    DUNE_BINARY_OP(+);
+    DUNE_BINARY_OP(-);
+    DUNE_BINARY_OP(*);
+    DUNE_BINARY_OP(/);
+
+#undef DUNE_BINARY_OP
+
+    // logical operators:
+    // For symmetry provide overloads with arithmetic types
+    // in the first or second argument.
+#define DUNE_BINARY_BOOL_OP(OP)                                         \
+    constexpr bool operator OP(const Float128& t,                       \
+                               const Float128& u) noexcept              \
+    {                                                                   \
+      return float128_t(t) OP float128_t(u);                            \
+    }                                                                   \
+    template <class T,                                                  \
+      std::enable_if_t<std::is_arithmetic<T>::value, int> = 0>          \
+    constexpr bool operator OP(const T& t,                              \
+                               const Float128& u) noexcept              \
+    {                                                                   \
+      return float128_t(t) OP float128_t(u);                            \
+    }                                                                   \
+    template <class U,                                                  \
+      std::enable_if_t<std::is_arithmetic<U>::value, int> = 0>          \
+    constexpr bool operator OP(const Float128& t,                       \
+                               const U& u) noexcept                     \
+    {                                                                   \
+      return float128_t(t) OP float128_t(u);                            \
+    }                                                                   \
+    static_assert(true, "Require semicolon to unconfuse editors")
+
+    DUNE_BINARY_BOOL_OP(==);
+    DUNE_BINARY_BOOL_OP(!=);
+    DUNE_BINARY_BOOL_OP(<);
+    DUNE_BINARY_BOOL_OP(>);
+    DUNE_BINARY_BOOL_OP(<=);
+    DUNE_BINARY_BOOL_OP(>=);
+
+#undef DUNE_BINARY_BOOL_OP
+
+    // Overloads for the cmath functions
+
+    // function with name `name` redirects to quadmath function `func`
+#define DUNE_UNARY_FUNC(name,func)                                         \
+    inline Float128 name(const Float128& u) noexcept                       \
+    {                                                                      \
+      return Float128{func (float128_t(u))};                               \
+    }                                                                      \
+    static_assert(true, "Require semicolon to unconfuse editors")
+
+    // like DUNE_UNARY_FUNC but with cutom return type
+#define DUNE_CUSTOM_UNARY_FUNC(type,name,func)                             \
+    inline type name(const Float128& u) noexcept                           \
+    {                                                                      \
+      return (type)(func (float128_t(u)));                                 \
+    }                                                                      \
+    static_assert(true, "Require semicolon to unconfuse editors")
+
+    // redirects to quadmath function with two arguments
+#define DUNE_BINARY_FUNC(name,func)                                        \
+    inline Float128 name(const Float128& t,                                \
+                         const Float128& u) noexcept                       \
+    {                                                                      \
+      return Float128{func (float128_t(t), float128_t(u))};                \
+    }                                                                      \
+    static_assert(true, "Require semicolon to unconfuse editors")
+
+    DUNE_UNARY_FUNC(abs, fabsq);
+    DUNE_UNARY_FUNC(acos, acosq);
+    DUNE_UNARY_FUNC(acosh, acoshq);
+    DUNE_UNARY_FUNC(asin, asinq);
+    DUNE_UNARY_FUNC(asinh, asinhq);
+    DUNE_UNARY_FUNC(atan, atanq);
+    DUNE_UNARY_FUNC(atanh, atanhq);
+    DUNE_UNARY_FUNC(cbrt, cbrtq);
+    DUNE_UNARY_FUNC(ceil, ceilq);
+    DUNE_UNARY_FUNC(cos, cosq);
+    DUNE_UNARY_FUNC(cosh, coshq);
+    DUNE_UNARY_FUNC(erf, erfq);
+    DUNE_UNARY_FUNC(erfc, erfcq);
+    DUNE_UNARY_FUNC(exp, expq);
+    DUNE_UNARY_FUNC(expm1, expm1q);
+    DUNE_UNARY_FUNC(fabs, fabsq);
+    DUNE_UNARY_FUNC(floor, floorq);
+    DUNE_CUSTOM_UNARY_FUNC(int, ilogb, ilogbq);
+    DUNE_UNARY_FUNC(lgamma, lgammaq);
+    DUNE_CUSTOM_UNARY_FUNC(long long int, llrint, llrintq);
+    DUNE_CUSTOM_UNARY_FUNC(long long int, llround, llroundq);
+    DUNE_UNARY_FUNC(log, logq);
+    DUNE_UNARY_FUNC(log10, log10q);
+    DUNE_UNARY_FUNC(log1p, log1pq);
+    DUNE_UNARY_FUNC(log2, log2q);
+    // DUNE_UNARY_FUNC(logb, logbq); // not available in gcc5
+    DUNE_CUSTOM_UNARY_FUNC(long int, lrint, lrintq);
+    DUNE_CUSTOM_UNARY_FUNC(long int, lround, lroundq);
+    DUNE_UNARY_FUNC(nearbyint, nearbyintq);
+    DUNE_BINARY_FUNC(nextafter, nextafterq);
+    DUNE_UNARY_FUNC(rint, rintq);
+    DUNE_UNARY_FUNC(round, roundq);
+    DUNE_UNARY_FUNC(sin, sinq);
+    DUNE_UNARY_FUNC(sinh, sinhq);
+    DUNE_UNARY_FUNC(sqrt, sqrtq);
+    DUNE_UNARY_FUNC(tan, tanq);
+    DUNE_UNARY_FUNC(tanh, tanhq);
+    DUNE_UNARY_FUNC(tgamma, tgammaq);
+    DUNE_UNARY_FUNC(trunc, truncq);
+
+    DUNE_CUSTOM_UNARY_FUNC(bool, isfinite, finiteq);
+    DUNE_CUSTOM_UNARY_FUNC(bool, isinf, isinfq);
+    DUNE_CUSTOM_UNARY_FUNC(bool, isnan, isnanq);
+    DUNE_CUSTOM_UNARY_FUNC(bool, signbit, signbitq);
+
+#undef DUNE_UNARY_FUNC
+#undef DUNE_CUSTOM_UNARY_FUNC
+#undef DUNE_BINARY_FUNC
+
+    // like DUNE_BINARY_FUNC but provide overloads with arithmetic
+    // types in the first or second argument.
+#define DUNE_BINARY_ARITHMETIC_FUNC(name,func)                          \
+    inline Float128 name(const Float128& t,                             \
+                         const Float128& u) noexcept                    \
+    {                                                                   \
+      return Float128{func (float128_t(t), float128_t(u))};             \
+    }                                                                   \
+    template <class T,                                                  \
+      std::enable_if_t<std::is_arithmetic<T>::value, int> = 0>          \
+    inline Float128 name(const T& t,                                    \
+                         const Float128& u) noexcept                    \
+    {                                                                   \
+      return Float128{func (float128_t(t), float128_t(u))};             \
+    }                                                                   \
+    template <class U,                                                  \
+      std::enable_if_t<std::is_arithmetic<U>::value, int> = 0>          \
+    inline Float128 name(const Float128& t,                             \
+                         const U& u) noexcept                           \
+    {                                                                   \
+      return Float128{func (float128_t(t), float128_t(u))};             \
+    }                                                                   \
+    static_assert(true, "Require semicolon to unconfuse editors")
+
+    DUNE_BINARY_ARITHMETIC_FUNC(atan2,atan2q);
+    DUNE_BINARY_ARITHMETIC_FUNC(copysign,copysignq);
+    DUNE_BINARY_ARITHMETIC_FUNC(fdim,fdimq);
+    DUNE_BINARY_ARITHMETIC_FUNC(fmax,fmaxq);
+    DUNE_BINARY_ARITHMETIC_FUNC(fmin,fminq);
+    DUNE_BINARY_ARITHMETIC_FUNC(fmod,fmodq);
+    DUNE_BINARY_ARITHMETIC_FUNC(hypot,hypotq);
+    DUNE_BINARY_ARITHMETIC_FUNC(pow,powq);
+    DUNE_BINARY_ARITHMETIC_FUNC(remainder,remainderq);
+
+#undef DUNE_BINARY_ARITHMETIC_FUNC
+
+    // some more cmath functions with special signature
+
+    inline Float128 fma(const Float128& t, const Float128& u, const Float128& v)
+    {
+      return Float128{fmaq(float128_t(t),float128_t(u),float128_t(v))};
+    }
+
+    inline Float128 frexp(const Float128& u, int* p)
+    {
+      return Float128{frexpq(float128_t(u), p)};
+    }
+
+    inline Float128 ldexp(const Float128& u, int p)
+    {
+      return Float128{ldexpq(float128_t(u), p)};
+    }
+
+    inline Float128 remquo(const Float128& t, const Float128& u, int* quo)
+    {
+      return Float128{remquoq(float128_t(t), float128_t(u), quo)};
+    }
+
+    inline Float128 scalbln(const Float128& u, long int e)
+    {
+      return Float128{scalblnq(float128_t(u), e)};
+    }
+
+    inline Float128 scalbn(const Float128& u, int e)
+    {
+      return Float128{scalbnq(float128_t(u), e)};
+    }
+
+  } // end namespace Impl
+
+  template <>
+  struct IsNumber<Impl::Float128>
+      : public std::true_type {};
+
+} // end namespace Dune
+
+namespace std
+{
+#ifndef NO_STD_NUMERIC_LIMITS_SPECIALIZATION
+  template <>
+  class numeric_limits<Dune::Impl::Float128>
+  {
+    using Float128 = Dune::Impl::Float128;
+    using float128_t = Dune::Impl::float128_t;
+
+  public:
+    static constexpr bool is_specialized = true;
+    static constexpr Float128 min() noexcept { return FLT128_MIN; }
+    static constexpr Float128 max() noexcept { return FLT128_MAX; }
+    static constexpr Float128 lowest() noexcept { return -FLT128_MAX; }
+    static constexpr int digits = FLT128_MANT_DIG;
+    static constexpr int digits10 = FLT128_DIG;
+    static constexpr int max_digits10 = 0;
+    static constexpr bool is_signed = true;
+    static constexpr bool is_integer = false;
+    static constexpr bool is_exact = true;
+    static constexpr int radix = 2;
+    static constexpr Float128 epsilon() noexcept { return FLT128_EPSILON; }
+    static constexpr Float128 round_error() noexcept { return float128_t{0.5}; }
+    static constexpr int min_exponent = FLT128_MIN_EXP;
+    static constexpr int min_exponent10 = FLT128_MIN_10_EXP;
+    static constexpr int max_exponent = FLT128_MAX_EXP;
+    static constexpr int max_exponent10 = FLT128_MAX_10_EXP;
+    static constexpr bool has_infinity = false;
+    static constexpr bool has_quiet_NaN = true;
+    static constexpr bool has_signaling_NaN = false;
+    static constexpr float_denorm_style has_denorm = denorm_present;
+    static constexpr bool has_denorm_loss = false;
+    static constexpr Float128 infinity() noexcept { return float128_t{}; }
+    static Float128 quiet_NaN() noexcept { return nanq(""); }
+    static constexpr Float128 signaling_NaN() noexcept { return float128_t{}; }
+    static constexpr Float128 denorm_min() noexcept { return FLT128_DENORM_MIN; }
+    static constexpr bool is_iec559 = false;
+    static constexpr bool is_bounded = false;
+    static constexpr bool is_modulo = false;
+    static constexpr bool traps = false;
+    static constexpr bool tinyness_before = false;
+    static constexpr float_round_style round_style = round_toward_zero;
+  };
+#endif
+} // end namespace std
+
+#endif // HAVE_QUADMATH
+#endif // DUNE_QUADMATH_HH

dune/common/test/CMakeLists.txt

@@ -246,6 +246,10 @@ dune_add_test(SOURCES pathtest.cc
 dune_add_test(SOURCES poolallocatortest.cc
               LABELS quick)
 
+dune_add_test(SOURCES quadmathtest.cc
+              LINK_LIBRARIES dunecommon
+              CMAKE_GUARD HAVE_QUADMATH)
+
 dune_add_test(SOURCES rangeutilitiestest.cc
               LINK_LIBRARIES dunecommon
               LABELS quick)

dune/common/test/fmatrixtest.cc

@@ -18,6 +18,7 @@
 #include <dune/common/classname.hh>
 #include <dune/common/fmatrix.hh>
 #include <dune/common/ftraits.hh>
+#include <dune/common/quadmath.hh>
 #include <dune/common/rangeutilities.hh>
 #include <dune/common/simd/loop.hh>
 #include <dune/common/simd/simd.hh>
@@ -760,14 +761,25 @@ int main()
     ScalarOperatorTest<float>();
     test_matrix<double, double, double, 1, 1>();
     ScalarOperatorTest<double>();
+#if HAVE_QUADMATH
+    test_matrix<Dune::Float128, Dune::Float128, Dune::Float128, 1, 1>();
+    ScalarOperatorTest<Dune::Float128>();
+#endif
     // test n x m matrices
     test_interface<int, int, 10, 5>();
     test_matrix<int, int, int, 10, 5>();
     test_matrix<double, double, double, 5, 10>();
     test_interface<double, double, 5, 10>();
+#if HAVE_QUADMATH
+    test_matrix<Dune::Float128, Dune::Float128, Dune::Float128, 5, 10>();
+    test_interface<Dune::Float128, Dune::Float128, 5, 10>();
+#endif
     // mixed precision
     test_interface<float, float, 5, 10>();
     test_matrix<float, double, float, 5, 10>();
+#if HAVE_QUADMATH
+    test_matrix<float, double, Dune::Float128, 5, 10>();
+#endif
     // test complex matrices
     test_matrix<std::complex<float>, std::complex<float>, std::complex<float>, 1, 1>();
     test_matrix<std::complex<double>, std::complex<double>, std::complex<double>, 5, 10>();

dune/common/test/fvectortest.cc

@@ -14,6 +14,7 @@
 #include <dune/common/exceptions.hh>
 #include <dune/common/fvector.hh>
 #include <dune/common/gmpfield.hh>
+#include <dune/common/quadmath.hh>
 #include <dune/common/typetraits.hh>
 
 using Dune::FieldVector;
@@ -545,17 +546,33 @@ int main()
     FieldVectorTest<double, 3>();
     FieldVectorTest<long double, 3>();
 #if HAVE_GMP
-    // we skip the complex test and the int test, as these will be very hard to implement with GMPField
-    typedef Dune::GMPField<128u> ft;
-    FieldVectorMainTest<ft,ft,3>();
-    FieldVectorMainTest<ft,ft,2>();
-    FieldVectorMainTest<ft,ft,1>();
-    FieldVectorMainTest<ft,ft,0>();
-    ScalarOperatorTest<ft>();
-    ScalarOrderingTest<ft>();
-    DotProductTest<ft,3>();
+    {
+      // we skip the complex test and the int test, as these will be very hard to implement with GMPField
+      typedef Dune::GMPField<128u> ft;
+      FieldVectorMainTest<ft,ft,3>();
+      FieldVectorMainTest<ft,ft,2>();
+      FieldVectorMainTest<ft,ft,1>();
+      FieldVectorMainTest<ft,ft,0>();
+      ScalarOperatorTest<ft>();
+      ScalarOrderingTest<ft>();
+      DotProductTest<ft,3>();
+    }
 #endif // HAVE_GMP
 
+#if HAVE_QUADMATH
+    {
+      // we skip the int test, as these will be very hard to implement with Float128
+      typedef Dune::Float128 ft;
+      FieldVectorMainTest<ft,ft,3>();
+      FieldVectorMainTest<ft,ft,2>();
+      FieldVectorMainTest<ft,ft,1>();
+      FieldVectorMainTest<ft,ft,0>();
+      ScalarOperatorTest<ft>();
+      ScalarOrderingTest<ft>();
+      DotProductTest<ft,3>();
+    }
+#endif
+
     //test the mathclassifiers Dune::isNaN, Dune::isInf, Dune::isFinite
     fieldvectorMathclassifiersTest();
 

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:
+
+#include "config.h"
+
+#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);
+  auto M4 = M.rightmultiplyany(M2);
+
+  using namespace FMatrixHelp;
+
+  invertMatrix(M,A);
+
+  // test cmath functions for Float128 type
+  using T = Float128;
+  test.check(cmp(abs(T{-1}),T{1}), "abs");
+  test.check(cmp(fabs(T{-1}),T{1}), "fabs");
+
+  test.check(cmp(cos(acos(T{0.5})),T{0.5}), "cos(acos)");
+  test.check(cmp(cosh(acosh(T{1.5})),T{1.5}), "cosh(acosh)");
+  test.check(cmp(sin(asin(T{0.5})),T{0.5}), "sin(asin)");
+  test.check(cmp(sinh(asinh(T{0.5})),T{0.5}), "sinh(asinh)");
+  test.check(cmp(tan(atan(T{0.5})),T{0.5}), "tan(atan)");
+  test.check(cmp(atan2(T{1},T{2}), atan(T{0.5})), "atan2");
+  test.check(cmp(tanh(atanh(T{0.5})),T{0.5}), "tanh(atanh)");
+
+  test.check(cmp(fdim(T{4},T{1}),T{3}), "fdim"); // a > b ? a - b : +0
+  test.check(cmp(fma(T{0.5},T{0.4},T{1.8}),(T{0.5} * T{0.4}) + T{1.8}), "fma");
+  test.check(cmp(fmax(T{0.6},T{0.4}),T{0.6}), "fmax");
+  test.check(cmp(fmin(T{0.6},T{0.4}),T{0.4}), "fmin");
+  test.check(cmp(hypot(T{1.6}, T{2.3}), sqrt(T{1.6}*T{1.6} + T{2.3}*T{2.3})), "hypot");
+  // ilogb
+  test.check(cmp(llrint(T{2.3}),(long long int)(2)), "llrint");
+  test.check(cmp(lrint(T{2.3}),(long int)(2)), "lrint");
+  test.check(cmp(rint(T{2.3}),T{2}), "lrint");
+  test.check(cmp(llround(T{2.3}),(long long int)(2)), "llround");
+  test.check(cmp(lround(T{2.3}),(long int)(2)), "lround");
+  test.check(cmp(round(T{2.3}),T{2}), "round");
+  test.check(cmp(nearbyint(T{2.3}),T{2}), "nearbyint");
+  test.check(cmp(trunc(T{2.7}),T{2}), "trunc");
+  test.check(cmp(ceil(T{1.6}),T{2}), "ceil");
+  test.check(cmp(floor(T{1.6}),T{1}), "floor");
+
+  test.check(cmp(log(exp(T{1.5})),T{1.5}), "log(exp)");
+  test.check(cmp(exp(T{0.2}+T{0.4}), exp(T{0.2})*exp(T{0.4})), "exp"); // exp(a+b) = exp(a)*exp(b)
+  test.check(cmp(expm1(T{0.6}),exp(T{0.6})-T{1}), "expm1");
+  test.check(cmp(log10(T{1000}),T{3}), "log10");
+  test.check(cmp(log2(T{8}),T{3}), "log2");
+  test.check(cmp(log1p(T{1.6}),log(T{1} + T{1.6})), "log1p");
+  // nextafter
+
+  // these two functions produce larger errors
+  test.check(weakcmp(fmod(T{5.1},T{3}),T{2.1}), "fmod");
+  test.check(weakcmp(remainder(T{5.1},T{3}),T{-0.9}), "remainder");
+
+  test.check(cmp(pow(T{2},T{3}),T{8}), "pow");
+  test.check(cmp(cbrt(T{0.5*0.5*0.5}),T{0.5}), "cbrt");
+  test.check(cmp(sqrt(T{4}),T{2}), "sqrt");
+
+  test.check(cmp(erf(T{0}),T{0}), "erf");
+  test.check(cmp(erfc(T{0.6}), T{1}-erf(T{0.6})), "erfc");
+  test.check(cmp(lgamma(T{3}),log(T{2})), "lgamma");
+  test.check(cmp(tgamma(T{3}),T{2}), "tgamma");
+}