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 "reduction.h"
17 #include "flang/Common/long-double.h"
18 #include <cinttypes>
19 #include <complex>
20
21 namespace Fortran::runtime {
22
23 template <typename INTERMEDIATE> class IntegerSumAccumulator {
24 public:
IntegerSumAccumulator(const Descriptor & array)25 explicit IntegerSumAccumulator(const Descriptor &array) : array_{array} {}
Reinitialize()26 void Reinitialize() { sum_ = 0; }
GetResult(A * p,int=-1) const27 template <typename A> void GetResult(A *p, int /*zeroBasedDim*/ = -1) const {
28 *p = static_cast<A>(sum_);
29 }
AccumulateAt(const SubscriptValue at[])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:
RealSumAccumulator(const Descriptor & array)42 explicit RealSumAccumulator(const Descriptor &array) : array_{array} {}
Reinitialize()43 void Reinitialize() { sum_ = correction_ = 0; }
Result() const44 template <typename A> A Result() const { return sum_; }
GetResult(A * p,int=-1) const45 template <typename A> void GetResult(A *p, int /*zeroBasedDim*/ = -1) const {
46 *p = Result<A>();
47 }
Accumulate(A x)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 }
AccumulateAt(const SubscriptValue at[])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:
ComplexSumAccumulator(const Descriptor & array)67 explicit ComplexSumAccumulator(const Descriptor &array) : array_{array} {}
Reinitialize()68 void Reinitialize() {
69 reals_.Reinitialize();
70 imaginaries_.Reinitialize();
71 }
GetResult(A * p,int=-1) const72 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 }
Accumulate(const A & z)77 template <typename A> bool Accumulate(const A &z) {
78 reals_.Accumulate(z.real());
79 imaginaries_.Accumulate(z.imag());
80 return true;
81 }
AccumulateAt(const SubscriptValue at[])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" {
RTNAME(SumInteger1)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 }
RTNAME(SumInteger2)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 }
RTNAME(SumInteger4)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 }
RTNAME(SumInteger8)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__
RTNAME(SumInteger16)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)
RTNAME(SumReal4)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 }
RTNAME(SumReal8)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
RTNAME(SumReal10)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
RTNAME(SumReal16)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
RTNAME(CppSumComplex4)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 }
RTNAME(CppSumComplex8)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
RTNAME(CppSumComplex10)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
RTNAME(CppSumComplex16)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
RTNAME(SumDim)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