Merge branch 'KahanSum' into 'master'

Kahan sum

See merge request ogs/ogs!5065

MathLib/KahanSum.h

+/**
+ * \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

+/**
+ * \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);
+}