Merge branch 'KahanSum' into 'master'

Kahan sum

See merge request ogs/ogs!5065

MathLib/KahanSum.h

@@ -0,0 +1,87 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#pragma once
+
+#include <spdlog/fmt/bundled/ostream.h>
+
+#include <iomanip>
+#include <ostream>
+#include <range/v3/range/concepts.hpp>
+#include <range/v3/range/traits.hpp>
+
+namespace MathLib
+{
+class KahanSum;
+}
+template <>
+struct fmt::formatter<MathLib::KahanSum> : fmt::ostream_formatter
+{
+};
+
+namespace MathLib
+{
+
+class KahanSum
+{
+public:
+ explicit constexpr KahanSum(double const value = 0) : value_(value) {}
+
+ explicit constexpr KahanSum(ranges::range auto const& range) : value_(0)
+ {
+ for (auto const v : range)
+ {
+ *this += v;
+ }
+ }
+
+ constexpr KahanSum operator+(double const increment) const
+ {
+ KahanSum result = *this;
+ return result += increment;
+ }
+
+ constexpr KahanSum operator-(double const increment) const
+ {
+ KahanSum result = *this;
+ return result += -increment;
+ }
+
+ constexpr KahanSum& operator-=(double const increment)
+ {
+ return *this += -increment;
+ }
+
+ constexpr KahanSum& operator+=(double const increment)
+ {
+ double const y = increment - correction_;
+ double const t = value_ + y;
+ correction_ = (t - value_) - y;
+ value_ = t;
+ return *this;
+ }
+
+ constexpr double value() const { return value_; }
+ constexpr double operator()() const { return value_; }
+
+ friend inline std::ostream& operator<<(std::ostream& os, KahanSum const& x)
+ {
+ auto const precision = os.precision();
+ return os << std::setprecision(
+ std::numeric_limits<double>::max_digits10)
+ << x.value() << " (± " << x.correction_ << ')'
+ << std::setprecision(precision);
+ }
+
+private:
+ double value_;
+ double correction_ = 0.;
+};
+
+} // namespace MathLib

Tests/MathLib/KahanSum.cpp

@@ -0,0 +1,71 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#include "MathLib/KahanSum.h"
+
+#include <gtest/gtest.h>
+
+#include <range/v3/view/reverse.hpp>
+
+#include "Tests/AutoCheckTools.h"
+
+using namespace MathLib;
+namespace ac = autocheck;
+
+struct NonNegativeDoubleListGenerator
+{
+ using result_type = std::vector<double>;
+ result_type operator()(std::size_t const size) const
+ {
+ result_type r;
+ r.reserve(size);
+ auto positive_double_gen =
+ ac::map(&ac::absoluteValue, ac::generator<double>{});
+ std::generate_n(std::back_inserter(r), size,
+ ac::fix(size, positive_double_gen));
+ return r;
+ }
+};
+
+struct MathLibKahanSumPropertyTest : public ::testing::Test
+{
+ ac::gtest_reporter gtest_reporter;
+};
+
+TEST_F(MathLibKahanSumPropertyTest, SortedSequences)
+{
+ auto property = [](const std::vector<double>& sequence)
+ {
+ // Sort to get an optimal behaviour for sum computation.
+ std::vector<double> sorted{sequence};
+ std::sort(sorted.begin(), sorted.end(),
+ [](double const a, double const b)
+ { return std::abs(a) < std::abs(b); });
+ KahanSum const a{sorted};
+
+ // Reverse sort yields worst sum behaviour resulting in maximum error.
+ KahanSum const b{sorted | ranges::views::reverse};
+
+ if (a() == b())
+ {
+ return true;
+ }
+ // Testing the relative error to be within the theoretical limits. Keep
+ // in mind, that the sequence is constructed with non-negative numbers
+ // and the condition number is 1.
+ return std::abs(a() - b()) / ((a() + b()) / 2.) <
+ 2. * std::numeric_limits<double>::epsilon();
+ };
+
+ // The non-negative double list generator is used to keep the condition
+ // number, which is defined as sum|a_i|/|sum a_i|, equal to 1.
+ autocheck::check<std::vector<double>>(
+ property, 1000, ac::make_arbitrary(NonNegativeDoubleListGenerator{}),
+ gtest_reporter);
+}