flang/runtime/sum.cpp

   1 //===-- runtime/sum.cpp ---------------------------------------------------===//
   2 //
   3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
   4 // See https://llvm.org/LICENSE.txt for license information.
   5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
   6 //
   7 //===----------------------------------------------------------------------===//
   8
   9 // Implements SUM for all required operand types and shapes.
  10 //
  11 // Real and complex SUM reductions attempt to reduce floating-point
  12 // cancellation on intermediate results by using "Kahan summation"
  13 // (basically the same as manual "double-double").
  14
  15 #include "reduction-templates.h"
  16 #include "flang/Common/float128.h"
  17 #include "flang/Runtime/reduction.h"
  18 #include <cfloat>
  19 #include <cinttypes>
  20 #include <complex>
  21
  22 namespace Fortran::runtime {
  23
  24 template <typename INTERMEDIATE> class IntegerSumAccumulator {
  25 public:
  26   explicit IntegerSumAccumulator(const Descriptor &array) : array_{array} {}
  27   void Reinitialize() { sum_ = 0; }
  28   template <typename A> void GetResult(A *p, int /*zeroBasedDim*/ = -1) const {
  29     *p = static_cast<A>(sum_);
  30   }
  31   template <typename A> bool AccumulateAt(const SubscriptValue at[]) {
  32     sum_ += *array_.Element<A>(at);
  33     return true;
  34   }
  35
  36 private:
  37   const Descriptor &array_;
  38   INTERMEDIATE sum_{0};
  39 };
  40
  41 template <typename INTERMEDIATE> class RealSumAccumulator {
  42 public:
  43   explicit RealSumAccumulator(const Descriptor &array) : array_{array} {}
  44   void Reinitialize() { sum_ = correction_ = 0; }
  45   template <typename A> A Result() const { return sum_; }
  46   template <typename A> void GetResult(A *p, int /*zeroBasedDim*/ = -1) const {
  47     *p = Result<A>();
  48   }
  49   template <typename A> bool Accumulate(A x) {
  50     // Kahan summation
  51     auto next{x + correction_};
  52     auto oldSum{sum_};
  53     sum_ += next;
  54     correction_ = (sum_ - oldSum) - next; // algebraically zero
  55     return true;
  56   }
  57   template <typename A> bool AccumulateAt(const SubscriptValue at[]) {
  58     return Accumulate(*array_.Element<A>(at));
  59   }
  60
  61 private:
  62   const Descriptor &array_;
  63   INTERMEDIATE sum_{0.0}, correction_{0.0};
  64 };
  65
  66 template <typename PART> class ComplexSumAccumulator {
  67 public:
  68   explicit ComplexSumAccumulator(const Descriptor &array) : array_{array} {}
  69   void Reinitialize() {
  70     reals_.Reinitialize();
  71     imaginaries_.Reinitialize();
  72   }
  73   template <typename A> void GetResult(A *p, int /*zeroBasedDim*/ = -1) const {
  74     using ResultPart = typename A::value_type;
  75     *p = {reals_.template Result<ResultPart>(),
  76         imaginaries_.template Result<ResultPart>()};
  77   }
  78   template <typename A> bool Accumulate(const A &z) {
  79     reals_.Accumulate(z.real());
  80     imaginaries_.Accumulate(z.imag());
  81     return true;
  82   }
  83   template <typename A> bool AccumulateAt(const SubscriptValue at[]) {
  84     return Accumulate(*array_.Element<A>(at));
  85   }
  86
  87 private:
  88   const Descriptor &array_;
  89   RealSumAccumulator<PART> reals_{array_}, imaginaries_{array_};
  90 };
  91
  92 extern "C" {
  93 CppTypeFor<TypeCategory::Integer, 1> RTNAME(SumInteger1)(const Descriptor &x,
  94     const char *source, int line, int dim, const Descriptor *mask) {
  95   return GetTotalReduction<TypeCategory::Integer, 1>(x, source, line, dim, mask,
  96       IntegerSumAccumulator<CppTypeFor<TypeCategory::Integer, 4>>{x}, "SUM");
  97 }
  98 CppTypeFor<TypeCategory::Integer, 2> RTNAME(SumInteger2)(const Descriptor &x,
  99     const char *source, int line, int dim, const Descriptor *mask) {
 100   return GetTotalReduction<TypeCategory::Integer, 2>(x, source, line, dim, mask,
 101       IntegerSumAccumulator<CppTypeFor<TypeCategory::Integer, 4>>{x}, "SUM");
 102 }
 103 CppTypeFor<TypeCategory::Integer, 4> RTNAME(SumInteger4)(const Descriptor &x,
 104     const char *source, int line, int dim, const Descriptor *mask) {
 105   return GetTotalReduction<TypeCategory::Integer, 4>(x, source, line, dim, mask,
 106       IntegerSumAccumulator<CppTypeFor<TypeCategory::Integer, 4>>{x}, "SUM");
 107 }
 108 CppTypeFor<TypeCategory::Integer, 8> RTNAME(SumInteger8)(const Descriptor &x,
 109     const char *source, int line, int dim, const Descriptor *mask) {
 110   return GetTotalReduction<TypeCategory::Integer, 8>(x, source, line, dim, mask,
 111       IntegerSumAccumulator<CppTypeFor<TypeCategory::Integer, 8>>{x}, "SUM");
 112 }
 113 #ifdef __SIZEOF_INT128__
 114 CppTypeFor<TypeCategory::Integer, 16> RTNAME(SumInteger16)(const Descriptor &x,
 115     const char *source, int line, int dim, const Descriptor *mask) {
 116   return GetTotalReduction<TypeCategory::Integer, 16>(x, source, line, dim,
 117       mask, IntegerSumAccumulator<CppTypeFor<TypeCategory::Integer, 16>>{x},
 118       "SUM");
 119 }
 120 #endif
 121
 122 // TODO: real/complex(2 & 3)
 123 CppTypeFor<TypeCategory::Real, 4> RTNAME(SumReal4)(const Descriptor &x,
 124     const char *source, int line, int dim, const Descriptor *mask) {
 125   return GetTotalReduction<TypeCategory::Real, 4>(
 126       x, source, line, dim, mask, RealSumAccumulator<double>{x}, "SUM");
 127 }
 128 CppTypeFor<TypeCategory::Real, 8> RTNAME(SumReal8)(const Descriptor &x,
 129     const char *source, int line, int dim, const Descriptor *mask) {
 130   return GetTotalReduction<TypeCategory::Real, 8>(
 131       x, source, line, dim, mask, RealSumAccumulator<double>{x}, "SUM");
 132 }
 133 #if LDBL_MANT_DIG == 64
 134 CppTypeFor<TypeCategory::Real, 10> RTNAME(SumReal10)(const Descriptor &x,
 135     const char *source, int line, int dim, const Descriptor *mask) {
 136   return GetTotalReduction<TypeCategory::Real, 10>(
 137       x, source, line, dim, mask, RealSumAccumulator<long double>{x}, "SUM");
 138 }
 139 #endif
 140 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 141 CppTypeFor<TypeCategory::Real, 16> RTNAME(SumReal16)(const Descriptor &x,
 142     const char *source, int line, int dim, const Descriptor *mask) {
 143   return GetTotalReduction<TypeCategory::Real, 16>(
 144       x, source, line, dim, mask, RealSumAccumulator<long double>{x}, "SUM");
 145 }
 146 #endif
 147
 148 void RTNAME(CppSumComplex4)(CppTypeFor<TypeCategory::Complex, 4> &result,
 149     const Descriptor &x, const char *source, int line, int dim,
 150     const Descriptor *mask) {
 151   result = GetTotalReduction<TypeCategory::Complex, 4>(
 152       x, source, line, dim, mask, ComplexSumAccumulator<double>{x}, "SUM");
 153 }
 154 void RTNAME(CppSumComplex8)(CppTypeFor<TypeCategory::Complex, 8> &result,
 155     const Descriptor &x, const char *source, int line, int dim,
 156     const Descriptor *mask) {
 157   result = GetTotalReduction<TypeCategory::Complex, 8>(
 158       x, source, line, dim, mask, ComplexSumAccumulator<double>{x}, "SUM");
 159 }
 160 #if LDBL_MANT_DIG == 64
 161 void RTNAME(CppSumComplex10)(CppTypeFor<TypeCategory::Complex, 10> &result,
 162     const Descriptor &x, const char *source, int line, int dim,
 163     const Descriptor *mask) {
 164   result = GetTotalReduction<TypeCategory::Complex, 10>(
 165       x, source, line, dim, mask, ComplexSumAccumulator<long double>{x}, "SUM");
 166 }
 167 #elif LDBL_MANT_DIG == 113
 168 void RTNAME(CppSumComplex16)(CppTypeFor<TypeCategory::Complex, 16> &result,
 169     const Descriptor &x, const char *source, int line, int dim,
 170     const Descriptor *mask) {
 171   result = GetTotalReduction<TypeCategory::Complex, 16>(
 172       x, source, line, dim, mask, ComplexSumAccumulator<long double>{x}, "SUM");
 173 }
 174 #endif
 175
 176 void RTNAME(SumDim)(Descriptor &result, const Descriptor &x, int dim,
 177     const char *source, int line, const Descriptor *mask) {
 178   TypedPartialNumericReduction<IntegerSumAccumulator, RealSumAccumulator,
 179       ComplexSumAccumulator>(result, x, dim, source, line, mask, "SUM");
 180 }
 181 } // extern "C"
 182 } // namespace Fortran::runtime