From 759523d9848581fdc0961696dda20c269377e7ec Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Fri, 8 Jun 2018 00:54:11 +0200
Subject: [PATCH] wrapper for quad precision floating point type
---
cmake/modules/DuneCommonMacros.cmake | 1 +
cmake/modules/FindQuadMath.cmake | 54 +++++
config.h.cmake | 4 +
dune/common/CMakeLists.txt | 1 +
dune/common/float_cmp.cc | 15 +-
dune/common/fvector.hh | 4 +
dune/common/quadmath.hh | 348 +++++++++++++++++++++++++++
dune/common/test/CMakeLists.txt | 4 +
dune/common/test/fmatrixtest.cc | 12 +
dune/common/test/fvectortest.cc | 2 +
dune/common/test/quadmathtest.cc | 64 +++++
11 files changed, 506 insertions(+), 3 deletions(-)
create mode 100644 cmake/modules/FindQuadMath.cmake
create mode 100644 dune/common/quadmath.hh
create mode 100644 dune/common/test/quadmathtest.cc
diff --git a/cmake/modules/DuneCommonMacros.cmake b/cmake/modules/DuneCommonMacros.cmake
index da58528ed..495610fa5 100644
--- a/cmake/modules/DuneCommonMacros.cmake
+++ b/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)
diff --git a/cmake/modules/FindQuadMath.cmake b/cmake/modules/FindQuadMath.cmake
new file mode 100644
index 000000000..fedb7a86c
--- /dev/null
+++ b/cmake/modules/FindQuadMath.cmake
@@ -0,0 +1,54 @@
+# .. 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()
diff --git a/config.h.cmake b/config.h.cmake
index e06c65b73..09dd40fb9 100644
--- a/config.h.cmake
+++ b/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
diff --git a/dune/common/CMakeLists.txt b/dune/common/CMakeLists.txt
index 2a8884210..521ff5f1b 100644
--- a/dune/common/CMakeLists.txt
+++ b/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
diff --git a/dune/common/float_cmp.cc b/dune/common/float_cmp.cc
index 582da7225..1dc18a89b 100644
--- a/dune/common/float_cmp.cc
+++ b/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 {
diff --git a/dune/common/fvector.hh b/dune/common/fvector.hh
index d614cf8b3..579eaa6f9 100644
--- a/dune/common/fvector.hh
+++ b/dune/common/fvector.hh
@@ -339,6 +339,10 @@ namespace Dune {
/** \brief Const conversion operator */
operator const K& () const { return _data; }
+
+ template <class T,
+ std::enable_if_t<std::is_convertible<K,T>::value && std::is_arithmetic<T>::value, int> = 0>
+ operator T () const { return _data; }
};
/* ----- FV / FV ----- */
diff --git a/dune/common/quadmath.hh b/dune/common/quadmath.hh
new file mode 100644
index 000000000..fc18bc664
--- /dev/null
+++ b/dune/common/quadmath.hh
@@ -0,0 +1,348 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_FLOAT128_HH
+#define DUNE_FLOAT128_HH
+
+#if HAVE_QUADMATH
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
+#include <cstdlib> // abs
+#include <istream>
+#include <ostream>
+#include <type_traits>
+#include <utility>
+#include <quadmath.h>
+
+#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
+ {
+ /// Wrapper for quad-precision type __float128
+ class Float128
+ {
+ __float128 value_ = 0.q;
+
+ public:
+ Float128() = default;
+ Float128(const __float128& value)
+ : value_(value)
+ {}
+
+ // constructor from any floating-point or integer type
+ template <class T,
+ std::enable_if_t<std::is_arithmetic<T>::value, int> = 0>
+ Float128(const T& number)
+ : value_(number)
+ {}
+
+ // accessors
+ operator __float128() const { return value_; }
+
+ __float128 const& value() const { return value_; }
+ __float128& value() { 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
+ Float128& operator++() { ++value_; return *this; }
+ Float128& operator--() { --value_; return *this; }
+
+ Float128 operator++(int) { Float128 tmp{*this}; ++value_; return tmp; }
+ Float128 operator--(int) { Float128 tmp{*this}; --value_; return tmp; }
+
+ // unary operators
+ Float128 operator+() const { return Float128{+value_}; }
+ Float128 operator-() const { return Float128{-value_}; }
+
+ // assignment operators
+#define DUNE_ASSIGN_OP(OP) \
+ Float128& operator OP(const Float128& u) \
+ { \
+ value_ OP u.value(); \
+ 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
+ };
+
+ // binary operators
+#define DUNE_BINARY_OP(OP) \
+ Float128 operator OP(const Float128& t, \
+ const Float128& u) \
+ { \
+ return Float128{__float128(t) OP __float128(u)}; \
+ } \
+ template <class T, \
+ std::enable_if_t<std::is_arithmetic<T>::value, int> = 0> \
+ Float128 operator OP(const T& t, \
+ const Float128& u) \
+ { \
+ return Float128{__float128(t) OP __float128(u)}; \
+ } \
+ template <class U, \
+ std::enable_if_t<std::is_arithmetic<U>::value, int> = 0> \
+ Float128 operator OP(const Float128& t, \
+ const U& u) \
+ { \
+ return Float128{__float128(t) OP __float128(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
+
+#define DUNE_BINARY_BOOL_OP(OP) \
+ bool operator OP(const Float128& t, \
+ const Float128& u) \
+ { \
+ return __float128(t) OP __float128(u); \
+ } \
+ template <class T, \
+ std::enable_if_t<std::is_arithmetic<T>::value, int> = 0> \
+ bool operator OP(const T& t, \
+ const Float128& u) \
+ { \
+ return __float128(t) OP __float128(u); \
+ } \
+ template <class U, \
+ std::enable_if_t<std::is_arithmetic<U>::value, int> = 0> \
+ bool operator OP(const Float128& t, \
+ const U& u) \
+ { \
+ return __float128(t) OP __float128(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
+#define DUNE_UNARY_FUNC(name,func) \
+ Float128 name(const Float128& u) \
+ { \
+ return Float128{func(__float128(u))}; \
+ } \
+ static_assert(true, "Require semicolon to unconfuse editors")
+
+#define DUNE_CUSTOM_UNARY_FUNC(type,name,func) \
+ type name(const Float128& u) \
+ { \
+ return (type)(func(__float128(u))); \
+ } \
+ static_assert(true, "Require semicolon to unconfuse editors")
+
+#define DUNE_BINARY_FUNC(name,func) \
+ Float128 name(const Float128& t, const Float128& u) \
+ { \
+ return Float128{func(__float128(t), __float128(u))}; \
+ } \
+ static_assert(true, "Require semicolon to unconfuse editors")
+
+#define DUNE_TERTIARY_FUNC(name,func) \
+ Float128 name(const Float128&t, const Float128& u, const Float128& v) \
+ { \
+ return Float128{func(__float128(t),__float128(u),__float128(v))}; \
+ } \
+ 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_BINARY_FUNC(atan2, atan2q);
+ DUNE_UNARY_FUNC(atanh, atanhq);
+ DUNE_UNARY_FUNC(cbrt, cbrtq);
+ DUNE_UNARY_FUNC(ceil, ceilq);
+ DUNE_BINARY_FUNC(copysign, copysignq);
+ 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_BINARY_FUNC(fdim, fdimq);
+ DUNE_UNARY_FUNC(floor, floorq);
+ DUNE_TERTIARY_FUNC(fma, fmaq);
+ DUNE_BINARY_FUNC(fmax, fmaxq);
+ DUNE_BINARY_FUNC(fmin, fminq);
+ DUNE_BINARY_FUNC(fmod, fmodq);
+ DUNE_BINARY_FUNC(hypot, hypotq);
+ 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);
+ 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_BINARY_FUNC(pow, powq);
+ DUNE_BINARY_FUNC(remainder, remainderq);
+ 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_BINARY_FUNC
+#undef DUNE_TERTIARY_FUNC
+#undef DUNE_CUSTOM_UNARY_FUNC
+
+ Float128 frexp(const Float128& u, int* p)
+ {
+ return Float128{frexpq(__float128(u),p)};
+ }
+
+ Float128 ldexp(const Float128& u, int p)
+ {
+ return Float128{ldexpq(__float128(u),p)};
+ }
+
+ Float128 nan(const char* arg)
+ {
+ return Float128{nanq(arg)};
+ }
+
+ Float128 remquo(const Float128& t, const Float128& u, int* quo)
+ {
+ return Float128{remquoq(__float128(t),__float128(u),quo)};
+ }
+
+ Float128 scalbln(const Float128& u, long int exp)
+ {
+ return Float128{scalblnq(__float128(u),exp)};
+ }
+
+ Float128 scalbn(const Float128& u, int exp)
+ {
+ return Float128{scalbnq(__float128(u),exp)};
+ }
+
+ } // end namespace Impl
+
+ template <>
+ struct IsNumber<Impl::Float128>
+ : public std::true_type {};
+
+} // namespace Dune
+
+namespace std
+{
+#ifndef NO_STD_NUMERIC_LIMITS_SPECIALIZATION
+ template <>
+ class numeric_limits<Dune::Impl::Float128>
+ {
+ using Float128 = Dune::Impl::Float128;
+
+ public:
+ static const bool is_specialized = true;
+ static Float128 min() { return FLT128_MIN; }
+ static Float128 max() { return FLT128_MAX; }
+ static Float128 lowest() { return -FLT128_MAX; }
+ static const int digits = FLT128_MANT_DIG;
+ static const int digits10 = FLT128_DIG;
+ static const int max_digits10 = 0;
+ static const bool is_signed = true;
+ static const bool is_integer = false;
+ static const bool is_exact = true;
+ static const int radix = 2;
+ static Float128 epsilon() { return FLT128_EPSILON; }
+ static Float128 round_error() { return __float128{0.5}; }
+ static const int min_exponent = FLT128_MIN_EXP;
+ static const int min_exponent10 = FLT128_MIN_10_EXP;
+ static const int max_exponent = FLT128_MAX_EXP;
+ static const int max_exponent10 = FLT128_MAX_10_EXP;
+ static const bool has_infinity = false;
+ static const bool has_quiet_NaN = true;
+ static const bool has_signaling_NaN = false;
+ static const float_denorm_style has_denorm = denorm_present;
+ static const bool has_denorm_loss = false;
+ static Float128 infinity() { return __float128{}; }
+ static Float128 quiet_NaN() { return nanq(""); }
+ static Float128 signaling_NaN() { return __float128{}; }
+ static Float128 denorm_min() { return FLT128_DENORM_MIN; }
+ static const bool is_iec559 = false;
+ static const bool is_bounded = false;
+ static const bool is_modulo = false;
+ static const bool traps = false;
+ static const bool tinyness_before = false;
+ static const float_round_style round_style = round_toward_zero;
+ };
+#endif
+} // end namespace std
+
+
+#endif // HAVE_QUADMATH
+#endif // DUNE_FLOAT128_HH
diff --git a/dune/common/test/CMakeLists.txt b/dune/common/test/CMakeLists.txt
index 0f39de56c..b95fbecc1 100644
--- a/dune/common/test/CMakeLists.txt
+++ b/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)
diff --git a/dune/common/test/fmatrixtest.cc b/dune/common/test/fmatrixtest.cc
index f62f53323..648907f4f 100644
--- a/dune/common/test/fmatrixtest.cc
+++ b/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>();
diff --git a/dune/common/test/fvectortest.cc b/dune/common/test/fvectortest.cc
index f7680901d..f1476fc1e 100644
--- a/dune/common/test/fvectortest.cc
+++ b/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;
@@ -544,6 +545,7 @@ int main()
FieldVectorTest<float, 3>();
FieldVectorTest<double, 3>();
FieldVectorTest<long double, 3>();
+ FieldVectorTest<Dune::Float128, 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;
diff --git a/dune/common/test/quadmathtest.cc b/dune/common/test/quadmathtest.cc
new file mode 100644
index 000000000..973e66f63
--- /dev/null
+++ b/dune/common/test/quadmathtest.cc
@@ -0,0 +1,64 @@
+// -*- 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
+{
+ bool operator()(T const& x, T const& y) { return Dune::FloatCmp::eq(x, y); }
+};
+
+
+int main()
+{
+ // check vector and matrix type with Float128 field type
+ TestSuite test{};
+ Comparator<Float128> cmp{};
+
+ 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);
+}
--
GitLab