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/Runtime/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 183