[flang] Implement NORM2 in the runtime

[klausler]

Summary

Implement the reduction transformational intrinsic function NORM2 in
the runtime, using infrastructure already in place for MAXVAL & al.

[PeteSteinfeld]

All builds, tests, and looks good.

flang/runtime/extrema.cpp

 //===-- runtime/extrema.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
 //
 //===----------------------------------------------------------------------===//
 
 // Implements MAXLOC, MINLOC, MAXVAL, & MINVAL for all required operand types
-// and shapes and (for MAXLOC & MINLOC) result integer kinds.
+// and shapes and (for MAXLOC & MINLOC) result integer kinds.  Also implements
+// NORM2 using common infrastructure.
 
 #include "character.h"
 #include "reduction-templates.h"
 #include "reduction.h"
 #include "flang/Common/long-double.h"
+#include <algorithm>
 #include <cinttypes>
+#include <cmath>
+#include <optional>
 
 namespace Fortran::runtime {
+
 // MAXLOC & MINLOC
 
 template <typename T, bool IS_MAX, bool BACK> struct NumericCompare {
   using Type = T;
   explicit NumericCompare(std::size_t /*elemLen; ignored*/) {}
   bool operator()(const T &value, const T &previous) const {
     if (value == previous) {
       return BACK;

 template <TypeCategory CAT, int KIND, bool IS_MAXVAL>
 inline CppTypeFor<CAT, KIND> TotalNumericMaxOrMin(const Descriptor &x,
     const char *source, int line, int dim, const Descriptor *mask,
     const char *intrinsic) {
   return GetTotalReduction<CAT, KIND>(x, source, line, dim, mask,
       NumericExtremumAccumulator<CAT, KIND, IS_MAXVAL>{x}, intrinsic);
 }
 
-template <TypeCategory CAT, int KIND, bool IS_MAXVAL,
-    template <TypeCategory, int, bool> class ACCUMULATOR>
-static void DoMaxOrMin(Descriptor &result, const Descriptor &x, int dim,
+template <TypeCategory CAT, int KIND, typename ACCUMULATOR>
+static void DoMaxMinNorm2(Descriptor &result, const Descriptor &x, int dim,
     const Descriptor *mask, const char *intrinsic, Terminator &terminator) {
   using Type = CppTypeFor<CAT, KIND>;
+  ACCUMULATOR accumulator{x};
   if (dim == 0 || x.rank() == 1) {
     // Total reduction
     result.Establish(x.type(), x.ElementBytes(), nullptr, 0, nullptr,
         CFI_attribute_allocatable);
     if (int stat{result.Allocate()}) {
       terminator.Crash(
           "%s: could not allocate memory for result; STAT=%d", intrinsic, stat);
     }
-    ACCUMULATOR<CAT, KIND, IS_MAXVAL> accumulator{x};
     DoTotalReduction<Type>(x, dim, mask, accumulator, intrinsic, terminator);
     accumulator.GetResult(result.OffsetElement<Type>());
   } else {
     // Partial reduction
-    using Accumulator = ACCUMULATOR<CAT, KIND, IS_MAXVAL>;
-    Accumulator accumulator{x};
-    PartialReduction<Accumulator, CAT, KIND>(
+    PartialReduction<ACCUMULATOR, CAT, KIND>(
         result, x, dim, mask, terminator, intrinsic, accumulator);
   }
 }
 
 template <TypeCategory CAT, bool IS_MAXVAL> struct MaxOrMinHelper {
   template <int KIND> struct Functor {
     void operator()(Descriptor &result, const Descriptor &x, int dim,
         const Descriptor *mask, const char *intrinsic,
         Terminator &terminator) const {
-      DoMaxOrMin<CAT, KIND, IS_MAXVAL, NumericExtremumAccumulator>(
+      DoMaxMinNorm2<CAT, KIND,
+          NumericExtremumAccumulator<CAT, KIND, IS_MAXVAL>>(
           result, x, dim, mask, intrinsic, terminator);
     }
   };
 };
 
 template <bool IS_MAXVAL>
 inline void NumericMaxOrMin(Descriptor &result, const Descriptor &x, int dim,
     const char *source, int line, const Descriptor *mask,
ApplyFloatingPointKind<
         MaxOrMinHelper<TypeCategory::Real, IS_MAXVAL>::template Functor, void>(
         type->second, terminator, result, x, dim, mask, intrinsic, terminator);
     break;
   default:
     terminator.Crash("%s: bad type code %d", intrinsic, x.type().raw());
   }
 }
 
-template <TypeCategory, int KIND, bool IS_MAXVAL>
-class CharacterExtremumAccumulator {
+template <int KIND, bool IS_MAXVAL> class CharacterExtremumAccumulator {
 public:
   using Type = CppTypeFor<TypeCategory::Character, KIND>;
   explicit CharacterExtremumAccumulator(const Descriptor &array)
       : array_{array}, charLen_{array_.ElementBytes() / KIND} {}
   void Reinitialize() { extremum_ = nullptr; }
   template <typename A> void GetResult(A *p, int /*zeroBasedDim*/ = -1) const {
     static_assert(std::is_same_v<A, Type>);
     if (extremum_) {
private:
   const Type *extremum_{nullptr};
 };
 
 template <bool IS_MAXVAL> struct CharacterMaxOrMinHelper {
   template <int KIND> struct Functor {
     void operator()(Descriptor &result, const Descriptor &x, int dim,
         const Descriptor *mask, const char *intrinsic,
         Terminator &terminator) const {
-      DoMaxOrMin<TypeCategory::Character, KIND, IS_MAXVAL,
-          CharacterExtremumAccumulator>(
+      DoMaxMinNorm2<TypeCategory::Character, KIND,
+          CharacterExtremumAccumulator<KIND, IS_MAXVAL>>(
           result, x, dim, mask, intrinsic, terminator);
     }
   };
 };
 
 template <bool IS_MAXVAL>
 inline void CharacterMaxOrMin(Descriptor &result, const Descriptor &x, int dim,
     const char *source, int line, const Descriptor *mask,
void RTNAME(MinvalDim)(Descriptor &result, const Descriptor &x, int dim,
     const char *source, int line, const Descriptor *mask) {
   if (x.type().IsCharacter()) {
     CharacterMaxOrMin<false>(result, x, dim, source, line, mask, "MINVAL");
   } else {
     NumericMaxOrMin<false>(result, x, dim, source, line, mask, "MINVAL");
   }
 }
 } // extern "C"
+
+// NORM2
+
+template <int KIND> class Norm2Accumulator {
+public:
+  using Type = CppTypeFor<TypeCategory::Real, KIND>;
+  // Use at least double precision for accumulators
+  using AccumType = CppTypeFor<TypeCategory::Real, std::max(KIND, 8)>;
+  explicit Norm2Accumulator(const Descriptor &array) : array_{array} {}
+  void Reinitialize() { max_ = sum_ = 0; }
+  template <typename A> void GetResult(A *p, int /*zeroBasedDim*/ = -1) const {
+    // m * sqrt(1 + sum((others(:)/m)**2))
+    *p = static_cast<Type>(max_ * std::sqrt(1 + sum_));
+  }
+  bool Accumulate(Type x) {
+    auto absX{AccumType{std::abs(x)}};
+    if (!max_) {
+      max_ = x;
+    } else if (absX > max_) {
+      auto t{max_ / absX}; // < 1.0
+      auto tsq{t * t};
+      sum_ *= tsq; // scale sum to reflect change to the max
+      sum_ += tsq; // include a term for the previous max
+      max_ = absX;
+    } else { // absX <= max_
+      auto t{absX / max_};
+      sum_ += t * t;
+    }
+    return true;
+  }
+  template <typename A> bool AccumulateAt(const SubscriptValue at[]) {
+    return Accumulate(*array_.Element<A>(at));
+  }
+
+private:
+  const Descriptor &array_;
+  AccumType max_{0}; // value (m) with largest magnitude
+  AccumType sum_{0}; // sum((others(:)/m)**2)
+};
+
+template <int KIND> struct Norm2Helper {
+  void operator()(Descriptor &result, const Descriptor &x, int dim,
+      const Descriptor *mask, Terminator &terminator) const {
+    DoMaxMinNorm2<TypeCategory::Real, KIND, Norm2Accumulator<KIND>>(
+        result, x, dim, mask, "NORM2", terminator);
+  }
+};
+
+extern "C" {
+// TODO: REAL(2 & 3)
+CppTypeFor<TypeCategory::Real, 4> RTNAME(Norm2_4)(const Descriptor &x,
+    const char *source, int line, int dim, const Descriptor *mask) {
+  return GetTotalReduction<TypeCategory::Real, 4>(
+      x, source, line, dim, mask, Norm2Accumulator<4>{x}, "NORM2");
+}
+CppTypeFor<TypeCategory::Real, 8> RTNAME(Norm2_8)(const Descriptor &x,
+    const char *source, int line, int dim, const Descriptor *mask) {
+  return GetTotalReduction<TypeCategory::Real, 8>(
+      x, source, line, dim, mask, Norm2Accumulator<8>{x}, "NORM2");
+}
+#if LONG_DOUBLE == 80
+CppTypeFor<TypeCategory::Real, 10> RTNAME(Norm2_10)(const Descriptor &x,
+    const char *source, int line, int dim, const Descriptor *mask) {
+  return GetTotalReduction<TypeCategory::Real, 10>(
+      x, source, line, dim, mask, Norm2Accumulator<10>{x}, "NORM2");
+}
+#elif LONG_DOUBLE == 128
+CppTypeFor<TypeCategory::Real, 16> RTNAME(Norm2_16)(const Descriptor &x,
+    const char *source, int line, int dim, const Descriptor *mask) {
+  return GetTotalReduction<TypeCategory::Real, 16>(
+      x, source, line, dim, mask, Norm2Accumulator<16>{x}, "NORM2");
+}
+#endif
+
+void RTNAME(Norm2Dim)(Descriptor &result, const Descriptor &x, int dim,
+    const char *source, int line, const Descriptor *mask) {
+  Terminator terminator{source, line};
+  auto type{x.type().GetCategoryAndKind()};
+  RUNTIME_CHECK(terminator, type);
+  if (type->first == TypeCategory::Real) {
+    ApplyFloatingPointKind<Norm2Helper, void>(
+        type->second, terminator, result, x, dim, mask, terminator);
+  } else {
+    terminator.Crash("NORM2: bad type code %d", x.type().raw());
+  }
+}
+} // extern "C"
 } // namespace Fortran::runtime

flang/runtime/reduction.h

 void RTNAME(MinvalCharacter)(Descriptor &, const Descriptor &,
     const char *source, int line, const Descriptor *mask = nullptr);
 
 void RTNAME(MaxvalDim)(Descriptor &, const Descriptor &, int dim,
     const char *source, int line, const Descriptor *mask = nullptr);
 void RTNAME(MinvalDim)(Descriptor &, const Descriptor &, int dim,
     const char *source, int line, const Descriptor *mask = nullptr);
 
+// NORM2
+float RTNAME(Norm2_2)(const Descriptor &, const char *source, int line,
+    int dim = 0, const Descriptor *mask = nullptr);
+float RTNAME(Norm2_3)(const Descriptor &, const char *source, int line,
+    int dim = 0, const Descriptor *mask = nullptr);
+float RTNAME(Norm2_4)(const Descriptor &, const char *source, int line,
+    int dim = 0, const Descriptor *mask = nullptr);
+double RTNAME(Norm2_8)(const Descriptor &, const char *source, int line,
+    int dim = 0, const Descriptor *mask = nullptr);
+long double RTNAME(Norm2_10)(const Descriptor &, const char *source, int line,
+    int dim = 0, const Descriptor *mask = nullptr);
+long double RTNAME(Norm2_16)(const Descriptor &, const char *source, int line,
+    int dim = 0, const Descriptor *mask = nullptr);
+void RTNAME(Norm2Dim)(Descriptor &, const Descriptor &, int dim,
+    const char *source, int line, const Descriptor *mask = nullptr);
+
 // ALL, ANY, COUNT, & PARITY logical reductions
 bool RTNAME(All)(const Descriptor &, const char *source, int line, int dim = 0);
 void RTNAME(AllDim)(Descriptor &result, const Descriptor &, int dim,
     const char *source, int line);
 bool RTNAME(Any)(const Descriptor &, const char *source, int line, int dim = 0);
 void RTNAME(AnyDim)(Descriptor &result, const Descriptor &, int dim,
     const char *source, int line);
 std::int64_t RTNAME(Count)(

flang/runtime/reduction.cpp

 //===-- runtime/reduction.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
 //
 //===----------------------------------------------------------------------===//
 
 // Implements ALL, ANY, COUNT, IPARITY, & PARITY for all required operand
 // types and shapes.
 //
 // FINDLOC, SUM, and PRODUCT are in their own eponymous source files;
-// MAXLOC, MINLOC, MAXVAL, and MINVAL are in extrema.cpp.
+// NORM2, MAXLOC, MINLOC, MAXVAL, and MINVAL are in extrema.cpp.
 
 #include "reduction.h"
 #include "reduction-templates.h"
 #include <cinttypes>
 
 namespace Fortran::runtime {
 
 // IPARITY()

flang/unittests/RuntimeGTest/Reduction.cpp

TEST(Reductions, DimMaskProductInt4) {
   EXPECT_EQ(prod.GetDimension(0).Extent(), 3);
   EXPECT_EQ(*prod.ZeroBasedIndexedElement<std::int32_t>(0), 1);
   EXPECT_EQ(*prod.ZeroBasedIndexedElement<std::int32_t>(1), 4);
   EXPECT_EQ(*prod.ZeroBasedIndexedElement<std::int32_t>(2), 30);
   EXPECT_EQ(RTNAME(SumInteger4)(prod, __FILE__, __LINE__), 35);
   prod.Destroy();
 }
 
-TEST(Reductions, DoubleMaxMin) {
+TEST(Reductions, DoubleMaxMinNorm2) {
   std::vector<int> shape{3, 4, 2}; // rows, columns, planes
   //   0  -3   6  -9     12 -15  18 -21
   //  -1   4  -7  10    -13  16 -19  22
   //   2  -5   8 -11     14 -17  20  22   <- note last two are equal to test
   //   BACK=
-  auto array{MakeArray<TypeCategory::Real, 8>(shape,
-      std::vector<double>{0, -1, 2, -3, 4, -5, 6, -7, 8, -9, 10, -11, 12, -13,
-          14, -15, 16, -17, 18, -19, 20, -21, 22, 22})};
+  std::vector<double> rawData{0, -1, 2, -3, 4, -5, 6, -7, 8, -9, 10, -11, 12,
+      -13, 14, -15, 16, -17, 18, -19, 20, -21, 22, 22};
+  auto array{MakeArray<TypeCategory::Real, 8>(shape, rawData)};
   EXPECT_EQ(RTNAME(MaxvalReal8)(*array, __FILE__, __LINE__), 22.0);
   EXPECT_EQ(RTNAME(MinvalReal8)(*array, __FILE__, __LINE__), -21.0);
+  double naiveNorm2{0};
+  for (auto x : rawData) {
+    naiveNorm2 += x * x;
+  }
+  naiveNorm2 = std::sqrt(naiveNorm2);
+  double norm2Error{
+      std::abs(naiveNorm2 - RTNAME(Norm2_8)(*array, __FILE__, __LINE__))};
+  EXPECT_LE(norm2Error, 0.000001 * naiveNorm2);
   StaticDescriptor<2> statDesc;
   Descriptor &loc{statDesc.descriptor()};
   RTNAME(Maxloc)
   (loc, *array, /*KIND=*/8, __FILE__, __LINE__, /*MASK=*/nullptr,
       /*BACK=*/false);
   EXPECT_EQ(loc.rank(), 1);
   EXPECT_EQ(loc.type().raw(), (TypeCode{TypeCategory::Integer, 8}.raw()));
   EXPECT_EQ(loc.GetDimension(0).LowerBound(), 1);