[flang] Fold more reduction intrinsic function calls

[klausler]

Summary

Refactor the recently-implemented MAXVAL/MINVAL folding so
that the parts that can be used to implement other reduction
transformational intrinsic function folding are exposed.

Use them to implement folding of IALL, IANY, IPARITY,
SUM. and PRODUCT. Replace the folding of ALL & ANY to
use the new infrastructure and become able to handle DIM=
arguments.

[jeanPerier]

Looks good

flang/include/flang/Evaluate/real.h

public:
   static constexpr int bits{Word::bits};
   static_assert(bits >= Details::bits);
   using Fraction = Integer<binaryPrecision>; // all bits made explicit
 
   template <typename W, int P> friend class Real;
 
   constexpr Real() {} // +0.0
   constexpr Real(const Real &) = default;
+  constexpr Real(Real &&) = default;
   constexpr Real(const Word &bits) : word_{bits} {}
   constexpr Real &operator=(const Real &) = default;
   constexpr Real &operator=(Real &&) = default;
 
   constexpr bool operator==(const Real &that) const {
     return word_ == that.word_;
   }
 

flang/lib/Evaluate/fold-character.cpp

if (name == "achar" || name == "char") {
     }
   } else if (name == "trim") { // not elemental
     if (auto scalar{
             GetScalarConstantArguments<T>(context, funcRef.arguments())}) {
       return Expr<T>{Constant<T>{
           CharacterUtils<KIND>::TRIM(std::get<Scalar<T>>(*scalar))}};
     }
   }
-  // TODO: cshift, eoshift, maxloc, minloc, pack, reduce,
-  // spread, transfer, transpose, unpack
+  // TODO: cshift, eoshift, maxloc, minloc, pack, spread, transfer,
+  // transpose, unpack
   return Expr<T>{std::move(funcRef)};
 }
 
 template <int KIND>
 Expr<Type<TypeCategory::Character, KIND>> FoldOperation(
     FoldingContext &context, Concat<KIND> &&x) {
   if (auto array{ApplyElementwise(context, x)}) {
     return *array;

flang/lib/Evaluate/fold-complex.cpp

 //===-- lib/Evaluate/fold-complex.cpp -------------------------------------===//
 //
 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
 // See https://llvm.org/LICENSE.txt for license information.
 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
 //
 //===----------------------------------------------------------------------===//
 
 #include "fold-implementation.h"
+#include "fold-reduction.h"
 
 namespace Fortran::evaluate {
 
 template <int KIND>
 Expr<Type<TypeCategory::Complex, KIND>> FoldIntrinsicFunction(
     FoldingContext &context,
     FunctionRef<Type<TypeCategory::Complex, KIND>> &&funcRef) {
   using T = Type<TypeCategory::Complex, KIND>;
+  using Part = typename T::Part;
   ActualArguments &args{funcRef.arguments()};
   auto *intrinsic{std::get_if<SpecificIntrinsic>(&funcRef.proc().u)};
   CHECK(intrinsic);
   std::string name{intrinsic->name};
   if (name == "acos" || name == "acosh" || name == "asin" || name == "asinh" ||
       name == "atan" || name == "atanh" || name == "cos" || name == "cosh" ||
       name == "exp" || name == "log" || name == "sin" || name == "sinh" ||
       name == "sqrt" || name == "tan" || name == "tanh") {
return FoldElementalIntrinsic<T, T>(
         context, std::move(funcRef), &Scalar<T>::CONJG);
   } else if (name == "cmplx") {
     if (args.size() > 0 && args[0].has_value()) {
       if (auto *x{UnwrapExpr<Expr<SomeComplex>>(args[0])}) {
         // CMPLX(X [, KIND]) with complex X
         return Fold(context, ConvertToType<T>(std::move(*x)));
       } else {
         // CMPLX(X [, Y [, KIND]]) with non-complex X
-        using Part = typename T::Part;
         Expr<SomeType> re{std::move(*args[0].value().UnwrapExpr())};
         Expr<SomeType> im{args.size() >= 2 && args[1].has_value()
                 ? std::move(*args[1]->UnwrapExpr())
                 : AsGenericExpr(Constant<Part>{Scalar<Part>{}})};
         return Fold(context,
             Expr<T>{
                 ComplexConstructor<KIND>{ToReal<KIND>(context, std::move(re)),
                     ToReal<KIND>(context, std::move(im))}});
       }
     }
   } else if (name == "merge") {
     return FoldMerge<T>(context, std::move(funcRef));
+  } else if (name == "product") {
+    auto one{Scalar<Part>::FromInteger(value::Integer<8>{1}).value};
+    return FoldProduct<T>(context, std::move(funcRef), Scalar<T>{one});
+  } else if (name == "sum") {
+    return FoldSum<T>(context, std::move(funcRef));
   }
-  // TODO: cshift, dot_product, eoshift, matmul, pack, product,
-  // reduce, spread, sum, transfer, transpose, unpack
+  // TODO: cshift, dot_product, eoshift, matmul, pack, spread, transfer,
+  // transpose, unpack
   return Expr<T>{std::move(funcRef)};
 }
 
 template <int KIND>
 Expr<Type<TypeCategory::Complex, KIND>> FoldOperation(
     FoldingContext &context, ComplexConstructor<KIND> &&x) {
   if (auto array{ApplyElementwise(context, x)}) {
     return *array;

flang/lib/Evaluate/fold-integer.cpp

if (int rank{array->Rank()}; rank > 0) {
           }
         }
       }
     }
   }
   return Expr<T>{std::move(funcRef)};
 }
 
