233 lines · cpp
1//===-- lib/runtime/sum.cpp -------------------------------------*- C++ -*-===//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-exception6//7//===----------------------------------------------------------------------===//8 9// Implements SUM for all required operand types and shapes.10//11// Real and complex SUM reductions attempt to reduce floating-point12// cancellation on intermediate results by using "Kahan summation"13// (basically the same as manual "double-double").14 15#include "flang-rt/runtime/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 22namespace Fortran::runtime {23 24template <typename INTERMEDIATE> class IntegerSumAccumulator {25public:26 explicit RT_API_ATTRS IntegerSumAccumulator(const Descriptor &array)27 : array_{array} {}28 void RT_API_ATTRS Reinitialize() { sum_ = 0; }29 template <typename A>30 RT_API_ATTRS void GetResult(A *p, int /*zeroBasedDim*/ = -1) const {31 *p = static_cast<A>(sum_);32 }33 template <typename A>34 RT_API_ATTRS bool AccumulateAt(const SubscriptValue at[]) {35 sum_ += *array_.Element<A>(at);36 return true;37 }38 39private:40 const Descriptor &array_;41 INTERMEDIATE sum_{0};42};43 44template <typename INTERMEDIATE> class RealSumAccumulator {45public:46 explicit RT_API_ATTRS RealSumAccumulator(const Descriptor &array)47 : array_{array} {}48 void RT_API_ATTRS Reinitialize() { sum_ = correction_ = 0; }49 template <typename A> RT_API_ATTRS A Result() const { return sum_; }50 template <typename A>51 RT_API_ATTRS void GetResult(A *p, int /*zeroBasedDim*/ = -1) const {52 *p = Result<A>();53 }54 template <typename A> RT_API_ATTRS bool Accumulate(A x) {55 // Kahan summation56 auto next{x - correction_};57 auto oldSum{sum_};58 sum_ += next;59 correction_ = (sum_ - oldSum) - next; // algebraically zero60 return true;61 }62 template <typename A>63 RT_API_ATTRS bool AccumulateAt(const SubscriptValue at[]) {64 return Accumulate(*array_.Element<A>(at));65 }66 67private:68 const Descriptor &array_;69 INTERMEDIATE sum_{0.0}, correction_{0.0};70};71 72template <typename PART> class ComplexSumAccumulator {73public:74 explicit RT_API_ATTRS ComplexSumAccumulator(const Descriptor &array)75 : array_{array} {}76 void RT_API_ATTRS Reinitialize() {77 reals_.Reinitialize();78 imaginaries_.Reinitialize();79 }80 template <typename A>81 RT_API_ATTRS void GetResult(A *p, int /*zeroBasedDim*/ = -1) const {82 using ResultPart = typename A::value_type;83 *p = {reals_.template Result<ResultPart>(),84 imaginaries_.template Result<ResultPart>()};85 }86 template <typename A> RT_API_ATTRS bool Accumulate(const A &z) {87 reals_.Accumulate(z.real());88 imaginaries_.Accumulate(z.imag());89 return true;90 }91 template <typename A>92 RT_API_ATTRS bool AccumulateAt(const SubscriptValue at[]) {93 return Accumulate(*array_.Element<A>(at));94 }95 96private:97 const Descriptor &array_;98 RealSumAccumulator<PART> reals_{array_}, imaginaries_{array_};99};100 101extern "C" {102RT_EXT_API_GROUP_BEGIN103 104CppTypeFor<TypeCategory::Integer, 1> RTDEF(SumInteger1)(const Descriptor &x,105 const char *source, int line, int dim, const Descriptor *mask) {106 return GetTotalReduction<TypeCategory::Integer, 1>(x, source, line, dim, mask,107 IntegerSumAccumulator<CppTypeFor<TypeCategory::Integer, 4>>{x}, "SUM");108}109CppTypeFor<TypeCategory::Integer, 2> RTDEF(SumInteger2)(const Descriptor &x,110 const char *source, int line, int dim, const Descriptor *mask) {111 return GetTotalReduction<TypeCategory::Integer, 2>(x, source, line, dim, mask,112 IntegerSumAccumulator<CppTypeFor<TypeCategory::Integer, 4>>{x}, "SUM");113}114CppTypeFor<TypeCategory::Integer, 4> RTDEF(SumInteger4)(const Descriptor &x,115 const char *source, int line, int dim, const Descriptor *mask) {116 return GetTotalReduction<TypeCategory::Integer, 4>(x, source, line, dim, mask,117 IntegerSumAccumulator<CppTypeFor<TypeCategory::Integer, 4>>{x}, "SUM");118}119CppTypeFor<TypeCategory::Integer, 8> RTDEF(SumInteger8)(const Descriptor &x,120 const char *source, int line, int dim, const Descriptor *mask) {121 return GetTotalReduction<TypeCategory::Integer, 8>(x, source, line, dim, mask,122 IntegerSumAccumulator<CppTypeFor<TypeCategory::Integer, 8>>{x}, "SUM");123}124#ifdef __SIZEOF_INT128__125CppTypeFor<TypeCategory::Integer, 16> RTDEF(SumInteger16)(const Descriptor &x,126 const char *source, int line, int dim, const Descriptor *mask) {127 return GetTotalReduction<TypeCategory::Integer, 16>(x, source, line, dim,128 mask, IntegerSumAccumulator<CppTypeFor<TypeCategory::Integer, 16>>{x},129 "SUM");130}131#endif132 133CppTypeFor<TypeCategory::Unsigned, 1> RTDEF(SumUnsigned1)(const Descriptor &x,134 const char *source, int line, int dim, const Descriptor *mask) {135 return GetTotalReduction<TypeCategory::Unsigned, 1>(x, source, line, dim,136 mask, IntegerSumAccumulator<CppTypeFor<TypeCategory::Unsigned, 4>>{x},137 "SUM");138}139CppTypeFor<TypeCategory::Unsigned, 2> RTDEF(SumUnsigned2)(const Descriptor &x,140 const char *source, int line, int dim, const Descriptor *mask) {141 return GetTotalReduction<TypeCategory::Unsigned, 2>(x, source, line, dim,142 mask, IntegerSumAccumulator<CppTypeFor<TypeCategory::Unsigned, 4>>{x},143 "SUM");144}145CppTypeFor<TypeCategory::Unsigned, 4> RTDEF(SumUnsigned4)(const Descriptor &x,146 const char *source, int line, int dim, const Descriptor *mask) {147 return GetTotalReduction<TypeCategory::Unsigned, 4>(x, source, line, dim,148 mask, IntegerSumAccumulator<CppTypeFor<TypeCategory::Unsigned, 4>>{x},149 "SUM");150}151CppTypeFor<TypeCategory::Unsigned, 8> RTDEF(SumUnsigned8)(const Descriptor &x,152 const char *source, int line, int dim, const Descriptor *mask) {153 return GetTotalReduction<TypeCategory::Unsigned, 8>(x, source, line, dim,154 mask, IntegerSumAccumulator<CppTypeFor<TypeCategory::Unsigned, 8>>{x},155 "SUM");156}157#ifdef __SIZEOF_INT128__158CppTypeFor<TypeCategory::Unsigned, 16> RTDEF(SumUnsigned16)(const Descriptor &x,159 const char *source, int line, int dim, const Descriptor *mask) {160 return GetTotalReduction<TypeCategory::Unsigned, 16>(x, source, line, dim,161 mask, IntegerSumAccumulator<CppTypeFor<TypeCategory::Unsigned, 16>>{x},162 "SUM");163}164#endif165 166// TODO: real/complex(2 & 3)167CppTypeFor<TypeCategory::Real, 4> RTDEF(SumReal4)(const Descriptor &x,168 const char *source, int line, int dim, const Descriptor *mask) {169 return GetTotalReduction<TypeCategory::Real, 4>(170 x, source, line, dim, mask, RealSumAccumulator<float>{x}, "SUM");171}172CppTypeFor<TypeCategory::Real, 8> RTDEF(SumReal8)(const Descriptor &x,173 const char *source, int line, int dim, const Descriptor *mask) {174 return GetTotalReduction<TypeCategory::Real, 8>(175 x, source, line, dim, mask, RealSumAccumulator<double>{x}, "SUM");176}177#if HAS_FLOAT80178CppTypeFor<TypeCategory::Real, 10> RTDEF(SumReal10)(const Descriptor &x,179 const char *source, int line, int dim, const Descriptor *mask) {180 return GetTotalReduction<TypeCategory::Real, 10>(x, source, line, dim, mask,181 RealSumAccumulator<CppTypeFor<TypeCategory::Real, 10>>{x}, "SUM");182}183#endif184#if HAS_LDBL128 || HAS_FLOAT128185CppTypeFor<TypeCategory::Real, 16> RTDEF(SumReal16)(const Descriptor &x,186 const char *source, int line, int dim, const Descriptor *mask) {187 return GetTotalReduction<TypeCategory::Real, 16>(x, source, line, dim, mask,188 RealSumAccumulator<CppTypeFor<TypeCategory::Real, 16>>{x}, "SUM");189}190#endif191 192void RTDEF(CppSumComplex4)(CppTypeFor<TypeCategory::Complex, 4> &result,193 const Descriptor &x, const char *source, int line, int dim,194 const Descriptor *mask) {195 result = GetTotalReduction<TypeCategory::Complex, 4>(196 x, source, line, dim, mask, ComplexSumAccumulator<float>{x}, "SUM");197}198void RTDEF(CppSumComplex8)(CppTypeFor<TypeCategory::Complex, 8> &result,199 const Descriptor &x, const char *source, int line, int dim,200 const Descriptor *mask) {201 result = GetTotalReduction<TypeCategory::Complex, 8>(202 x, source, line, dim, mask, ComplexSumAccumulator<double>{x}, "SUM");203}204#if HAS_FLOAT80205void RTDEF(CppSumComplex10)(CppTypeFor<TypeCategory::Complex, 10> &result,206 const Descriptor &x, const char *source, int line, int dim,207 const Descriptor *mask) {208 result =209 GetTotalReduction<TypeCategory::Complex, 10>(x, source, line, dim, mask,210 ComplexSumAccumulator<CppTypeFor<TypeCategory::Real, 10>>{x}, "SUM");211}212#endif213#if HAS_LDBL128 || HAS_FLOAT128214void RTDEF(CppSumComplex16)(CppTypeFor<TypeCategory::Complex, 16> &result,215 const Descriptor &x, const char *source, int line, int dim,216 const Descriptor *mask) {217 result =218 GetTotalReduction<TypeCategory::Complex, 16>(x, source, line, dim, mask,219 ComplexSumAccumulator<CppTypeFor<TypeCategory::Real, 16>>{x}, "SUM");220}221#endif222 223void RTDEF(SumDim)(Descriptor &result, const Descriptor &x, int dim,224 const char *source, int line, const Descriptor *mask) {225 TypedPartialNumericReduction<IntegerSumAccumulator, RealSumAccumulator,226 ComplexSumAccumulator, /*MIN_REAL_KIND=*/4>(227 result, x, dim, source, line, mask, "SUM");228}229 230RT_EXT_API_GROUP_END231} // extern "C"232} // namespace Fortran::runtime233