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