+// for IALL, IANY, & IPARITY
+template <typename T>
+static Expr<T> FoldBitReduction(FoldingContext &context, FunctionRef<T> &&ref,
+    Scalar<T> (Scalar<T>::*operation)(const Scalar<T> &) const,
+    Scalar<T> identity) {
+  static_assert(T::category == TypeCategory::Integer);
+  using Element = Scalar<T>;
+  std::optional<ConstantSubscript> dim;
+  if (std::optional<Constant<T>> array{
+          ProcessReductionArgs<T>(context, ref.arguments(), dim, identity,
+              /*ARRAY=*/0, /*DIM=*/1, /*MASK=*/2)}) {
+    auto accumulator{[&](Element &element, const ConstantSubscripts &at) {
+      element = (element.*operation)(array->At(at));
+    }};
+    return Expr<T>{DoReduction(*array, dim, identity, accumulator)};
+  }
+  return Expr<T>{std::move(ref)};
+}
+
 template <int KIND>
 Expr<Type<TypeCategory::Integer, KIND>> FoldIntrinsicFunction(
     FoldingContext &context,
     FunctionRef<Type<TypeCategory::Integer, KIND>> &&funcRef) {
   using T = Type<TypeCategory::Integer, KIND>;
   using Int4 = Type<TypeCategory::Integer, 4>;
   ActualArguments &args{funcRef.arguments()};
   auto *intrinsic{std::get_if<SpecificIntrinsic>(&funcRef.proc().u)};
if (name == "iand") { // done in fptr declaration
       fptr = &Scalar<T>::IOR;
     } else if (name == "ieor") {
       fptr = &Scalar<T>::IEOR;
     } else {
       common::die("missing case to fold intrinsic function %s", name.c_str());
     }
     return FoldElementalIntrinsic<T, T, T>(
         context, std::move(funcRef), ScalarFunc<T, T, T>(fptr));
+  } else if (name == "iall") {
+    return FoldBitReduction(
+        context, std::move(funcRef), &Scalar<T>::IAND, Scalar<T>{}.NOT());
+  } else if (name == "iany") {
+    return FoldBitReduction(
+        context, std::move(funcRef), &Scalar<T>::IOR, Scalar<T>{});
   } else if (name == "ibclr" || name == "ibset" || name == "ishft" ||
       name == "shifta" || name == "shiftr" || name == "shiftl") {
     // Second argument can be of any kind. However, it must be smaller or
     // equal than BIT_SIZE. It can be converted to Int4 to simplify.
     auto fptr{&Scalar<T>::IBCLR};
     if (name == "ibclr") { // done in fprt definition
     } else if (name == "ibset") {
       fptr = &Scalar<T>::IBSET;
Expr<Type<TypeCategory::Integer, KIND>> FoldIntrinsicFunction(
   } else if (name == "int_ptr_kind") {
     return Expr<T>{8};
   } else if (name == "kind") {
     if constexpr (common::HasMember<T, IntegerTypes>) {
       return Expr<T>{args[0].value().GetType()->kind()};
     } else {
       DIE("kind() result not integral");
     }
+  } else if (name == "iparity") {
+    return FoldBitReduction(
+        context, std::move(funcRef), &Scalar<T>::IEOR, Scalar<T>{});
   } else if (name == "lbound") {
     return LBOUND(context, std::move(funcRef));
   } else if (name == "leadz" || name == "trailz" || name == "poppar" ||
       name == "popcnt") {
     if (auto *sn{UnwrapExpr<Expr<SomeInteger>>(args[0])}) {
       return std::visit(
           [&funcRef, &context, &name](const auto &n) -> Expr<T> {
             using TI = typename std::decay_t<decltype(n)>::Result;
if (const auto *cx{UnwrapExpr<Expr<SomeReal>>(args[0])}) {
           cx->u)};
     } else if (const auto *cx{UnwrapExpr<Expr<SomeComplex>>(args[0])}) {
       return Expr<T>{std::visit(
           [](const auto &kx) {
             return Scalar<typename ResultType<decltype(kx)>::Part>::PRECISION;
           },
           cx->u)};
     }
+  } else if (name == "product") {
+    return FoldProduct<T>(context, std::move(funcRef), Scalar<T>{1});
   } else if (name == "radix") {
     return Expr<T>{2};
   } else if (name == "range") {
     if (const auto *cx{UnwrapExpr<Expr<SomeInteger>>(args[0])}) {
       return Expr<T>{std::visit(
           [](const auto &kx) {
             return Scalar<ResultType<decltype(kx)>>::RANGE;
           },
Expr<Type<TypeCategory::Integer, KIND>> FoldIntrinsicFunction(
   } else if (name == "storage_size") { // in bits
     if (auto info{
             characteristics::TypeAndShape::Characterize(args[0], context)}) {
       if (auto bytes{info->MeasureElementSizeInBytes(context, true)}) {
         return Expr<T>{
             Fold(context, Expr<T>{8} * ConvertToType<T>(std::move(*bytes)))};
       }
     }
+  } else if (name == "sum") {
+    return FoldSum<T>(context, std::move(funcRef));
   } else if (name == "ubound") {
     return UBOUND(context, std::move(funcRef));
   }
   // TODO:
-  // cshift, dot_product, eoshift,
-  // findloc, iall, iany, iparity, ibits, image_status, ishftc,
-  // matmul, maxloc, minloc, pack, product, reduce,
-  // sign, spread, sum, transfer, transpose, unpack
+  // cshift, dot_product, eoshift, findloc, ibits, image_status, ishftc,
+  // matmul, maxloc, minloc, not, pack, sign, spread, transfer, transpose,
+  // unpack
   return Expr<T>{std::move(funcRef)};
 }
 
 // Substitutes a bare type parameter reference with its value if it has one now
 // in an instantiation.  Bare LEN type parameters are substituted only when
 // the known value is constant.
 Expr<TypeParamInquiry::Result> FoldOperation(
     FoldingContext &context, TypeParamInquiry &&inquiry) {

flang/lib/Evaluate/fold-logical.cpp

 //===-- lib/Evaluate/fold-logical.cpp -------------------------------------===//
 //
 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
 // See https://llvm.org/LICENSE.txt for license information.
 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
 //
 //===----------------------------------------------------------------------===//
 
 #include "fold-implementation.h"
+#include "fold-reduction.h"
 #include "flang/Evaluate/check-expression.h"
 
 namespace Fortran::evaluate {
 
+// for ALL & ANY
+template <typename T>
+static Expr<T> FoldAllAny(FoldingContext &context, FunctionRef<T> &&ref,
+    Scalar<T> (Scalar<T>::*operation)(const Scalar<T> &) const,
+    Scalar<T> identity) {
+  static_assert(T::category == TypeCategory::Logical);
+  using Element = Scalar<T>;
+  std::optional<ConstantSubscript> dim;
+  if (std::optional<Constant<T>> array{
+          ProcessReductionArgs<T>(context, ref.arguments(), dim, identity,
+              /*ARRAY(MASK)=*/0, /*DIM=*/1)}) {
+    auto accumulator{[&](Element &element, const ConstantSubscripts &at) {
+      element = (element.*operation)(array->At(at));
+    }};
+    return Expr<T>{DoReduction(*array, dim, identity, accumulator)};
+  }
+  return Expr<T>{std::move(ref)};
+}
+
 template <int KIND>
 Expr<Type<TypeCategory::Logical, KIND>> FoldIntrinsicFunction(
     FoldingContext &context,
     FunctionRef<Type<TypeCategory::Logical, KIND>> &&funcRef) {
   using T = Type<TypeCategory::Logical, KIND>;
   ActualArguments &args{funcRef.arguments()};
   auto *intrinsic{std::get_if<SpecificIntrinsic>(&funcRef.proc().u)};
   CHECK(intrinsic);
   std::string name{intrinsic->name};
   if (name == "all") {
-    if (!args[1]) { // TODO: ALL(x,DIM=d)
-      if (const auto *constant{UnwrapConstantValue<T>(args[0])}) {
-        bool result{true};
-        for (const auto &element : constant->values()) {
-          if (!element.IsTrue()) {
-            result = false;
-            break;
-          }
-        }
-        return Expr<T>{result};
-      }
-    }
+    return FoldAllAny(
+        context, std::move(funcRef), &Scalar<T>::AND, Scalar<T>{true});
   } else if (name == "any") {
-    if (!args[1]) { // TODO: ANY(x,DIM=d)
-      if (const auto *constant{UnwrapConstantValue<T>(args[0])}) {
-        bool result{false};
-        for (const auto &element : constant->values()) {
-          if (element.IsTrue()) {
-            result = true;
-            break;
-          }
-        }
-        return Expr<T>{result};
-      }
-    }
+    return FoldAllAny(
+        context, std::move(funcRef), &Scalar<T>::OR, Scalar<T>{false});
   } else if (name == "associated") {
     bool gotConstant{true};
     const Expr<SomeType> *firstArgExpr{args[0]->UnwrapExpr()};
     if (!firstArgExpr || !IsNullPointer(*firstArgExpr)) {
       gotConstant = false;
     } else if (args[1]) { // There's a second argument
       const Expr<SomeType> *secondArgExpr{args[1]->UnwrapExpr()};
       if (!secondArgExpr || !IsNullPointer(*secondArgExpr)) {
return FoldMerge<T>(context, std::move(funcRef));
       name == "__builtin_ieee_support_sqrt" ||
       name == "__builtin_ieee_support_standard" ||
       name == "__builtin_ieee_support_subnormal" ||
       name == "__builtin_ieee_support_underflow_control") {
     return Expr<T>{true};
   }
   // TODO: btest, cshift, dot_product, eoshift, is_iostat_end,
   // is_iostat_eor, lge, lgt, lle, llt, logical, matmul, out_of_range,
-  // pack, parity, reduce, spread, transfer, transpose, unpack,
-  // extends_type_of, same_type_as
+  // pack, parity, spread, transfer, transpose, unpack, extends_type_of,
+  // same_type_as
   return Expr<T>{std::move(funcRef)};
 }
 
 template <typename T>
 Expr<LogicalResult> FoldOperation(
     FoldingContext &context, Relational<T> &&relation) {
   if (auto array{ApplyElementwise(context, relation,
           std::function<Expr<LogicalResult>(Expr<T> &&, Expr<T> &&)>{

flang/lib/Evaluate/fold-real.cpp

return FoldMaxvalMinval<T>(context, std::move(funcRef),
         RelationalOperator::GT, T::Scalar::HUGE().Negate());
   } else if (name == "merge") {
     return FoldMerge<T>(context, std::move(funcRef));
   } else if (name == "min") {
     return FoldMINorMAX(context, std::move(funcRef), Ordering::Less);
   } else if (name == "minval") {
     return FoldMaxvalMinval<T>(
         context, std::move(funcRef), RelationalOperator::LT, T::Scalar::HUGE());
+  } else if (name == "product") {
+    auto one{Scalar<T>::FromInteger(value::Integer<8>{1}).value};
+    return FoldProduct<T>(context, std::move(funcRef), one);
   } else if (name == "real") {
     if (auto *expr{args[0].value().UnwrapExpr()}) {
       return ToReal<KIND>(context, std::move(*expr));
     }
   } else if (name == "sign") {
     return FoldElementalIntrinsic<T, T, T>(
         context, std::move(funcRef), &Scalar<T>::SIGN);
+  } else if (name == "sum") {
+    return FoldSum<T>(context, std::move(funcRef));
   } else if (name == "tiny") {
     return Expr<T>{Scalar<T>::TINY()};
   }
   // TODO: cshift, dim, dot_product, eoshift, fraction, matmul,
-  // maxloc, minloc, modulo, nearest, norm2, pack, product,
-  // reduce, rrspacing, scale, set_exponent, spacing, spread,
-  // sum, transfer, transpose, unpack, bessel_jn (transformational) and
-  // bessel_yn (transformational)
+  // maxloc, minloc, modulo, nearest, norm2, pack, rrspacing, scale,
+  // set_exponent, spacing, spread, transfer, transpose, unpack,
+  // bessel_jn (transformational) and bessel_yn (transformational)
   return Expr<T>{std::move(funcRef)};
 }
 
 template <int KIND>
 Expr<Type<TypeCategory::Real, KIND>> FoldOperation(
     FoldingContext &context, ComplexComponent<KIND> &&x) {
   using Operand = Type<TypeCategory::Complex, KIND>;
   using Result = Type<TypeCategory::Real, KIND>;

flang/lib/Evaluate/fold-reduction.h

 
 #ifndef FORTRAN_EVALUATE_FOLD_REDUCTION_H_
 #define FORTRAN_EVALUATE_FOLD_REDUCTION_H_
 
 #include "fold-implementation.h"
 
 namespace Fortran::evaluate {
 
-// MAXVAL & MINVAL
+// Common preprocessing for reduction transformational intrinsic function
+// folding.  If the intrinsic can have DIM= &/or MASK= arguments, extract
+// and check them.  If a MASK= is present, apply it to the array data and
+// substitute identity values for elements corresponding to .FALSE. in
+// the mask.  If the result is present, the intrinsic call can be folded.
 template <typename T>
-Expr<T> FoldMaxvalMinval(FoldingContext &context, FunctionRef<T> &&ref,
-    RelationalOperator opr, Scalar<T> identity) {
-  static_assert(T::category == TypeCategory::Integer ||
-      T::category == TypeCategory::Real ||
-      T::category == TypeCategory::Character);
-  using Element = typename Constant<T>::Element;
-  auto &arg{ref.arguments()};
-  CHECK(arg.size() <= 3);
+static std::optional<Constant<T>> ProcessReductionArgs(FoldingContext &context,
+    ActualArguments &arg, std::optional<ConstantSubscript> &dim,
+    const Scalar<T> &identity, int arrayIndex,
+    std::optional<std::size_t> dimIndex = std::nullopt,
+    std::optional<std::size_t> maskIndex = std::nullopt) {
   if (arg.empty()) {
-    return Expr<T>{std::move(ref)};
-  }
-  Constant<T> *array{Folder<T>{context}.Folding(arg[0])};
-  if (!array || array->Rank() < 1) {
-    return Expr<T>{std::move(ref)};
+    return std::nullopt;
   }
-  std::optional<ConstantSubscript> dim;
-  if (arg.size() >= 2 && arg[1]) {
-    if (auto *dimConst{Folder<SubscriptInteger>{context}.Folding(arg[1])}) {
+  Constant<T> *folded{Folder<T>{context}.Folding(arg[arrayIndex])};
+  if (!folded || folded->Rank() < 1) {
+    return std::nullopt;
+  }
+  if (dimIndex && arg.size() >= *dimIndex + 1 && arg[*dimIndex]) {
+    if (auto *dimConst{
+            Folder<SubscriptInteger>{context}.Folding(arg[*dimIndex])}) {
       if (auto dimScalar{dimConst->GetScalarValue()}) {
         dim.emplace(dimScalar->ToInt64());
-        if (*dim < 1 || *dim > array->Rank()) {
+        if (*dim < 1 || *dim > folded->Rank()) {
           context.messages().Say(
               "DIM=%jd is not valid for an array of rank %d"_err_en_US,
-              static_cast<std::intmax_t>(*dim), array->Rank());
+              static_cast<std::intmax_t>(*dim), folded->Rank());
           dim.reset();
         }
       }
     }
     if (!dim) {
-      return Expr<T>{std::move(ref)};
-    }
+      return std::nullopt;
     }
-  Constant<LogicalResult> *mask{};
-  if (arg.size() >= 3 && arg[2]) {
-    mask = Folder<LogicalResult>{context}.Folding(arg[2]);
-    if (!mask) {
-      return Expr<T>{std::move(ref)};
-    }
-    if (!CheckConformance(context.messages(), AsShape(array->shape()),
+  if (maskIndex && arg.size() >= *maskIndex + 1 && arg[*maskIndex]) {
+    if (Constant<LogicalResult> *
+        mask{Folder<LogicalResult>{context}.Folding(arg[*maskIndex])}) {
+      if (CheckConformance(context.messages(), AsShape(folded->shape()),
               AsShape(mask->shape()),
               CheckConformanceFlags::RightScalarExpandable, "ARRAY=", "MASK=")
               .value_or(false)) {
-      return Expr<T>{std::move(ref)};
+        // Apply the mask in place to the array
+        std::size_t n{folded->size()};
+        std::vector<typename Constant<T>::Element> elements;
+        if (auto scalarMask{mask->GetScalarValue()}) {
+          if (scalarMask->IsTrue()) {
+            return Constant<T>{*folded};
+          } else { // MASK=.FALSE.
+            elements = std::vector<typename Constant<T>::Element>(n, identity);
+          }
+        } else { // mask is an array; test its elements
+          elements = std::vector<typename Constant<T>::Element>(n, identity);
+          ConstantSubscripts at{folded->lbounds()};
+          for (std::size_t j{0}; j < n; ++j, folded->IncrementSubscripts(at)) {
+            if (mask->values()[j].IsTrue()) {
+              elements[j] = folded->At(at);
+            }
           }
         }
-  // Do it
-  ConstantSubscripts at{array->lbounds()}, maskAt;
-  bool maskAllFalse{false};
-  if (mask) {
-    if (auto scalar{mask->GetScalarValue()}) {
-      if (scalar->IsTrue()) {
-        mask = nullptr; // all .TRUE.
+        if constexpr (T::category == TypeCategory::Character) {
+          return Constant<T>{static_cast<ConstantSubscript>(identity.size()),
+              std::move(elements), ConstantSubscripts{folded->shape()}};
         } else {
-        maskAllFalse = true;
+          return Constant<T>{
+              std::move(elements), ConstantSubscripts{folded->shape()}};
         }
       } else {
-      maskAt = mask->lbounds();
+        return std::nullopt;
       }
+    } else {
+      return std::nullopt;
     }
-  std::vector<Element> result;
-  ConstantSubscripts resultShape; // empty -> scalar
-  // Internal function to accumulate into result.back().
-  auto Accumulate{[&]() {
-    if (!maskAllFalse && (maskAt.empty() || mask->At(maskAt).IsTrue())) {
-      Expr<LogicalResult> test{
-          PackageRelation(opr, Expr<T>{Constant<T>{array->At(at)}},
-              Expr<T>{Constant<T>{result.back()}})};
-      auto folded{GetScalarConstantValue<LogicalResult>(
-          test.Rewrite(context, std::move(test)))};
-      CHECK(folded.has_value());
-      if (folded->IsTrue()) {
-        result.back() = array->At(at);
-      }
-    }
-  }};
+  } else {
+    return Constant<T>{*folded};
+
+// Generalized reduction to an array of one dimension fewer (w/ DIM=)
+// or to a scalar (w/o DIM=).
+template <typename T, typename ACCUMULATOR>
+static Constant<T> DoReduction(const Constant<T> &array,
+    std::optional<ConstantSubscript> &dim, const Scalar<T> &identity,
+    ACCUMULATOR &accumulator) {
+  ConstantSubscripts at{array.lbounds()};
+  std::vector<typename Constant<T>::Element> elements;
+  ConstantSubscripts resultShape; // empty -> scalar
   if (dim) { // DIM= is present, so result is an array
-    resultShape = array->shape();
+    resultShape = array.shape();
     resultShape.erase(resultShape.begin() + (*dim - 1));
-    ConstantSubscript dimExtent{array->shape().at(*dim - 1)};
+    ConstantSubscript dimExtent{array.shape().at(*dim - 1)};
     ConstantSubscript &dimAt{at[*dim - 1]};
     ConstantSubscript dimLbound{dimAt};
-    ConstantSubscript *maskDimAt{maskAt.empty() ? nullptr : &maskAt[*dim - 1]};
-    ConstantSubscript maskLbound{maskDimAt ? *maskDimAt : 0};
     for (auto n{GetSize(resultShape)}; n-- > 0;
-         IncrementSubscripts(at, array->shape())) {
+         IncrementSubscripts(at, array.shape())) {
       dimAt = dimLbound;
-      if (maskDimAt) {
-        *maskDimAt = maskLbound;
+      elements.push_back(identity);
+      for (ConstantSubscript j{0}; j < dimExtent; ++j, ++dimAt) {
+        accumulator(elements.back(), at);
       }
-      result.push_back(identity);
-      for (ConstantSubscript j{0}; j < dimExtent;
-           ++j, ++dimAt, maskDimAt && ++*maskDimAt) {
-        Accumulate();
-      }
-      if (maskDimAt) {
-        IncrementSubscripts(maskAt, mask->shape());
-      }
-    }
-  } else { // no DIM=, result is scalar
-    result.push_back(identity);
-    for (auto n{array->size()}; n-- > 0;
-         IncrementSubscripts(at, array->shape())) {
-      Accumulate();
-      if (!maskAt.empty()) {
-        IncrementSubscripts(maskAt, mask->shape());
-      }
-    }
+  } else { // no DIM=, result is scalar
+    elements.push_back(identity);
+    for (auto n{array.size()}; n-- > 0;
+         IncrementSubscripts(at, array.shape())) {
+      accumulator(elements.back(), at);
+  if constexpr (T::category == TypeCategory::Character) {
+    return {static_cast<ConstantSubscript>(identity.size()),
+        std::move(elements), std::move(resultShape)};
+  } else {
+    return {std::move(elements), std::move(resultShape)};
+  }
+}
+
+// MAXVAL & MINVAL
+template <typename T>
+static Expr<T> FoldMaxvalMinval(FoldingContext &context, FunctionRef<T> &&ref,
+    RelationalOperator opr, const Scalar<T> &identity) {
+  static_assert(T::category == TypeCategory::Integer ||
+      T::category == TypeCategory::Real ||
+      T::category == TypeCategory::Character);
+  using Element = Scalar<T>; // pmk: was typename Constant<T>::Element;
+  std::optional<ConstantSubscript> dim;
+  if (std::optional<Constant<T>> array{
+          ProcessReductionArgs<T>(context, ref.arguments(), dim, identity,
+              /*ARRAY=*/0, /*DIM=*/1, /*MASK=*/2)}) {
+    auto accumulator{[&](Element &element, const ConstantSubscripts &at) {
+      Expr<LogicalResult> test{PackageRelation(opr,
+          Expr<T>{Constant<T>{array->At(at)}}, Expr<T>{Constant<T>{element}})};
+      auto folded{GetScalarConstantValue<LogicalResult>(
+          test.Rewrite(context, std::move(test)))};
+      CHECK(folded.has_value());
+      if (folded->IsTrue()) {
+        element = array->At(at);
+      }
+    }};
+    return Expr<T>{DoReduction(*array, dim, identity, accumulator)};
+  return Expr<T>{std::move(ref)};
+
+// PRODUCT
+template <typename T>
+static Expr<T> FoldProduct(
+    FoldingContext &context, FunctionRef<T> &&ref, Scalar<T> identity) {
+  static_assert(T::category == TypeCategory::Integer ||
+      T::category == TypeCategory::Real ||
+      T::category == TypeCategory::Complex);
+  using Element = typename Constant<T>::Element;
+  std::optional<ConstantSubscript> dim;
+  if (std::optional<Constant<T>> array{
+          ProcessReductionArgs<T>(context, ref.arguments(), dim, identity,
+              /*ARRAY=*/0, /*DIM=*/1, /*MASK=*/2)}) {
+    bool overflow{false};
+    auto accumulator{[&](Element &element, const ConstantSubscripts &at) {
+      if constexpr (T::category == TypeCategory::Integer) {
+        auto prod{element.MultiplySigned(array->At(at))};
+        overflow |= prod.SignedMultiplicationOverflowed();
+        element = prod.lower;
+      } else { // Real & Complex
+        auto prod{element.Multiply(array->At(at))};
+        overflow |= prod.flags.test(RealFlag::Overflow);
+        element = prod.value;
       }
-  if constexpr (T::category == TypeCategory::Character) {
-    return Expr<T>{Constant<T>{static_cast<ConstantSubscript>(identity.size()),
-        std::move(result), std::move(resultShape)}};
+    }};
+    if (overflow) {
+      context.messages().Say(
+          "PRODUCT() of %s data overflowed"_en_US, T::AsFortran());
     } else {
-    return Expr<T>{Constant<T>{std::move(result), std::move(resultShape)}};
+      return Expr<T>{DoReduction(*array, dim, identity, accumulator)};
+    }
   }
+  return Expr<T>{std::move(ref)};
+}
+
+// SUM
+template <typename T>
+static Expr<T> FoldSum(FoldingContext &context, FunctionRef<T> &&ref) {
+  static_assert(T::category == TypeCategory::Integer ||
+      T::category == TypeCategory::Real ||
+      T::category == TypeCategory::Complex);
+  using Element = typename Constant<T>::Element;
+  std::optional<ConstantSubscript> dim;
+  Element identity{}, correction{};
+  if (std::optional<Constant<T>> array{
+          ProcessReductionArgs<T>(context, ref.arguments(), dim, identity,
+              /*ARRAY=*/0, /*DIM=*/1, /*MASK=*/2)}) {
+    bool overflow{false};
+    auto accumulator{[&](Element &element, const ConstantSubscripts &at) {
+      if constexpr (T::category == TypeCategory::Integer) {
+        auto sum{element.AddSigned(array->At(at))};
+        overflow |= sum.overflow;
+        element = sum.value;
+      } else { // Real & Complex: use Kahan summation
+        auto next{array->At(at).Add(correction)};
+        overflow |= next.flags.test(RealFlag::Overflow);
+        auto sum{element.Add(next.value)};
+        overflow |= sum.flags.test(RealFlag::Overflow);
+        // correction = (sum - element) - next; algebraically zero
+        correction =
+            sum.value.Subtract(element).value.Subtract(next.value).value;
+        element = sum.value;
+      }
+    }};
+    if (overflow) {
+      context.messages().Say(
+          "SUM() of %s data overflowed"_en_US, T::AsFortran());
+    } else {
+      return Expr<T>{DoReduction(*array, dim, identity, accumulator)};
+    }
+  }
+  return Expr<T>{std::move(ref)};
 }
 
 } // namespace Fortran::evaluate
 #endif // FORTRAN_EVALUATE_FOLD_REDUCTION_H_

flang/test/Evaluate/folding20.f90

 ! RUN: %S/test_folding.sh %s %t %flang_fc1
 ! REQUIRES: shell
-! Tests intrinsic MAXVAL/MINVAL function folding
+! Tests reduction intrinsic function folding
 module m
+  implicit none
+  integer, parameter :: intmatrix(*,*) = reshape([1, 2, 3, 4, 5, 6], [2, 3])
+  logical, parameter :: odds(2,3) = mod(intmatrix, 2) == 1
+  character(*), parameter :: chmatrix(*,*) = reshape(['abc', 'def', 'ghi', 'jkl', 'mno', 'pqr'], [2, 3])
+
+  logical, parameter :: test_allidentity = all([Logical::])
+  logical, parameter :: test_all = .not. all(odds)
+  logical, parameter :: test_alldim1 = all(.not. all(odds,1))
+  logical, parameter :: test_alldim2 = all(all(odds,2) .eqv. [.true., .false.])
+  logical, parameter :: test_anyidentity = .not. any([Logical::])
+  logical, parameter :: test_any = any(odds)
+  logical, parameter :: test_anydim1 = all(any(odds,1))
+  logical, parameter :: test_anydim2 = all(any(odds,2) .eqv. [.true., .false.])
+
+  logical, parameter :: test_iallidentity = iall([integer::]) == -1
+  logical, parameter :: test_iall = iall(intmatrix) == 0
+  logical, parameter :: test_iall_masked = iall(intmatrix,odds) == 1
+  logical, parameter :: test_ialldim1 = all(iall(intmatrix,dim=1) == [0, 0, 4])
+  logical, parameter :: test_ialldim2 = all(iall(intmatrix,dim=2) == [1, 0])
+  logical, parameter :: test_ianyidentity = iany([integer::]) == 0
+  logical, parameter :: test_iany = iany(intmatrix) == 7
+  logical, parameter :: test_iany_masked = iany(intmatrix,odds) == 7
+  logical, parameter :: test_ianydim1 = all(iany(intmatrix,dim=1) == [3, 7, 7])
+  logical, parameter :: test_ianydim2 = all(iany(intmatrix,dim=2) == [7, 6])
+  logical, parameter :: test_iparityidentity = iparity([integer::]) == 0
+  logical, parameter :: test_iparity = iparity(intmatrix) == 7
+  logical, parameter :: test_iparity_masked = iparity(intmatrix,odds) == 7
+  logical, parameter :: test_iparitydim1 = all(iparity(intmatrix,dim=1) == [3, 7, 3])
+  logical, parameter :: test_iparitydim2 = all(iparity(intmatrix,dim=2) == [7, 0])
+
   logical, parameter :: test_imaxidentity = maxval([integer::]) == -huge(0) - 1
   logical, parameter :: test_iminidentity = minval([integer::]) == huge(0)
-  integer, parameter :: intmatrix(*,*) = reshape([1, 2, 3, 4, 5, 6], [2, 3])
   logical, parameter :: test_imaxval = maxval(intmatrix) == 6
   logical, parameter :: test_iminval = minval(intmatrix) == 1
-  logical, parameter :: odds(2,3) = mod(intmatrix, 2) == 1
   logical, parameter :: test_imaxval_masked = maxval(intmatrix,odds) == 5
   logical, parameter :: test_iminval_masked = minval(intmatrix,.not.odds) == 2
   logical, parameter :: test_rmaxidentity = maxval([real::]) == -huge(0.0)
   logical, parameter :: test_rminidentity = minval([real::]) == huge(0.0)
   logical, parameter :: test_rmaxval = maxval(real(intmatrix)) == 6.0
   logical, parameter :: test_rminval = minval(real(intmatrix)) == 1.0
   logical, parameter :: test_rmaxval_scalar_mask = maxval(real(intmatrix), .true.) == 6.0
   logical, parameter :: test_rminval_scalar_mask = minval(real(intmatrix), .false.) == huge(0.0)
-  character(*), parameter :: chmatrix(*,*) = reshape(['abc', 'def', 'ghi', 'jkl', 'mno', 'pqr'], [2, 3])
   logical, parameter :: test_cmaxlen = len(maxval([character*4::])) == 4
   logical, parameter :: test_cmaxidentity = maxval([character*4::]) == repeat(char(0), 4)
   logical, parameter :: test_cminidentity = minval([character*4::]) == repeat(char(127), 4)
   logical, parameter :: test_cmaxval = maxval(chmatrix) == 'pqr'
   logical, parameter :: test_cminval = minval(chmatrix) == 'abc'
-  logical, parameter :: test_dim1 = all(maxval(intmatrix,dim=1) == [2, 4, 6])
-  logical, parameter :: test_dim2 = all(minval(intmatrix,dim=2,mask=odds) == [1, huge(0)])
+  logical, parameter :: test_maxvaldim1 = all(maxval(intmatrix,dim=1) == [2, 4, 6])
+  logical, parameter :: test_minvaldim2 = all(minval(intmatrix,dim=2,mask=odds) == [1, huge(0)])
+
+  logical, parameter :: test_iproductidentity = product([integer::]) == 1
+  logical, parameter :: test_iproduct = product(intmatrix) == 720
+  logical, parameter :: test_iproduct_masked = product(intmatrix,odds) == 15
+  logical, parameter :: test_productdim1 = all(product(intmatrix,dim=1) == [2, 12, 30])
+  logical, parameter :: test_productdim2 = all(product(intmatrix,dim=2) == [15, 48])
+  logical, parameter :: test_rproductidentity = product([real::]) == 1.
+  logical, parameter :: test_rproduct = product(real(intmatrix)) == 720.
+  logical, parameter :: test_cproductidentity = product([complex::]) == (1.,0.)
+  logical, parameter :: test_cproduct = product(cmplx(intmatrix,-intmatrix)) == (0.,5760.)
+
+  logical, parameter :: test_isumidentity = sum([integer::]) == 0
+  logical, parameter :: test_isum = sum(intmatrix) == 21
+  logical, parameter :: test_isum_masked = sum(intmatrix,odds) == 9
+  logical, parameter :: test_sumdim1 = all(sum(intmatrix,dim=1) == [3, 7, 11])
+  logical, parameter :: test_sumdim2 = all(sum(intmatrix,dim=2) == [9, 12])
+  logical, parameter :: test_rsumidentity = sum([real::]) == 0.
+  logical, parameter :: test_rsum = sum(real(intmatrix)) == 21.
+  logical, parameter :: test_csumidentity = sum([complex::]) == (0.,0.)
+  logical, parameter :: test_csum = sum(cmplx(intmatrix,-intmatrix)) == (21.,-21.)
 end