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:
IntegerSumAccumulator(const Descriptor & array)26 explicit IntegerSumAccumulator(const Descriptor &array) : array_{array} {}
Reinitialize()27 void Reinitialize() { sum_ = 0; }
GetResult(A * p,int=-1) const28 template <typename A> void GetResult(A *p, int /*zeroBasedDim*/ = -1) const {
29 *p = static_cast<A>(sum_);
30 }
AccumulateAt(const SubscriptValue at[])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:
RealSumAccumulator(const Descriptor & array)43 explicit RealSumAccumulator(const Descriptor &array) : array_{array} {}
Reinitialize()44 void Reinitialize() { sum_ = correction_ = 0; }
Result() const45 template <typename A> A Result() const { return sum_; }
GetResult(A * p,int=-1) const46 template <typename A> void GetResult(A *p, int /*zeroBasedDim*/ = -1) const {
47 *p = Result<A>();
48 }
Accumulate(A x)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 }
AccumulateAt(const SubscriptValue at[])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:
ComplexSumAccumulator(const Descriptor & array)68 explicit ComplexSumAccumulator(const Descriptor &array) : array_{array} {}
Reinitialize()69 void Reinitialize() {
70 reals_.Reinitialize();
71 imaginaries_.Reinitialize();
72 }
GetResult(A * p,int=-1) const73 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 }
Accumulate(const A & z)78 template <typename A> bool Accumulate(const A &z) {
79 reals_.Accumulate(z.real());
80 imaginaries_.Accumulate(z.imag());
81 return true;
82 }
AccumulateAt(const SubscriptValue at[])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" {
RTNAME(SumInteger1)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 }
RTNAME(SumInteger2)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 }
RTNAME(SumInteger4)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 }
RTNAME(SumInteger8)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__
RTNAME(SumInteger16)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)
RTNAME(SumReal4)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 }
RTNAME(SumReal8)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
RTNAME(SumReal10)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
RTNAME(SumReal16)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
RTNAME(CppSumComplex4)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 }
RTNAME(CppSumComplex8)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
RTNAME(CppSumComplex10)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
RTNAME(CppSumComplex16)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
RTNAME(SumDim)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