[flang-commits] [flang] c375ec8 - [flang] Use a "double-double" accumulator in SUM

[peter klausler via flang-commits]

Author: peter klausler
Date: 2021-06-16T14:29:39-07:00
New Revision: c375ec86132984bc3915323fbe3f8c71f4d5503f

URL: https://github.com/llvm/llvm-project/commit/c375ec86132984bc3915323fbe3f8c71f4d5503f
DIFF: https://github.com/llvm/llvm-project/commit/c375ec86132984bc3915323fbe3f8c71f4d5503f.diff

LOG: [flang] Use a "double-double" accumulator in SUM

Use a "double-double" accumulator, a/k/a Kahan summation,
in the SUM intrinsic in the runtime for real & complex.
This seems to be the best-recommended technique for reducing
error, as opposed to the initial implementation of SUM's
distinct accumulators for positive and negative items.

Differential Revision: https://reviews.llvm.org/D104338

Added: 
    

Modified: 
    flang/runtime/sum.cpp

Removed: 
    


################################################################################
diff  --git a/flang/runtime/sum.cpp b/flang/runtime/sum.cpp
index db808b2b4c32..fb3250541953 100644
--- a/flang/runtime/sum.cpp
+++ b/flang/runtime/sum.cpp
@@ -9,8 +9,8 @@
 // Implements SUM for all required operand types and shapes.
 //
 // Real and complex SUM reductions attempt to reduce floating-point
-// cancellation on intermediate results by adding up partial sums
-// for positive and negative elements independently.
+// cancellation on intermediate results by using "Kahan summation"
+// (basically the same as manual "double-double").
 
 #include "reduction-templates.h"
 #include "reduction.h"
@@ -40,24 +40,17 @@ template <typename INTERMEDIATE> class IntegerSumAccumulator {
 template <typename INTERMEDIATE> class RealSumAccumulator {
 public:
   explicit RealSumAccumulator(const Descriptor &array) : array_{array} {}
-  void Reinitialize() { positives_ = negatives_ = inOrder_ = 0; }
-  template <typename A> A Result() const {
-    auto sum{static_cast<A>(positives_ + negatives_)};
-    return std::isfinite(sum) ? sum : static_cast<A>(inOrder_);
-  }
+  void Reinitialize() { sum_ = correction_ = 0; }
+  template <typename A> A Result() const { return sum_; }
   template <typename A> void GetResult(A *p, int /*zeroBasedDim*/ = -1) const {
     *p = Result<A>();
   }
   template <typename A> bool Accumulate(A x) {
-    // Accumulate the nonnegative and negative elements independently
-    // to reduce cancellation; also record an in-order sum for use
-    // in case of overflow.
-    if (x >= 0) {
-      positives_ += x;
-    } else {
-      negatives_ += x;
-    }
-    inOrder_ += x;
+    // Kahan summation
+    auto next{x + correction_};
+    auto oldSum{sum_};
+    sum_ += next;
+    correction_ = (sum_ - oldSum) - next; // algebraically zero
     return true;
   }
   template <typename A> bool AccumulateAt(const SubscriptValue at[]) {
@@ -66,7 +59,7 @@ template <typename INTERMEDIATE> class RealSumAccumulator {
 
 private:
   const Descriptor &array_;
-  INTERMEDIATE positives_{0.0}, negatives_{0.0}, inOrder_{0.0};
+  INTERMEDIATE sum_{0.0}, correction_{0.0};
 };
 
 template <typename PART> class ComplexSumAccumulator {