433 lines · c
1//===-- lib/Evaluate/fold-reduction.h -------------------------------------===//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#ifndef FORTRAN_EVALUATE_FOLD_REDUCTION_H_10#define FORTRAN_EVALUATE_FOLD_REDUCTION_H_11 12#include "fold-implementation.h"13 14namespace Fortran::evaluate {15 16// DOT_PRODUCT17template <typename T>18static Expr<T> FoldDotProduct(19 FoldingContext &context, FunctionRef<T> &&funcRef) {20 using Element = typename Constant<T>::Element;21 auto args{funcRef.arguments()};22 CHECK(args.size() == 2);23 Folder<T> folder{context};24 Constant<T> *va{folder.Folding(args[0])};25 Constant<T> *vb{folder.Folding(args[1])};26 if (va && vb) {27 CHECK(va->Rank() == 1 && vb->Rank() == 1);28 if (va->size() != vb->size()) {29 context.messages().Say(30 "Vector arguments to DOT_PRODUCT have distinct extents %zd and %zd"_err_en_US,31 va->size(), vb->size());32 return MakeInvalidIntrinsic(std::move(funcRef));33 }34 Element sum{};35 bool overflow{false};36 if constexpr (T::category == TypeCategory::Complex) {37 std::vector<Element> conjugates;38 for (const Element &x : va->values()) {39 conjugates.emplace_back(x.CONJG());40 }41 Constant<T> conjgA{42 std::move(conjugates), ConstantSubscripts{va->shape()}};43 Expr<T> products{Fold(44 context, Expr<T>{std::move(conjgA)} * Expr<T>{Constant<T>{*vb}})};45 Constant<T> &cProducts{DEREF(UnwrapConstantValue<T>(products))};46 [[maybe_unused]] Element correction{};47 const auto &rounding{context.targetCharacteristics().roundingMode()};48 for (const Element &x : cProducts.values()) {49 if constexpr (useKahanSummation) {50 auto next{x.Subtract(correction, rounding)};51 overflow |= next.flags.test(RealFlag::Overflow);52 auto added{sum.Add(next.value, rounding)};53 overflow |= added.flags.test(RealFlag::Overflow);54 correction = added.value.Subtract(sum, rounding)55 .value.Subtract(next.value, rounding)56 .value;57 sum = std::move(added.value);58 } else {59 auto added{sum.Add(x, rounding)};60 overflow |= added.flags.test(RealFlag::Overflow);61 sum = std::move(added.value);62 }63 }64 } else if constexpr (T::category == TypeCategory::Logical) {65 Expr<T> conjunctions{Fold(context,66 Expr<T>{LogicalOperation<T::kind>{LogicalOperator::And,67 Expr<T>{Constant<T>{*va}}, Expr<T>{Constant<T>{*vb}}}})};68 Constant<T> &cConjunctions{DEREF(UnwrapConstantValue<T>(conjunctions))};69 for (const Element &x : cConjunctions.values()) {70 if (x.IsTrue()) {71 sum = Element{true};72 break;73 }74 }75 } else if constexpr (T::category == TypeCategory::Integer) {76 Expr<T> products{77 Fold(context, Expr<T>{Constant<T>{*va}} * Expr<T>{Constant<T>{*vb}})};78 Constant<T> &cProducts{DEREF(UnwrapConstantValue<T>(products))};79 for (const Element &x : cProducts.values()) {80 auto next{sum.AddSigned(x)};81 overflow |= next.overflow;82 sum = std::move(next.value);83 }84 } else if constexpr (T::category == TypeCategory::Unsigned) {85 Expr<T> products{86 Fold(context, Expr<T>{Constant<T>{*va}} * Expr<T>{Constant<T>{*vb}})};87 Constant<T> &cProducts{DEREF(UnwrapConstantValue<T>(products))};88 for (const Element &x : cProducts.values()) {89 sum = sum.AddUnsigned(x).value;90 }91 } else {92 static_assert(T::category == TypeCategory::Real);93 Expr<T> products{94 Fold(context, Expr<T>{Constant<T>{*va}} * Expr<T>{Constant<T>{*vb}})};95 Constant<T> &cProducts{DEREF(UnwrapConstantValue<T>(products))};96 [[maybe_unused]] Element correction{};97 const auto &rounding{context.targetCharacteristics().roundingMode()};98 for (const Element &x : cProducts.values()) {99 if constexpr (useKahanSummation) {100 auto next{x.Subtract(correction, rounding)};101 overflow |= next.flags.test(RealFlag::Overflow);102 auto added{sum.Add(next.value, rounding)};103 overflow |= added.flags.test(RealFlag::Overflow);104 correction = added.value.Subtract(sum, rounding)105 .value.Subtract(next.value, rounding)106 .value;107 sum = std::move(added.value);108 } else {109 auto added{sum.Add(x, rounding)};110 overflow |= added.flags.test(RealFlag::Overflow);111 sum = std::move(added.value);112 }113 }114 }115 if (overflow) {116 context.Warn(common::UsageWarning::FoldingException,117 "DOT_PRODUCT of %s data overflowed during computation"_warn_en_US,118 T::AsFortran());119 }120 return Expr<T>{Constant<T>{std::move(sum)}};121 }122 return Expr<T>{std::move(funcRef)};123}124 125// Fold and validate a DIM= argument. Returns false on error.126bool CheckReductionDIM(std::optional<int> &dim, FoldingContext &,127 ActualArguments &, std::optional<int> dimIndex, int rank);128 129// Fold and validate a MASK= argument. Return null on error, absent MASK=, or130// non-constant MASK=.131Constant<LogicalResult> *GetReductionMASK(132 std::optional<ActualArgument> &maskArg, const ConstantSubscripts &shape,133 FoldingContext &);134 135// Common preprocessing for reduction transformational intrinsic function136// folding. If the intrinsic can have DIM= &/or MASK= arguments, extract137// and check them. If a MASK= is present, apply it to the array data and138// substitute replacement values for elements corresponding to .FALSE. in139// the mask. If the result is present, the intrinsic call can be folded.140template <typename T> struct ArrayAndMask {141 Constant<T> array;142 Constant<LogicalResult> mask;143};144template <typename T>145static std::optional<ArrayAndMask<T>> ProcessReductionArgs(146 FoldingContext &context, ActualArguments &arg, std::optional<int> &dim,147 int arrayIndex, std::optional<int> dimIndex = std::nullopt,148 std::optional<int> maskIndex = std::nullopt) {149 if (arg.empty()) {150 return std::nullopt;151 }152 Constant<T> *folded{Folder<T>{context}.Folding(arg[arrayIndex])};153 if (!folded || folded->Rank() < 1) {154 return std::nullopt;155 }156 if (!CheckReductionDIM(dim, context, arg, dimIndex, folded->Rank())) {157 return std::nullopt;158 }159 std::size_t n{folded->size()};160 std::vector<Scalar<LogicalResult>> maskElement;161 if (maskIndex && static_cast<std::size_t>(*maskIndex) < arg.size() &&162 arg[*maskIndex]) {163 if (const Constant<LogicalResult> *origMask{164 GetReductionMASK(arg[*maskIndex], folded->shape(), context)}) {165 if (auto scalarMask{origMask->GetScalarValue()}) {166 maskElement =167 std::vector<Scalar<LogicalResult>>(n, scalarMask->IsTrue());168 } else {169 maskElement = origMask->values();170 }171 } else {172 return std::nullopt;173 }174 } else {175 maskElement = std::vector<Scalar<LogicalResult>>(n, true);176 }177 return ArrayAndMask<T>{Constant<T>(*folded),178 Constant<LogicalResult>{179 std::move(maskElement), ConstantSubscripts{folded->shape()}}};180}181 182// Generalized reduction to an array of one dimension fewer (w/ DIM=)183// or to a scalar (w/o DIM=). The ACCUMULATOR type must define184// operator()(Scalar<T> &, const ConstantSubscripts &, bool first)185// and Done(Scalar<T> &).186template <typename T, typename ACCUMULATOR, typename ARRAY>187static Constant<T> DoReduction(const Constant<ARRAY> &array,188 const Constant<LogicalResult> &mask, std::optional<int> &dim,189 const Scalar<T> &identity, ACCUMULATOR &accumulator) {190 ConstantSubscripts at{array.lbounds()};191 ConstantSubscripts maskAt{mask.lbounds()};192 std::vector<typename Constant<T>::Element> elements;193 ConstantSubscripts resultShape; // empty -> scalar194 if (dim) { // DIM= is present, so result is an array195 resultShape = array.shape();196 resultShape.erase(resultShape.begin() + (*dim - 1));197 ConstantSubscript dimExtent{array.shape().at(*dim - 1)};198 CHECK(dimExtent == mask.shape().at(*dim - 1));199 ConstantSubscript &dimAt{at[*dim - 1]};200 ConstantSubscript dimLbound{dimAt};201 ConstantSubscript &maskDimAt{maskAt[*dim - 1]};202 ConstantSubscript maskDimLbound{maskDimAt};203 for (auto n{GetSize(resultShape)}; n-- > 0;204 array.IncrementSubscripts(at), mask.IncrementSubscripts(maskAt)) {205 elements.push_back(identity);206 if (dimExtent > 0) {207 dimAt = dimLbound;208 maskDimAt = maskDimLbound;209 bool firstUnmasked{true};210 for (ConstantSubscript j{0}; j < dimExtent; ++j, ++dimAt, ++maskDimAt) {211 if (mask.At(maskAt).IsTrue()) {212 accumulator(elements.back(), at, firstUnmasked);213 firstUnmasked = false;214 }215 }216 --dimAt, --maskDimAt;217 }218 accumulator.Done(elements.back());219 }220 } else { // no DIM=, result is scalar221 elements.push_back(identity);222 bool firstUnmasked{true};223 for (auto n{array.size()}; n-- > 0;224 array.IncrementSubscripts(at), mask.IncrementSubscripts(maskAt)) {225 if (mask.At(maskAt).IsTrue()) {226 accumulator(elements.back(), at, firstUnmasked);227 firstUnmasked = false;228 }229 }230 accumulator.Done(elements.back());231 }232 if constexpr (T::category == TypeCategory::Character) {233 return {static_cast<ConstantSubscript>(identity.size()),234 std::move(elements), std::move(resultShape)};235 } else {236 return {std::move(elements), std::move(resultShape)};237 }238}239 240// MAXVAL & MINVAL241template <typename T, bool ABS = false> class MaxvalMinvalAccumulator {242public:243 MaxvalMinvalAccumulator(244 RelationalOperator opr, FoldingContext &context, const Constant<T> &array)245 : opr_{opr}, context_{context}, array_{array} {};246 void operator()(Scalar<T> &element, const ConstantSubscripts &at,247 [[maybe_unused]] bool firstUnmasked) const {248 auto aAt{array_.At(at)};249 if constexpr (ABS) {250 aAt = aAt.ABS();251 }252 if constexpr (T::category == TypeCategory::Real) {253 if (firstUnmasked || element.IsNotANumber()) {254 // Return NaN if and only if all unmasked elements are NaNs and255 // at least one unmasked element is visible.256 element = aAt;257 return;258 }259 }260 Expr<LogicalResult> test{PackageRelation(261 opr_, Expr<T>{Constant<T>{aAt}}, Expr<T>{Constant<T>{element}})};262 auto folded{GetScalarConstantValue<LogicalResult>(263 test.Rewrite(context_, std::move(test)))};264 CHECK(folded.has_value());265 if (folded->IsTrue()) {266 element = aAt;267 }268 }269 void Done(Scalar<T> &) const {}270 271private:272 RelationalOperator opr_;273 FoldingContext &context_;274 const Constant<T> &array_;275};276 277template <typename T>278static Expr<T> FoldMaxvalMinval(FoldingContext &context, FunctionRef<T> &&ref,279 RelationalOperator opr, const Scalar<T> &identity) {280 static_assert(T::category == TypeCategory::Integer ||281 T::category == TypeCategory::Unsigned ||282 T::category == TypeCategory::Real ||283 T::category == TypeCategory::Character);284 std::optional<int> dim;285 if (std::optional<ArrayAndMask<T>> arrayAndMask{286 ProcessReductionArgs<T>(context, ref.arguments(), dim,287 /*ARRAY=*/0, /*DIM=*/1, /*MASK=*/2)}) {288 MaxvalMinvalAccumulator<T> accumulator{opr, context, arrayAndMask->array};289 return Expr<T>{DoReduction<T>(290 arrayAndMask->array, arrayAndMask->mask, dim, identity, accumulator)};291 }292 return Expr<T>{std::move(ref)};293}294 295// PRODUCT296template <typename T> class ProductAccumulator {297public:298 ProductAccumulator(const Constant<T> &array) : array_{array} {}299 void operator()(300 Scalar<T> &element, const ConstantSubscripts &at, bool /*first*/) {301 if constexpr (T::category == TypeCategory::Integer) {302 auto prod{element.MultiplySigned(array_.At(at))};303 overflow_ |= prod.SignedMultiplicationOverflowed();304 element = prod.lower;305 } else if constexpr (T::category == TypeCategory::Unsigned) {306 element = element.MultiplyUnsigned(array_.At(at)).lower;307 } else { // Real & Complex308 auto prod{element.Multiply(array_.At(at))};309 overflow_ |= prod.flags.test(RealFlag::Overflow);310 element = prod.value;311 }312 }313 bool overflow() const { return overflow_; }314 void Done(Scalar<T> &) const {}315 316private:317 const Constant<T> &array_;318 bool overflow_{false};319};320 321template <typename T>322static Expr<T> FoldProduct(323 FoldingContext &context, FunctionRef<T> &&ref, Scalar<T> identity) {324 static_assert(T::category == TypeCategory::Integer ||325 T::category == TypeCategory::Unsigned ||326 T::category == TypeCategory::Real ||327 T::category == TypeCategory::Complex);328 std::optional<int> dim;329 if (std::optional<ArrayAndMask<T>> arrayAndMask{330 ProcessReductionArgs<T>(context, ref.arguments(), dim,331 /*ARRAY=*/0, /*DIM=*/1, /*MASK=*/2)}) {332 ProductAccumulator accumulator{arrayAndMask->array};333 auto result{Expr<T>{DoReduction<T>(334 arrayAndMask->array, arrayAndMask->mask, dim, identity, accumulator)}};335 if (accumulator.overflow()) {336 context.Warn(common::UsageWarning::FoldingException,337 "PRODUCT() of %s data overflowed"_warn_en_US, T::AsFortran());338 }339 return result;340 }341 return Expr<T>{std::move(ref)};342}343 344// SUM345template <typename T> class SumAccumulator {346 using Element = typename Constant<T>::Element;347 348public:349 SumAccumulator(const Constant<T> &array, Rounding rounding)350 : array_{array}, rounding_{rounding} {}351 void operator()(352 Element &element, const ConstantSubscripts &at, bool /*first*/) {353 if constexpr (T::category == TypeCategory::Integer) {354 auto sum{element.AddSigned(array_.At(at))};355 overflow_ |= sum.overflow;356 element = sum.value;357 } else if constexpr (T::category == TypeCategory::Unsigned) {358 element = element.AddUnsigned(array_.At(at)).value;359 } else { // Real & Complex: use Kahan summation360 auto next{array_.At(at).Subtract(correction_, rounding_)};361 overflow_ |= next.flags.test(RealFlag::Overflow);362 auto sum{element.Add(next.value, rounding_)};363 overflow_ |= sum.flags.test(RealFlag::Overflow);364 // correction = (sum - element) - next; algebraically zero365 correction_ = sum.value.Subtract(element, rounding_)366 .value.Subtract(next.value, rounding_)367 .value;368 element = sum.value;369 }370 }371 bool overflow() const { return overflow_; }372 void Done([[maybe_unused]] Element &element) {373 if constexpr (T::category != TypeCategory::Integer &&374 T::category != TypeCategory::Unsigned) {375 auto corrected{element.Add(correction_, rounding_)};376 overflow_ |= corrected.flags.test(RealFlag::Overflow);377 correction_ = Scalar<T>{};378 element = corrected.value;379 }380 }381 382private:383 const Constant<T> &array_;384 Rounding rounding_;385 bool overflow_{false};386 Element correction_{};387};388 389template <typename T>390static Expr<T> FoldSum(FoldingContext &context, FunctionRef<T> &&ref) {391 static_assert(T::category == TypeCategory::Integer ||392 T::category == TypeCategory::Unsigned ||393 T::category == TypeCategory::Real ||394 T::category == TypeCategory::Complex);395 using Element = typename Constant<T>::Element;396 std::optional<int> dim;397 Element identity{};398 if (std::optional<ArrayAndMask<T>> arrayAndMask{399 ProcessReductionArgs<T>(context, ref.arguments(), dim,400 /*ARRAY=*/0, /*DIM=*/1, /*MASK=*/2)}) {401 SumAccumulator accumulator{402 arrayAndMask->array, context.targetCharacteristics().roundingMode()};403 auto result{Expr<T>{DoReduction<T>(404 arrayAndMask->array, arrayAndMask->mask, dim, identity, accumulator)}};405 if (accumulator.overflow()) {406 context.Warn(common::UsageWarning::FoldingException,407 "SUM() of %s data overflowed"_warn_en_US, T::AsFortran());408 }409 return result;410 }411 return Expr<T>{std::move(ref)};412}413 414// Utility for IALL, IANY, IPARITY, ALL, ANY, & PARITY415template <typename T> class OperationAccumulator {416public:417 OperationAccumulator(const Constant<T> &array,418 Scalar<T> (Scalar<T>::*operation)(const Scalar<T> &) const)419 : array_{array}, operation_{operation} {}420 void operator()(421 Scalar<T> &element, const ConstantSubscripts &at, bool /*first*/) {422 element = (element.*operation_)(array_.At(at));423 }424 void Done(Scalar<T> &) const {}425 426private:427 const Constant<T> &array_;428 Scalar<T> (Scalar<T>::*operation_)(const Scalar<T> &) const;429};430 431} // namespace Fortran::evaluate432#endif // FORTRAN_EVALUATE_FOLD_REDUCTION_H_433