flang/runtime/extrema.cpp

   1 //===-- runtime/extrema.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 MAXLOC, MINLOC, MAXVAL, & MINVAL for all required operand types
  10 // and shapes and (for MAXLOC & MINLOC) result integer kinds.  Also implements
  11 // NORM2 using common infrastructure.
  12
  13 #include "reduction-templates.h"
  14 #include "flang/Common/float128.h"
  15 #include "flang/Runtime/character.h"
  16 #include "flang/Runtime/reduction.h"
  17 #include <algorithm>
  18 #include <cfloat>
  19 #include <cinttypes>
  20 #include <cmath>
  21 #include <optional>
  22
  23 namespace Fortran::runtime {
  24
  25 // MAXLOC & MINLOC
  26
  27 template <typename T, bool IS_MAX, bool BACK> struct NumericCompare {
  28   using Type = T;
  29   explicit NumericCompare(std::size_t /*elemLen; ignored*/) {}
  30   bool operator()(const T &value, const T &previous) const {
  31     if (value == previous) {
  32       return BACK;
  33     } else if constexpr (IS_MAX) {
  34       return value > previous;
  35     } else {
  36       return value < previous;
  37     }
  38   }
  39 };
  40
  41 template <typename T, bool IS_MAX, bool BACK> class CharacterCompare {
  42 public:
  43   using Type = T;
  44   explicit CharacterCompare(std::size_t elemLen)
  45       : chars_{elemLen / sizeof(T)} {}
  46   bool operator()(const T &value, const T &previous) const {
  47     int cmp{CharacterScalarCompare<T>(&value, &previous, chars_, chars_)};
  48     if (cmp == 0) {
  49       return BACK;
  50     } else if constexpr (IS_MAX) {
  51       return cmp > 0;
  52     } else {
  53       return cmp < 0;
  54     }
  55   }
  56
  57 private:
  58   std::size_t chars_;
  59 };
  60
  61 template <typename COMPARE> class ExtremumLocAccumulator {
  62 public:
  63   using Type = typename COMPARE::Type;
  64   ExtremumLocAccumulator(const Descriptor &array)
  65       : array_{array}, argRank_{array.rank()}, compare_{array.ElementBytes()} {
  66     Reinitialize();
  67   }
  68   void Reinitialize() {
  69     // per standard: result indices are all zero if no data
  70     for (int j{0}; j < argRank_; ++j) {
  71       extremumLoc_[j] = 0;
  72     }
  73     previous_ = nullptr;
  74   }
  75   int argRank() const { return argRank_; }
  76   template <typename A> void GetResult(A *p, int zeroBasedDim = -1) {
  77     if (zeroBasedDim >= 0) {
  78       *p = extremumLoc_[zeroBasedDim] -
  79           array_.GetDimension(zeroBasedDim).LowerBound() + 1;
  80     } else {
  81       for (int j{0}; j < argRank_; ++j) {
  82         p[j] = extremumLoc_[j] - array_.GetDimension(j).LowerBound() + 1;
  83       }
  84     }
  85   }
  86   template <typename IGNORED> bool AccumulateAt(const SubscriptValue at[]) {
  87     const auto &value{*array_.Element<Type>(at)};
  88     if (!previous_ || compare_(value, *previous_)) {
  89       previous_ = &value;
  90       for (int j{0}; j < argRank_; ++j) {
  91         extremumLoc_[j] = at[j];
  92       }
  93     }
  94     return true;
  95   }
  96
  97 private:
  98   const Descriptor &array_;
  99   int argRank_;
 100   SubscriptValue extremumLoc_[maxRank];
 101   const Type *previous_{nullptr};
 102   COMPARE compare_;
 103 };
 104
 105 template <typename ACCUMULATOR, typename CPPTYPE>
 106 static void LocationHelper(const char *intrinsic, Descriptor &result,
 107     const Descriptor &x, int kind, const Descriptor *mask,
 108     Terminator &terminator) {
 109   ACCUMULATOR accumulator{x};
 110   DoTotalReduction<CPPTYPE>(x, 0, mask, accumulator, intrinsic, terminator);
 111   ApplyIntegerKind<LocationResultHelper<ACCUMULATOR>::template Functor, void>(
 112       kind, terminator, accumulator, result);
 113 }
 114
 115 template <TypeCategory CAT, int KIND, bool IS_MAX,
 116     template <typename, bool, bool> class COMPARE>
 117 inline void DoMaxOrMinLoc(const char *intrinsic, Descriptor &result,
 118     const Descriptor &x, int kind, const char *source, int line,
 119     const Descriptor *mask, bool back) {
 120   using CppType = CppTypeFor<CAT, KIND>;
 121   Terminator terminator{source, line};
 122   if (back) {
 123     LocationHelper<ExtremumLocAccumulator<COMPARE<CppType, IS_MAX, true>>,
 124         CppType>(intrinsic, result, x, kind, mask, terminator);
 125   } else {
 126     LocationHelper<ExtremumLocAccumulator<COMPARE<CppType, IS_MAX, false>>,
 127         CppType>(intrinsic, result, x, kind, mask, terminator);
 128   }
 129 }
 130
 131 template <bool IS_MAX> struct CharacterMaxOrMinLocHelper {
 132   template <int KIND> struct Functor {
 133     void operator()(const char *intrinsic, Descriptor &result,
 134         const Descriptor &x, int kind, const char *source, int line,
 135         const Descriptor *mask, bool back) const {
 136       DoMaxOrMinLoc<TypeCategory::Character, KIND, IS_MAX, NumericCompare>(
 137           intrinsic, result, x, kind, source, line, mask, back);
 138     }
 139   };
 140 };
 141
 142 template <bool IS_MAX>
 143 inline void CharacterMaxOrMinLoc(const char *intrinsic, Descriptor &result,
 144     const Descriptor &x, int kind, const char *source, int line,
 145     const Descriptor *mask, bool back) {
 146   int rank{x.rank()};
 147   SubscriptValue extent[1]{rank};
 148   result.Establish(TypeCategory::Integer, kind, nullptr, 1, extent,
 149       CFI_attribute_allocatable);
 150   result.GetDimension(0).SetBounds(1, extent[0]);
 151   Terminator terminator{source, line};
 152   if (int stat{result.Allocate()}) {
 153     terminator.Crash(
 154         "%s: could not allocate memory for result; STAT=%d", intrinsic, stat);
 155   }
 156   CheckIntegerKind(terminator, kind, intrinsic);
 157   auto catKind{x.type().GetCategoryAndKind()};
 158   RUNTIME_CHECK(terminator, catKind.has_value());
 159   switch (catKind->first) {
 160   case TypeCategory::Character:
 161     ApplyCharacterKind<CharacterMaxOrMinLocHelper<IS_MAX>::template Functor,
 162         void>(catKind->second, terminator, intrinsic, result, x, kind, source,
 163         line, mask, back);
 164     break;
 165   default:
 166     terminator.Crash(
 167         "%s: bad data type code (%d) for array", intrinsic, x.type().raw());
 168   }
 169 }
 170
 171 template <TypeCategory CAT, int KIND, bool IS_MAXVAL>
 172 inline void TotalNumericMaxOrMinLoc(const char *intrinsic, Descriptor &result,
 173     const Descriptor &x, int kind, const char *source, int line,
 174     const Descriptor *mask, bool back) {
 175   int rank{x.rank()};
 176   SubscriptValue extent[1]{rank};
 177   result.Establish(TypeCategory::Integer, kind, nullptr, 1, extent,
 178       CFI_attribute_allocatable);
 179   result.GetDimension(0).SetBounds(1, extent[0]);
 180   Terminator terminator{source, line};
 181   if (int stat{result.Allocate()}) {
 182     terminator.Crash(
 183         "%s: could not allocate memory for result; STAT=%d", intrinsic, stat);
 184   }
 185   CheckIntegerKind(terminator, kind, intrinsic);
 186   RUNTIME_CHECK(terminator, TypeCode(CAT, KIND) == x.type());
 187   DoMaxOrMinLoc<CAT, KIND, IS_MAXVAL, NumericCompare>(
 188       intrinsic, result, x, kind, source, line, mask, back);
 189 }
 190
 191 extern "C" {
 192 void RTNAME(MaxlocCharacter)(Descriptor &result, const Descriptor &x, int kind,
 193     const char *source, int line, const Descriptor *mask, bool back) {
 194   CharacterMaxOrMinLoc<true>(
 195       "MAXLOC", result, x, kind, source, line, mask, back);
 196 }
 197 void RTNAME(MaxlocInteger1)(Descriptor &result, const Descriptor &x, int kind,
 198     const char *source, int line, const Descriptor *mask, bool back) {
 199   TotalNumericMaxOrMinLoc<TypeCategory::Integer, 1, true>(
 200       "MAXLOC", result, x, kind, source, line, mask, back);
 201 }
 202 void RTNAME(MaxlocInteger2)(Descriptor &result, const Descriptor &x, int kind,
 203     const char *source, int line, const Descriptor *mask, bool back) {
 204   TotalNumericMaxOrMinLoc<TypeCategory::Integer, 2, true>(
 205       "MAXLOC", result, x, kind, source, line, mask, back);
 206 }
 207 void RTNAME(MaxlocInteger4)(Descriptor &result, const Descriptor &x, int kind,
 208     const char *source, int line, const Descriptor *mask, bool back) {
 209   TotalNumericMaxOrMinLoc<TypeCategory::Integer, 4, true>(
 210       "MAXLOC", result, x, kind, source, line, mask, back);
 211 }
 212 void RTNAME(MaxlocInteger8)(Descriptor &result, const Descriptor &x, int kind,
 213     const char *source, int line, const Descriptor *mask, bool back) {
 214   TotalNumericMaxOrMinLoc<TypeCategory::Integer, 8, true>(
 215       "MAXLOC", result, x, kind, source, line, mask, back);
 216 }
 217 #ifdef __SIZEOF_INT128__
 218 void RTNAME(MaxlocInteger16)(Descriptor &result, const Descriptor &x, int kind,
 219     const char *source, int line, const Descriptor *mask, bool back) {
 220   TotalNumericMaxOrMinLoc<TypeCategory::Integer, 16, true>(
 221       "MAXLOC", result, x, kind, source, line, mask, back);
 222 }
 223 #endif
 224 void RTNAME(MaxlocReal4)(Descriptor &result, const Descriptor &x, int kind,
 225     const char *source, int line, const Descriptor *mask, bool back) {
 226   TotalNumericMaxOrMinLoc<TypeCategory::Real, 4, true>(
 227       "MAXLOC", result, x, kind, source, line, mask, back);
 228 }
 229 void RTNAME(MaxlocReal8)(Descriptor &result, const Descriptor &x, int kind,
 230     const char *source, int line, const Descriptor *mask, bool back) {
 231   TotalNumericMaxOrMinLoc<TypeCategory::Real, 8, true>(
 232       "MAXLOC", result, x, kind, source, line, mask, back);
 233 }
 234 #if LDBL_MANT_DIG == 64
 235 void RTNAME(MaxlocReal10)(Descriptor &result, const Descriptor &x, int kind,
 236     const char *source, int line, const Descriptor *mask, bool back) {
 237   TotalNumericMaxOrMinLoc<TypeCategory::Real, 10, true>(
 238       "MAXLOC", result, x, kind, source, line, mask, back);
 239 }
 240 #endif
 241 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 242 void RTNAME(MaxlocReal16)(Descriptor &result, const Descriptor &x, int kind,
 243     const char *source, int line, const Descriptor *mask, bool back) {
 244   TotalNumericMaxOrMinLoc<TypeCategory::Real, 16, true>(
 245       "MAXLOC", result, x, kind, source, line, mask, back);
 246 }
 247 #endif
 248 void RTNAME(MinlocCharacter)(Descriptor &result, const Descriptor &x, int kind,
 249     const char *source, int line, const Descriptor *mask, bool back) {
 250   CharacterMaxOrMinLoc<false>(
 251       "MINLOC", result, x, kind, source, line, mask, back);
 252 }
 253 void RTNAME(MinlocInteger1)(Descriptor &result, const Descriptor &x, int kind,
 254     const char *source, int line, const Descriptor *mask, bool back) {
 255   TotalNumericMaxOrMinLoc<TypeCategory::Integer, 1, false>(
 256       "MINLOC", result, x, kind, source, line, mask, back);
 257 }
 258 void RTNAME(MinlocInteger2)(Descriptor &result, const Descriptor &x, int kind,
 259     const char *source, int line, const Descriptor *mask, bool back) {
 260   TotalNumericMaxOrMinLoc<TypeCategory::Integer, 2, false>(
 261       "MINLOC", result, x, kind, source, line, mask, back);
 262 }
 263 void RTNAME(MinlocInteger4)(Descriptor &result, const Descriptor &x, int kind,
 264     const char *source, int line, const Descriptor *mask, bool back) {
 265   TotalNumericMaxOrMinLoc<TypeCategory::Integer, 4, false>(
 266       "MINLOC", result, x, kind, source, line, mask, back);
 267 }
 268 void RTNAME(MinlocInteger8)(Descriptor &result, const Descriptor &x, int kind,
 269     const char *source, int line, const Descriptor *mask, bool back) {
 270   TotalNumericMaxOrMinLoc<TypeCategory::Integer, 8, false>(
 271       "MINLOC", result, x, kind, source, line, mask, back);
 272 }
 273 #ifdef __SIZEOF_INT128__
 274 void RTNAME(MinlocInteger16)(Descriptor &result, const Descriptor &x, int kind,
 275     const char *source, int line, const Descriptor *mask, bool back) {
 276   TotalNumericMaxOrMinLoc<TypeCategory::Integer, 16, false>(
 277       "MINLOC", result, x, kind, source, line, mask, back);
 278 }
 279 #endif
 280 void RTNAME(MinlocReal4)(Descriptor &result, const Descriptor &x, int kind,
 281     const char *source, int line, const Descriptor *mask, bool back) {
 282   TotalNumericMaxOrMinLoc<TypeCategory::Real, 4, false>(
 283       "MINLOC", result, x, kind, source, line, mask, back);
 284 }
 285 void RTNAME(MinlocReal8)(Descriptor &result, const Descriptor &x, int kind,
 286     const char *source, int line, const Descriptor *mask, bool back) {
 287   TotalNumericMaxOrMinLoc<TypeCategory::Real, 8, false>(
 288       "MINLOC", result, x, kind, source, line, mask, back);
 289 }
 290 #if LDBL_MANT_DIG == 64
 291 void RTNAME(MinlocReal10)(Descriptor &result, const Descriptor &x, int kind,
 292     const char *source, int line, const Descriptor *mask, bool back) {
 293   TotalNumericMaxOrMinLoc<TypeCategory::Real, 10, false>(
 294       "MINLOC", result, x, kind, source, line, mask, back);
 295 }
 296 #endif
 297 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 298 void RTNAME(MinlocReal16)(Descriptor &result, const Descriptor &x, int kind,
 299     const char *source, int line, const Descriptor *mask, bool back) {
 300   TotalNumericMaxOrMinLoc<TypeCategory::Real, 16, false>(
 301       "MINLOC", result, x, kind, source, line, mask, back);
 302 }
 303 #endif
 304 } // extern "C"
 305
 306 // MAXLOC/MINLOC with DIM=
 307
 308 template <TypeCategory CAT, int KIND, bool IS_MAX,
 309     template <typename, bool, bool> class COMPARE, bool BACK>
 310 static void DoPartialMaxOrMinLocDirection(const char *intrinsic,
 311     Descriptor &result, const Descriptor &x, int kind, int dim,
 312     const Descriptor *mask, Terminator &terminator) {
 313   using CppType = CppTypeFor<CAT, KIND>;
 314   using Accumulator = ExtremumLocAccumulator<COMPARE<CppType, IS_MAX, BACK>>;
 315   Accumulator accumulator{x};
 316   ApplyIntegerKind<PartialLocationHelper<Accumulator>::template Functor, void>(
 317       kind, terminator, result, x, dim, mask, terminator, intrinsic,
 318       accumulator);
 319 }
 320
 321 template <TypeCategory CAT, int KIND, bool IS_MAX,
 322     template <typename, bool, bool> class COMPARE>
 323 inline void DoPartialMaxOrMinLoc(const char *intrinsic, Descriptor &result,
 324     const Descriptor &x, int kind, int dim, const Descriptor *mask, bool back,
 325     Terminator &terminator) {
 326   if (back) {
 327     DoPartialMaxOrMinLocDirection<CAT, KIND, IS_MAX, COMPARE, true>(
 328         intrinsic, result, x, kind, dim, mask, terminator);
 329   } else {
 330     DoPartialMaxOrMinLocDirection<CAT, KIND, IS_MAX, COMPARE, false>(
 331         intrinsic, result, x, kind, dim, mask, terminator);
 332   }
 333 }
 334
 335 template <TypeCategory CAT, bool IS_MAX,
 336     template <typename, bool, bool> class COMPARE>
 337 struct DoPartialMaxOrMinLocHelper {
 338   template <int KIND> struct Functor {
 339     void operator()(const char *intrinsic, Descriptor &result,
 340         const Descriptor &x, int kind, int dim, const Descriptor *mask,
 341         bool back, Terminator &terminator) const {
 342       DoPartialMaxOrMinLoc<CAT, KIND, IS_MAX, COMPARE>(
 343           intrinsic, result, x, kind, dim, mask, back, terminator);
 344     }
 345   };
 346 };
 347
 348 template <bool IS_MAX>
 349 inline void TypedPartialMaxOrMinLoc(const char *intrinsic, Descriptor &result,
 350     const Descriptor &x, int kind, int dim, const char *source, int line,
 351     const Descriptor *mask, bool back) {
 352   Terminator terminator{source, line};
 353   CheckIntegerKind(terminator, kind, intrinsic);
 354   auto catKind{x.type().GetCategoryAndKind()};
 355   RUNTIME_CHECK(terminator, catKind.has_value());
 356   const Descriptor *maskToUse{mask};
 357   SubscriptValue maskAt[maxRank]; // contents unused
 358   if (mask && mask->rank() == 0) {
 359     if (IsLogicalElementTrue(*mask, maskAt)) {
 360       // A scalar MASK that's .TRUE.  In this case, just get rid of the MASK.
 361       maskToUse = nullptr;
 362     } else {
 363       // For scalar MASK arguments that are .FALSE., return all zeroes
 364
 365       // Element size of the destination descriptor is the size
 366       // of {TypeCategory::Integer, kind}.
 367       CreatePartialReductionResult(result, x,
 368           Descriptor::BytesFor(TypeCategory::Integer, kind), dim, terminator,
 369           intrinsic, TypeCode{TypeCategory::Integer, kind});
 370       std::memset(
 371           result.OffsetElement(), 0, result.Elements() * result.ElementBytes());
 372       return;
 373     }
 374   }
 375   switch (catKind->first) {
 376   case TypeCategory::Integer:
 377     ApplyIntegerKind<DoPartialMaxOrMinLocHelper<TypeCategory::Integer, IS_MAX,
 378                          NumericCompare>::template Functor,
 379         void>(catKind->second, terminator, intrinsic, result, x, kind, dim,
 380         maskToUse, back, terminator);
 381     break;
 382   case TypeCategory::Real:
 383     ApplyFloatingPointKind<DoPartialMaxOrMinLocHelper<TypeCategory::Real,
 384                                IS_MAX, NumericCompare>::template Functor,
 385         void>(catKind->second, terminator, intrinsic, result, x, kind, dim,
 386         maskToUse, back, terminator);
 387     break;
 388   case TypeCategory::Character:
 389     ApplyCharacterKind<DoPartialMaxOrMinLocHelper<TypeCategory::Character,
 390                            IS_MAX, CharacterCompare>::template Functor,
 391         void>(catKind->second, terminator, intrinsic, result, x, kind, dim,
 392         maskToUse, back, terminator);
 393     break;
 394   default:
 395     terminator.Crash(
 396         "%s: bad data type code (%d) for array", intrinsic, x.type().raw());
 397   }
 398 }
 399
 400 extern "C" {
 401 void RTNAME(MaxlocDim)(Descriptor &result, const Descriptor &x, int kind,
 402     int dim, const char *source, int line, const Descriptor *mask, bool back) {
 403   TypedPartialMaxOrMinLoc<true>(
 404       "MAXLOC", result, x, kind, dim, source, line, mask, back);
 405 }
 406 void RTNAME(MinlocDim)(Descriptor &result, const Descriptor &x, int kind,
 407     int dim, const char *source, int line, const Descriptor *mask, bool back) {
 408   TypedPartialMaxOrMinLoc<false>(
 409       "MINLOC", result, x, kind, dim, source, line, mask, back);
 410 }
 411 } // extern "C"
 412
 413 // MAXVAL and MINVAL
 414
 415 template <TypeCategory CAT, int KIND, bool IS_MAXVAL, typename Enable = void>
 416 struct MaxOrMinIdentity {
 417   using Type = CppTypeFor<CAT, KIND>;
 418   static constexpr Type Value() {
 419     return IS_MAXVAL ? std::numeric_limits<Type>::lowest()
 420                      : std::numeric_limits<Type>::max();
 421   }
 422 };
 423
 424 // std::numeric_limits<> may not know int128_t
 425 template <bool IS_MAXVAL>
 426 struct MaxOrMinIdentity<TypeCategory::Integer, 16, IS_MAXVAL> {
 427   using Type = CppTypeFor<TypeCategory::Integer, 16>;
 428   static constexpr Type Value() {
 429     return IS_MAXVAL ? Type{1} << 127 : ~Type{0} >> 1;
 430   }
 431 };
 432
 433 #if HAS_FLOAT128
 434 // std::numeric_limits<> may not support __float128.
 435 //
 436 // Usage of GCC quadmath.h's FLT128_MAX is complicated by the fact that
 437 // even GCC complains about 'Q' literal suffix under -Wpedantic.
 438 // We just recreate FLT128_MAX ourselves.
 439 //
 440 // This specialization must engage only when
 441 // CppTypeFor<TypeCategory::Real, 16> is __float128.
 442 template <bool IS_MAXVAL>
 443 struct MaxOrMinIdentity<TypeCategory::Real, 16, IS_MAXVAL,
 444     typename std::enable_if_t<
 445         std::is_same_v<CppTypeFor<TypeCategory::Real, 16>, __float128>>> {
 446   using Type = __float128;
 447   static Type Value() {
 448     // Create a buffer to store binary representation of __float128 constant.
 449     constexpr std::size_t alignment =
 450         std::max(alignof(Type), alignof(std::uint64_t));
 451     alignas(alignment) char data[sizeof(Type)];
 452
 453     // First, verify that our interpretation of __float128 format is correct,
 454     // e.g. by checking at least one known constant.
 455     *reinterpret_cast<Type *>(data) = Type(1.0);
 456     if (*reinterpret_cast<std::uint64_t *>(data) != 0 ||
 457         *(reinterpret_cast<std::uint64_t *>(data) + 1) != 0x3FFF000000000000) {
 458       Terminator terminator{__FILE__, __LINE__};
 459       terminator.Crash("not yet implemented: no full support for __float128");
 460     }
 461
 462     // Recreate FLT128_MAX.
 463     *reinterpret_cast<std::uint64_t *>(data) = 0xFFFFFFFFFFFFFFFF;
 464     *(reinterpret_cast<std::uint64_t *>(data) + 1) = 0x7FFEFFFFFFFFFFFF;
 465     Type max = *reinterpret_cast<Type *>(data);
 466     return IS_MAXVAL ? -max : max;
 467   }
 468 };
 469 #endif // HAS_FLOAT128
 470
 471 template <TypeCategory CAT, int KIND, bool IS_MAXVAL>
 472 class NumericExtremumAccumulator {
 473 public:
 474   using Type = CppTypeFor<CAT, KIND>;
 475   explicit NumericExtremumAccumulator(const Descriptor &array)
 476       : array_{array} {}
 477   void Reinitialize() {
 478     extremum_ = MaxOrMinIdentity<CAT, KIND, IS_MAXVAL>::Value();
 479   }
 480   template <typename A> void GetResult(A *p, int /*zeroBasedDim*/ = -1) const {
 481     *p = extremum_;
 482   }
 483   bool Accumulate(Type x) {
 484     if constexpr (IS_MAXVAL) {
 485       if (x > extremum_) {
 486         extremum_ = x;
 487       }
 488     } else if (x < extremum_) {
 489       extremum_ = x;
 490     }
 491     return true;
 492   }
 493   template <typename A> bool AccumulateAt(const SubscriptValue at[]) {
 494     return Accumulate(*array_.Element<A>(at));
 495   }
 496
 497 private:
 498   const Descriptor &array_;
 499   Type extremum_{MaxOrMinIdentity<CAT, KIND, IS_MAXVAL>::Value()};
 500 };
 501
 502 template <TypeCategory CAT, int KIND, bool IS_MAXVAL>
 503 inline CppTypeFor<CAT, KIND> TotalNumericMaxOrMin(const Descriptor &x,
 504     const char *source, int line, int dim, const Descriptor *mask,
 505     const char *intrinsic) {
 506   return GetTotalReduction<CAT, KIND>(x, source, line, dim, mask,
 507       NumericExtremumAccumulator<CAT, KIND, IS_MAXVAL>{x}, intrinsic);
 508 }
 509
 510 template <TypeCategory CAT, int KIND, typename ACCUMULATOR>
 511 static void DoMaxMinNorm2(Descriptor &result, const Descriptor &x, int dim,
 512     const Descriptor *mask, const char *intrinsic, Terminator &terminator) {
 513   using Type = CppTypeFor<CAT, KIND>;
 514   ACCUMULATOR accumulator{x};
 515   if (dim == 0 || x.rank() == 1) {
 516     // Total reduction
 517
 518     // Element size of the destination descriptor is the same
 519     // as the element size of the source.
 520     result.Establish(x.type(), x.ElementBytes(), nullptr, 0, nullptr,
 521         CFI_attribute_allocatable);
 522     if (int stat{result.Allocate()}) {
 523       terminator.Crash(
 524           "%s: could not allocate memory for result; STAT=%d", intrinsic, stat);
 525     }
 526     DoTotalReduction<Type>(x, dim, mask, accumulator, intrinsic, terminator);
 527     accumulator.GetResult(result.OffsetElement<Type>());
 528   } else {
 529     // Partial reduction
 530
 531     // Element size of the destination descriptor is the same
 532     // as the element size of the source.
 533     PartialReduction<ACCUMULATOR, CAT, KIND>(result, x, x.ElementBytes(), dim,
 534         mask, terminator, intrinsic, accumulator);
 535   }
 536 }
 537
 538 template <TypeCategory CAT, bool IS_MAXVAL> struct MaxOrMinHelper {
 539   template <int KIND> struct Functor {
 540     void operator()(Descriptor &result, const Descriptor &x, int dim,
 541         const Descriptor *mask, const char *intrinsic,
 542         Terminator &terminator) const {
 543       DoMaxMinNorm2<CAT, KIND,
 544           NumericExtremumAccumulator<CAT, KIND, IS_MAXVAL>>(
 545           result, x, dim, mask, intrinsic, terminator);
 546     }
 547   };
 548 };
 549
 550 template <bool IS_MAXVAL>
 551 inline void NumericMaxOrMin(Descriptor &result, const Descriptor &x, int dim,
 552     const char *source, int line, const Descriptor *mask,
 553     const char *intrinsic) {
 554   Terminator terminator{source, line};
 555   auto type{x.type().GetCategoryAndKind()};
 556   RUNTIME_CHECK(terminator, type);
 557   switch (type->first) {
 558   case TypeCategory::Integer:
 559     ApplyIntegerKind<
 560         MaxOrMinHelper<TypeCategory::Integer, IS_MAXVAL>::template Functor,
 561         void>(
 562         type->second, terminator, result, x, dim, mask, intrinsic, terminator);
 563     break;
 564   case TypeCategory::Real:
 565     ApplyFloatingPointKind<
 566         MaxOrMinHelper<TypeCategory::Real, IS_MAXVAL>::template Functor, void>(
 567         type->second, terminator, result, x, dim, mask, intrinsic, terminator);
 568     break;
 569   default:
 570     terminator.Crash("%s: bad type code %d", intrinsic, x.type().raw());
 571   }
 572 }
 573
 574 template <int KIND, bool IS_MAXVAL> class CharacterExtremumAccumulator {
 575 public:
 576   using Type = CppTypeFor<TypeCategory::Character, KIND>;
 577   explicit CharacterExtremumAccumulator(const Descriptor &array)
 578       : array_{array}, charLen_{array_.ElementBytes() / KIND} {}
 579   void Reinitialize() { extremum_ = nullptr; }
 580   template <typename A> void GetResult(A *p, int /*zeroBasedDim*/ = -1) const {
 581     static_assert(std::is_same_v<A, Type>);
 582     std::size_t byteSize{array_.ElementBytes()};
 583     if (extremum_) {
 584       std::memcpy(p, extremum_, byteSize);
 585     } else {
 586       // Empty array; fill with character 0 for MAXVAL.
 587       // For MINVAL, fill with 127 if ASCII as required
 588       // by the standard, otherwise set all of the bits.
 589       std::memset(p, IS_MAXVAL ? 0 : KIND == 1 ? 127 : 255, byteSize);
 590     }
 591   }
 592   bool Accumulate(const Type *x) {
 593     if (!extremum_) {
 594       extremum_ = x;
 595     } else {
 596       int cmp{CharacterScalarCompare(x, extremum_, charLen_, charLen_)};
 597       if (IS_MAXVAL == (cmp > 0)) {
 598         extremum_ = x;
 599       }
 600     }
 601     return true;
 602   }
 603   template <typename A> bool AccumulateAt(const SubscriptValue at[]) {
 604     return Accumulate(array_.Element<A>(at));
 605   }
 606
 607 private:
 608   const Descriptor &array_;
 609   std::size_t charLen_;
 610   const Type *extremum_{nullptr};
 611 };
 612
 613 template <bool IS_MAXVAL> struct CharacterMaxOrMinHelper {
 614   template <int KIND> struct Functor {
 615     void operator()(Descriptor &result, const Descriptor &x, int dim,
 616         const Descriptor *mask, const char *intrinsic,
 617         Terminator &terminator) const {
 618       DoMaxMinNorm2<TypeCategory::Character, KIND,
 619           CharacterExtremumAccumulator<KIND, IS_MAXVAL>>(
 620           result, x, dim, mask, intrinsic, terminator);
 621     }
 622   };
 623 };
 624
 625 template <bool IS_MAXVAL>
 626 inline void CharacterMaxOrMin(Descriptor &result, const Descriptor &x, int dim,
 627     const char *source, int line, const Descriptor *mask,
 628     const char *intrinsic) {
 629   Terminator terminator{source, line};
 630   auto type{x.type().GetCategoryAndKind()};
 631   RUNTIME_CHECK(terminator, type && type->first == TypeCategory::Character);
 632   ApplyCharacterKind<CharacterMaxOrMinHelper<IS_MAXVAL>::template Functor,
 633       void>(
 634       type->second, terminator, result, x, dim, mask, intrinsic, terminator);
 635 }
 636
 637 extern "C" {
 638 CppTypeFor<TypeCategory::Integer, 1> RTNAME(MaxvalInteger1)(const Descriptor &x,
 639     const char *source, int line, int dim, const Descriptor *mask) {
 640   return TotalNumericMaxOrMin<TypeCategory::Integer, 1, true>(
 641       x, source, line, dim, mask, "MAXVAL");
 642 }
 643 CppTypeFor<TypeCategory::Integer, 2> RTNAME(MaxvalInteger2)(const Descriptor &x,
 644     const char *source, int line, int dim, const Descriptor *mask) {
 645   return TotalNumericMaxOrMin<TypeCategory::Integer, 2, true>(
 646       x, source, line, dim, mask, "MAXVAL");
 647 }
 648 CppTypeFor<TypeCategory::Integer, 4> RTNAME(MaxvalInteger4)(const Descriptor &x,
 649     const char *source, int line, int dim, const Descriptor *mask) {
 650   return TotalNumericMaxOrMin<TypeCategory::Integer, 4, true>(
 651       x, source, line, dim, mask, "MAXVAL");
 652 }
 653 CppTypeFor<TypeCategory::Integer, 8> RTNAME(MaxvalInteger8)(const Descriptor &x,
 654     const char *source, int line, int dim, const Descriptor *mask) {
 655   return TotalNumericMaxOrMin<TypeCategory::Integer, 8, true>(
 656       x, source, line, dim, mask, "MAXVAL");
 657 }
 658 #ifdef __SIZEOF_INT128__
 659 CppTypeFor<TypeCategory::Integer, 16> RTNAME(MaxvalInteger16)(
 660     const Descriptor &x, const char *source, int line, int dim,
 661     const Descriptor *mask) {
 662   return TotalNumericMaxOrMin<TypeCategory::Integer, 16, true>(
 663       x, source, line, dim, mask, "MAXVAL");
 664 }
 665 #endif
 666
 667 // TODO: REAL(2 & 3)
 668 CppTypeFor<TypeCategory::Real, 4> RTNAME(MaxvalReal4)(const Descriptor &x,
 669     const char *source, int line, int dim, const Descriptor *mask) {
 670   return TotalNumericMaxOrMin<TypeCategory::Real, 4, true>(
 671       x, source, line, dim, mask, "MAXVAL");
 672 }
 673 CppTypeFor<TypeCategory::Real, 8> RTNAME(MaxvalReal8)(const Descriptor &x,
 674     const char *source, int line, int dim, const Descriptor *mask) {
 675   return TotalNumericMaxOrMin<TypeCategory::Real, 8, true>(
 676       x, source, line, dim, mask, "MAXVAL");
 677 }
 678 #if LDBL_MANT_DIG == 64
 679 CppTypeFor<TypeCategory::Real, 10> RTNAME(MaxvalReal10)(const Descriptor &x,
 680     const char *source, int line, int dim, const Descriptor *mask) {
 681   return TotalNumericMaxOrMin<TypeCategory::Real, 10, true>(
 682       x, source, line, dim, mask, "MAXVAL");
 683 }
 684 #endif
 685 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 686 CppTypeFor<TypeCategory::Real, 16> RTNAME(MaxvalReal16)(const Descriptor &x,
 687     const char *source, int line, int dim, const Descriptor *mask) {
 688   return TotalNumericMaxOrMin<TypeCategory::Real, 16, true>(
 689       x, source, line, dim, mask, "MAXVAL");
 690 }
 691 #endif
 692
 693 void RTNAME(MaxvalCharacter)(Descriptor &result, const Descriptor &x,
 694     const char *source, int line, const Descriptor *mask) {
 695   CharacterMaxOrMin<true>(result, x, 0, source, line, mask, "MAXVAL");
 696 }
 697
 698 CppTypeFor<TypeCategory::Integer, 1> RTNAME(MinvalInteger1)(const Descriptor &x,
 699     const char *source, int line, int dim, const Descriptor *mask) {
 700   return TotalNumericMaxOrMin<TypeCategory::Integer, 1, false>(
 701       x, source, line, dim, mask, "MINVAL");
 702 }
 703 CppTypeFor<TypeCategory::Integer, 2> RTNAME(MinvalInteger2)(const Descriptor &x,
 704     const char *source, int line, int dim, const Descriptor *mask) {
 705   return TotalNumericMaxOrMin<TypeCategory::Integer, 2, false>(
 706       x, source, line, dim, mask, "MINVAL");
 707 }
 708 CppTypeFor<TypeCategory::Integer, 4> RTNAME(MinvalInteger4)(const Descriptor &x,
 709     const char *source, int line, int dim, const Descriptor *mask) {
 710   return TotalNumericMaxOrMin<TypeCategory::Integer, 4, false>(
 711       x, source, line, dim, mask, "MINVAL");
 712 }
 713 CppTypeFor<TypeCategory::Integer, 8> RTNAME(MinvalInteger8)(const Descriptor &x,
 714     const char *source, int line, int dim, const Descriptor *mask) {
 715   return TotalNumericMaxOrMin<TypeCategory::Integer, 8, false>(
 716       x, source, line, dim, mask, "MINVAL");
 717 }
 718 #ifdef __SIZEOF_INT128__
 719 CppTypeFor<TypeCategory::Integer, 16> RTNAME(MinvalInteger16)(
 720     const Descriptor &x, const char *source, int line, int dim,
 721     const Descriptor *mask) {
 722   return TotalNumericMaxOrMin<TypeCategory::Integer, 16, false>(
 723       x, source, line, dim, mask, "MINVAL");
 724 }
 725 #endif
 726
 727 // TODO: REAL(2 & 3)
 728 CppTypeFor<TypeCategory::Real, 4> RTNAME(MinvalReal4)(const Descriptor &x,
 729     const char *source, int line, int dim, const Descriptor *mask) {
 730   return TotalNumericMaxOrMin<TypeCategory::Real, 4, false>(
 731       x, source, line, dim, mask, "MINVAL");
 732 }
 733 CppTypeFor<TypeCategory::Real, 8> RTNAME(MinvalReal8)(const Descriptor &x,
 734     const char *source, int line, int dim, const Descriptor *mask) {
 735   return TotalNumericMaxOrMin<TypeCategory::Real, 8, false>(
 736       x, source, line, dim, mask, "MINVAL");
 737 }
 738 #if LDBL_MANT_DIG == 64
 739 CppTypeFor<TypeCategory::Real, 10> RTNAME(MinvalReal10)(const Descriptor &x,
 740     const char *source, int line, int dim, const Descriptor *mask) {
 741   return TotalNumericMaxOrMin<TypeCategory::Real, 10, false>(
 742       x, source, line, dim, mask, "MINVAL");
 743 }
 744 #endif
 745 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 746 CppTypeFor<TypeCategory::Real, 16> RTNAME(MinvalReal16)(const Descriptor &x,
 747     const char *source, int line, int dim, const Descriptor *mask) {
 748   return TotalNumericMaxOrMin<TypeCategory::Real, 16, false>(
 749       x, source, line, dim, mask, "MINVAL");
 750 }
 751 #endif
 752
 753 void RTNAME(MinvalCharacter)(Descriptor &result, const Descriptor &x,
 754     const char *source, int line, const Descriptor *mask) {
 755   CharacterMaxOrMin<false>(result, x, 0, source, line, mask, "MINVAL");
 756 }
 757
 758 void RTNAME(MaxvalDim)(Descriptor &result, const Descriptor &x, int dim,
 759     const char *source, int line, const Descriptor *mask) {
 760   if (x.type().IsCharacter()) {
 761     CharacterMaxOrMin<true>(result, x, dim, source, line, mask, "MAXVAL");
 762   } else {
 763     NumericMaxOrMin<true>(result, x, dim, source, line, mask, "MAXVAL");
 764   }
 765 }
 766 void RTNAME(MinvalDim)(Descriptor &result, const Descriptor &x, int dim,
 767     const char *source, int line, const Descriptor *mask) {
 768   if (x.type().IsCharacter()) {
 769     CharacterMaxOrMin<false>(result, x, dim, source, line, mask, "MINVAL");
 770   } else {
 771     NumericMaxOrMin<false>(result, x, dim, source, line, mask, "MINVAL");
 772   }
 773 }
 774 } // extern "C"
 775
 776 // NORM2
 777
 778 template <int KIND> class Norm2Accumulator {
 779 public:
 780   using Type = CppTypeFor<TypeCategory::Real, KIND>;
 781   // Use at least double precision for accumulators.
 782   // Don't use __float128, it doesn't work with abs() or sqrt() yet.
 783   static constexpr int largestLDKind {
 784 #if LDBL_MANT_DIG == 113
 785     16
 786 #elif LDBL_MANT_DIG == 64
 787     10
 788 #else
 789     8
 790 #endif
 791   };
 792   using AccumType =
 793       CppTypeFor<TypeCategory::Real, std::clamp(KIND, 8, largestLDKind)>;
 794   explicit Norm2Accumulator(const Descriptor &array) : array_{array} {}
 795   void Reinitialize() { max_ = sum_ = 0; }
 796   template <typename A> void GetResult(A *p, int /*zeroBasedDim*/ = -1) const {
 797     // m * sqrt(1 + sum((others(:)/m)**2))
 798     *p = static_cast<Type>(max_ * std::sqrt(1 + sum_));
 799   }
 800   bool Accumulate(Type x) {
 801     auto absX{std::abs(static_cast<AccumType>(x))};
 802     if (!max_) {
 803       max_ = absX;
 804     } else if (absX > max_) {
 805       auto t{max_ / absX}; // < 1.0
 806       auto tsq{t * t};
 807       sum_ *= tsq; // scale sum to reflect change to the max
 808       sum_ += tsq; // include a term for the previous max
 809       max_ = absX;
 810     } else { // absX <= max_
 811       auto t{absX / max_};
 812       sum_ += t * t;
 813     }
 814     return true;
 815   }
 816   template <typename A> bool AccumulateAt(const SubscriptValue at[]) {
 817     return Accumulate(*array_.Element<A>(at));
 818   }
 819
 820 private:
 821   const Descriptor &array_;
 822   AccumType max_{0}; // value (m) with largest magnitude
 823   AccumType sum_{0}; // sum((others(:)/m)**2)
 824 };
 825
 826 template <int KIND> struct Norm2Helper {
 827   void operator()(Descriptor &result, const Descriptor &x, int dim,
 828       const Descriptor *mask, Terminator &terminator) const {
 829     DoMaxMinNorm2<TypeCategory::Real, KIND, Norm2Accumulator<KIND>>(
 830         result, x, dim, mask, "NORM2", terminator);
 831   }
 832 };
 833
 834 extern "C" {
 835 // TODO: REAL(2 & 3)
 836 CppTypeFor<TypeCategory::Real, 4> RTNAME(Norm2_4)(
 837     const Descriptor &x, const char *source, int line, int dim) {
 838   return GetTotalReduction<TypeCategory::Real, 4>(
 839       x, source, line, dim, nullptr, Norm2Accumulator<4>{x}, "NORM2");
 840 }
 841 CppTypeFor<TypeCategory::Real, 8> RTNAME(Norm2_8)(
 842     const Descriptor &x, const char *source, int line, int dim) {
 843   return GetTotalReduction<TypeCategory::Real, 8>(
 844       x, source, line, dim, nullptr, Norm2Accumulator<8>{x}, "NORM2");
 845 }
 846 #if LDBL_MANT_DIG == 64
 847 CppTypeFor<TypeCategory::Real, 10> RTNAME(Norm2_10)(
 848     const Descriptor &x, const char *source, int line, int dim) {
 849   return GetTotalReduction<TypeCategory::Real, 10>(
 850       x, source, line, dim, nullptr, Norm2Accumulator<10>{x}, "NORM2");
 851 }
 852 #endif
 853 #if LDBL_MANT_DIG == 113
 854 CppTypeFor<TypeCategory::Real, 16> RTNAME(Norm2_16)(
 855     const Descriptor &x, const char *source, int line, int dim) {
 856   return GetTotalReduction<TypeCategory::Real, 16>(
 857       x, source, line, dim, nullptr, Norm2Accumulator<16>{x}, "NORM2");
 858 }
 859 #endif
 860
 861 void RTNAME(Norm2Dim)(Descriptor &result, const Descriptor &x, int dim,
 862     const char *source, int line) {
 863   Terminator terminator{source, line};
 864   auto type{x.type().GetCategoryAndKind()};
 865   RUNTIME_CHECK(terminator, type);
 866   if (type->first == TypeCategory::Real) {
 867     ApplyFloatingPointKind<Norm2Helper, void>(
 868         type->second, terminator, result, x, dim, nullptr, terminator);
 869   } else {
 870     terminator.Crash("NORM2: bad type code %d", x.type().raw());
 871   }
 872 }
 873 } // extern "C"
 874 } // namespace Fortran::runtime