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/Runtime/character.h"
15 #include "flang/Runtime/float128.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 <TypeCategory CAT, bool IS_MAX> struct TypedMaxOrMinLocHelper {
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<CAT, 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 TypedMaxOrMinLoc(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::Integer:
161     ApplyIntegerKind<
162         TypedMaxOrMinLocHelper<TypeCategory::Integer, IS_MAX>::template Functor,
163         void>(catKind->second, terminator, intrinsic, result, x, kind, source,
164         line, mask, back);
165     break;
166   case TypeCategory::Real:
167     ApplyFloatingPointKind<
168         TypedMaxOrMinLocHelper<TypeCategory::Real, IS_MAX>::template Functor,
169         void>(catKind->second, terminator, intrinsic, result, x, kind, source,
170         line, mask, back);
171     break;
172   case TypeCategory::Character:
173     ApplyCharacterKind<TypedMaxOrMinLocHelper<TypeCategory::Character,
174                            IS_MAX>::template Functor,
175         void>(catKind->second, terminator, intrinsic, result, x, kind, source,
176         line, mask, back);
177     break;
178   default:
179     terminator.Crash(
180         "%s: bad data type code (%d) for array", intrinsic, x.type().raw());
181   }
182 }
183 
184 extern "C" {
185 void RTNAME(Maxloc)(Descriptor &result, const Descriptor &x, int kind,
186     const char *source, int line, const Descriptor *mask, bool back) {
187   TypedMaxOrMinLoc<true>("MAXLOC", result, x, kind, source, line, mask, back);
188 }
189 void RTNAME(Minloc)(Descriptor &result, const Descriptor &x, int kind,
190     const char *source, int line, const Descriptor *mask, bool back) {
191   TypedMaxOrMinLoc<false>("MINLOC", result, x, kind, source, line, mask, back);
192 }
193 } // extern "C"
194 
195 // MAXLOC/MINLOC with DIM=
196 
197 template <TypeCategory CAT, int KIND, bool IS_MAX,
198     template <typename, bool, bool> class COMPARE, bool BACK>
199 static void DoPartialMaxOrMinLocDirection(const char *intrinsic,
200     Descriptor &result, const Descriptor &x, int kind, int dim,
201     const Descriptor *mask, Terminator &terminator) {
202   using CppType = CppTypeFor<CAT, KIND>;
203   using Accumulator = ExtremumLocAccumulator<COMPARE<CppType, IS_MAX, BACK>>;
204   Accumulator accumulator{x};
205   ApplyIntegerKind<PartialLocationHelper<Accumulator>::template Functor, void>(
206       kind, terminator, result, x, dim, mask, terminator, intrinsic,
207       accumulator);
208 }
209 
210 template <TypeCategory CAT, int KIND, bool IS_MAX,
211     template <typename, bool, bool> class COMPARE>
212 inline void DoPartialMaxOrMinLoc(const char *intrinsic, Descriptor &result,
213     const Descriptor &x, int kind, int dim, const Descriptor *mask, bool back,
214     Terminator &terminator) {
215   if (back) {
216     DoPartialMaxOrMinLocDirection<CAT, KIND, IS_MAX, COMPARE, true>(
217         intrinsic, result, x, kind, dim, mask, terminator);
218   } else {
219     DoPartialMaxOrMinLocDirection<CAT, KIND, IS_MAX, COMPARE, false>(
220         intrinsic, result, x, kind, dim, mask, terminator);
221   }
222 }
223 
224 template <TypeCategory CAT, bool IS_MAX,
225     template <typename, bool, bool> class COMPARE>
226 struct DoPartialMaxOrMinLocHelper {
227   template <int KIND> struct Functor {
228     void operator()(const char *intrinsic, Descriptor &result,
229         const Descriptor &x, int kind, int dim, const Descriptor *mask,
230         bool back, Terminator &terminator) const {
231       DoPartialMaxOrMinLoc<CAT, KIND, IS_MAX, COMPARE>(
232           intrinsic, result, x, kind, dim, mask, back, terminator);
233     }
234   };
235 };
236 
237 template <bool IS_MAX>
238 inline void TypedPartialMaxOrMinLoc(const char *intrinsic, Descriptor &result,
239     const Descriptor &x, int kind, int dim, const char *source, int line,
240     const Descriptor *mask, bool back) {
241   Terminator terminator{source, line};
242   CheckIntegerKind(terminator, kind, intrinsic);
243   auto catKind{x.type().GetCategoryAndKind()};
244   RUNTIME_CHECK(terminator, catKind.has_value());
245   const Descriptor *maskToUse{mask};
246   SubscriptValue maskAt[maxRank]; // contents unused
247   if (mask && mask->rank() == 0) {
248     if (IsLogicalElementTrue(*mask, maskAt)) {
249       // A scalar MASK that's .TRUE.  In this case, just get rid of the MASK.
250       maskToUse = nullptr;
251     } else {
252       // For scalar MASK arguments that are .FALSE., return all zeroes
253       CreatePartialReductionResult(result, x, dim, terminator, intrinsic,
254           TypeCode{TypeCategory::Integer, kind});
255       std::memset(
256           result.OffsetElement(), 0, result.Elements() * result.ElementBytes());
257       return;
258     }
259   }
260   switch (catKind->first) {
261   case TypeCategory::Integer:
262     ApplyIntegerKind<DoPartialMaxOrMinLocHelper<TypeCategory::Integer, IS_MAX,
263                          NumericCompare>::template Functor,
264         void>(catKind->second, terminator, intrinsic, result, x, kind, dim,
265         maskToUse, back, terminator);
266     break;
267   case TypeCategory::Real:
268     ApplyFloatingPointKind<DoPartialMaxOrMinLocHelper<TypeCategory::Real,
269                                IS_MAX, NumericCompare>::template Functor,
270         void>(catKind->second, terminator, intrinsic, result, x, kind, dim,
271         maskToUse, back, terminator);
272     break;
273   case TypeCategory::Character:
274     ApplyCharacterKind<DoPartialMaxOrMinLocHelper<TypeCategory::Character,
275                            IS_MAX, CharacterCompare>::template Functor,
276         void>(catKind->second, terminator, intrinsic, result, x, kind, dim,
277         maskToUse, back, terminator);
278     break;
279   default:
280     terminator.Crash(
281         "%s: bad data type code (%d) for array", intrinsic, x.type().raw());
282   }
283 }
284 
285 extern "C" {
286 void RTNAME(MaxlocDim)(Descriptor &result, const Descriptor &x, int kind,
287     int dim, const char *source, int line, const Descriptor *mask, bool back) {
288   TypedPartialMaxOrMinLoc<true>(
289       "MAXLOC", result, x, kind, dim, source, line, mask, back);
290 }
291 void RTNAME(MinlocDim)(Descriptor &result, const Descriptor &x, int kind,
292     int dim, const char *source, int line, const Descriptor *mask, bool back) {
293   TypedPartialMaxOrMinLoc<false>(
294       "MINLOC", result, x, kind, dim, source, line, mask, back);
295 }
296 } // extern "C"
297 
298 // MAXVAL and MINVAL
299 
300 template <TypeCategory CAT, int KIND, bool IS_MAXVAL> struct MaxOrMinIdentity {
301   using Type = CppTypeFor<CAT, KIND>;
302   static constexpr Type Value() {
303     return IS_MAXVAL ? std::numeric_limits<Type>::lowest()
304                      : std::numeric_limits<Type>::max();
305   }
306 };
307 
308 // std::numeric_limits<> may not know int128_t
309 template <bool IS_MAXVAL>
310 struct MaxOrMinIdentity<TypeCategory::Integer, 16, IS_MAXVAL> {
311   using Type = CppTypeFor<TypeCategory::Integer, 16>;
312   static constexpr Type Value() {
313     return IS_MAXVAL ? Type{1} << 127 : ~Type{0} >> 1;
314   }
315 };
316 
317 template <TypeCategory CAT, int KIND, bool IS_MAXVAL>
318 class NumericExtremumAccumulator {
319 public:
320   using Type = CppTypeFor<CAT, KIND>;
321   explicit NumericExtremumAccumulator(const Descriptor &array)
322       : array_{array} {}
323   void Reinitialize() {
324     extremum_ = MaxOrMinIdentity<CAT, KIND, IS_MAXVAL>::Value();
325   }
326   template <typename A> void GetResult(A *p, int /*zeroBasedDim*/ = -1) const {
327     *p = extremum_;
328   }
329   bool Accumulate(Type x) {
330     if constexpr (IS_MAXVAL) {
331       if (x > extremum_) {
332         extremum_ = x;
333       }
334     } else if (x < extremum_) {
335       extremum_ = x;
336     }
337     return true;
338   }
339   template <typename A> bool AccumulateAt(const SubscriptValue at[]) {
340     return Accumulate(*array_.Element<A>(at));
341   }
342 
343 private:
344   const Descriptor &array_;
345   Type extremum_{MaxOrMinIdentity<CAT, KIND, IS_MAXVAL>::Value()};
346 };
347 
348 template <TypeCategory CAT, int KIND, bool IS_MAXVAL>
349 inline CppTypeFor<CAT, KIND> TotalNumericMaxOrMin(const Descriptor &x,
350     const char *source, int line, int dim, const Descriptor *mask,
351     const char *intrinsic) {
352   return GetTotalReduction<CAT, KIND>(x, source, line, dim, mask,
353       NumericExtremumAccumulator<CAT, KIND, IS_MAXVAL>{x}, intrinsic);
354 }
355 
356 template <TypeCategory CAT, int KIND, typename ACCUMULATOR>
357 static void DoMaxMinNorm2(Descriptor &result, const Descriptor &x, int dim,
358     const Descriptor *mask, const char *intrinsic, Terminator &terminator) {
359   using Type = CppTypeFor<CAT, KIND>;
360   ACCUMULATOR accumulator{x};
361   if (dim == 0 || x.rank() == 1) {
362     // Total reduction
363     result.Establish(x.type(), x.ElementBytes(), nullptr, 0, nullptr,
364         CFI_attribute_allocatable);
365     if (int stat{result.Allocate()}) {
366       terminator.Crash(
367           "%s: could not allocate memory for result; STAT=%d", intrinsic, stat);
368     }
369     DoTotalReduction<Type>(x, dim, mask, accumulator, intrinsic, terminator);
370     accumulator.GetResult(result.OffsetElement<Type>());
371   } else {
372     // Partial reduction
373     PartialReduction<ACCUMULATOR, CAT, KIND>(
374         result, x, dim, mask, terminator, intrinsic, accumulator);
375   }
376 }
377 
378 template <TypeCategory CAT, bool IS_MAXVAL> struct MaxOrMinHelper {
379   template <int KIND> struct Functor {
380     void operator()(Descriptor &result, const Descriptor &x, int dim,
381         const Descriptor *mask, const char *intrinsic,
382         Terminator &terminator) const {
383       DoMaxMinNorm2<CAT, KIND,
384           NumericExtremumAccumulator<CAT, KIND, IS_MAXVAL>>(
385           result, x, dim, mask, intrinsic, terminator);
386     }
387   };
388 };
389 
390 template <bool IS_MAXVAL>
391 inline void NumericMaxOrMin(Descriptor &result, const Descriptor &x, int dim,
392     const char *source, int line, const Descriptor *mask,
393     const char *intrinsic) {
394   Terminator terminator{source, line};
395   auto type{x.type().GetCategoryAndKind()};
396   RUNTIME_CHECK(terminator, type);
397   switch (type->first) {
398   case TypeCategory::Integer:
399     ApplyIntegerKind<
400         MaxOrMinHelper<TypeCategory::Integer, IS_MAXVAL>::template Functor,
401         void>(
402         type->second, terminator, result, x, dim, mask, intrinsic, terminator);
403     break;
404   case TypeCategory::Real:
405     ApplyFloatingPointKind<
406         MaxOrMinHelper<TypeCategory::Real, IS_MAXVAL>::template Functor, void>(
407         type->second, terminator, result, x, dim, mask, intrinsic, terminator);
408     break;
409   default:
410     terminator.Crash("%s: bad type code %d", intrinsic, x.type().raw());
411   }
412 }
413 
414 template <int KIND, bool IS_MAXVAL> class CharacterExtremumAccumulator {
415 public:
416   using Type = CppTypeFor<TypeCategory::Character, KIND>;
417   explicit CharacterExtremumAccumulator(const Descriptor &array)
418       : array_{array}, charLen_{array_.ElementBytes() / KIND} {}
419   void Reinitialize() { extremum_ = nullptr; }
420   template <typename A> void GetResult(A *p, int /*zeroBasedDim*/ = -1) const {
421     static_assert(std::is_same_v<A, Type>);
422     if (extremum_) {
423       std::memcpy(p, extremum_, charLen_);
424     } else {
425       // empty array: result is all zero-valued characters
426       std::memset(p, 0, charLen_);
427     }
428   }
429   bool Accumulate(const Type *x) {
430     if (!extremum_) {
431       extremum_ = x;
432     } else {
433       int cmp{CharacterScalarCompare(x, extremum_, charLen_, charLen_)};
434       if (IS_MAXVAL == (cmp > 0)) {
435         extremum_ = x;
436       }
437     }
438     return true;
439   }
440   template <typename A> bool AccumulateAt(const SubscriptValue at[]) {
441     return Accumulate(array_.Element<A>(at));
442   }
443 
444 private:
445   const Descriptor &array_;
446   std::size_t charLen_;
447   const Type *extremum_{nullptr};
448 };
449 
450 template <bool IS_MAXVAL> struct CharacterMaxOrMinHelper {
451   template <int KIND> struct Functor {
452     void operator()(Descriptor &result, const Descriptor &x, int dim,
453         const Descriptor *mask, const char *intrinsic,
454         Terminator &terminator) const {
455       DoMaxMinNorm2<TypeCategory::Character, KIND,
456           CharacterExtremumAccumulator<KIND, IS_MAXVAL>>(
457           result, x, dim, mask, intrinsic, terminator);
458     }
459   };
460 };
461 
462 template <bool IS_MAXVAL>
463 inline void CharacterMaxOrMin(Descriptor &result, const Descriptor &x, int dim,
464     const char *source, int line, const Descriptor *mask,
465     const char *intrinsic) {
466   Terminator terminator{source, line};
467   auto type{x.type().GetCategoryAndKind()};
468   RUNTIME_CHECK(terminator, type && type->first == TypeCategory::Character);
469   ApplyCharacterKind<CharacterMaxOrMinHelper<IS_MAXVAL>::template Functor,
470       void>(
471       type->second, terminator, result, x, dim, mask, intrinsic, terminator);
472 }
473 
474 extern "C" {
475 CppTypeFor<TypeCategory::Integer, 1> RTNAME(MaxvalInteger1)(const Descriptor &x,
476     const char *source, int line, int dim, const Descriptor *mask) {
477   return TotalNumericMaxOrMin<TypeCategory::Integer, 1, true>(
478       x, source, line, dim, mask, "MAXVAL");
479 }
480 CppTypeFor<TypeCategory::Integer, 2> RTNAME(MaxvalInteger2)(const Descriptor &x,
481     const char *source, int line, int dim, const Descriptor *mask) {
482   return TotalNumericMaxOrMin<TypeCategory::Integer, 2, true>(
483       x, source, line, dim, mask, "MAXVAL");
484 }
485 CppTypeFor<TypeCategory::Integer, 4> RTNAME(MaxvalInteger4)(const Descriptor &x,
486     const char *source, int line, int dim, const Descriptor *mask) {
487   return TotalNumericMaxOrMin<TypeCategory::Integer, 4, true>(
488       x, source, line, dim, mask, "MAXVAL");
489 }
490 CppTypeFor<TypeCategory::Integer, 8> RTNAME(MaxvalInteger8)(const Descriptor &x,
491     const char *source, int line, int dim, const Descriptor *mask) {
492   return TotalNumericMaxOrMin<TypeCategory::Integer, 8, true>(
493       x, source, line, dim, mask, "MAXVAL");
494 }
495 #ifdef __SIZEOF_INT128__
496 CppTypeFor<TypeCategory::Integer, 16> RTNAME(MaxvalInteger16)(
497     const Descriptor &x, const char *source, int line, int dim,
498     const Descriptor *mask) {
499   return TotalNumericMaxOrMin<TypeCategory::Integer, 16, true>(
500       x, source, line, dim, mask, "MAXVAL");
501 }
502 #endif
503 
504 // TODO: REAL(2 & 3)
505 CppTypeFor<TypeCategory::Real, 4> RTNAME(MaxvalReal4)(const Descriptor &x,
506     const char *source, int line, int dim, const Descriptor *mask) {
507   return TotalNumericMaxOrMin<TypeCategory::Real, 4, true>(
508       x, source, line, dim, mask, "MAXVAL");
509 }
510 CppTypeFor<TypeCategory::Real, 8> RTNAME(MaxvalReal8)(const Descriptor &x,
511     const char *source, int line, int dim, const Descriptor *mask) {
512   return TotalNumericMaxOrMin<TypeCategory::Real, 8, true>(
513       x, source, line, dim, mask, "MAXVAL");
514 }
515 #if LDBL_MANT_DIG == 64
516 CppTypeFor<TypeCategory::Real, 10> RTNAME(MaxvalReal10)(const Descriptor &x,
517     const char *source, int line, int dim, const Descriptor *mask) {
518   return TotalNumericMaxOrMin<TypeCategory::Real, 10, true>(
519       x, source, line, dim, mask, "MAXVAL");
520 }
521 #endif
522 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
523 CppTypeFor<TypeCategory::Real, 16> RTNAME(MaxvalReal16)(const Descriptor &x,
524     const char *source, int line, int dim, const Descriptor *mask) {
525   return TotalNumericMaxOrMin<TypeCategory::Real, 16, true>(
526       x, source, line, dim, mask, "MAXVAL");
527 }
528 #endif
529 
530 void RTNAME(MaxvalCharacter)(Descriptor &result, const Descriptor &x,
531     const char *source, int line, const Descriptor *mask) {
532   CharacterMaxOrMin<true>(result, x, 0, source, line, mask, "MAXVAL");
533 }
534 
535 CppTypeFor<TypeCategory::Integer, 1> RTNAME(MinvalInteger1)(const Descriptor &x,
536     const char *source, int line, int dim, const Descriptor *mask) {
537   return TotalNumericMaxOrMin<TypeCategory::Integer, 1, false>(
538       x, source, line, dim, mask, "MINVAL");
539 }
540 CppTypeFor<TypeCategory::Integer, 2> RTNAME(MinvalInteger2)(const Descriptor &x,
541     const char *source, int line, int dim, const Descriptor *mask) {
542   return TotalNumericMaxOrMin<TypeCategory::Integer, 2, false>(
543       x, source, line, dim, mask, "MINVAL");
544 }
545 CppTypeFor<TypeCategory::Integer, 4> RTNAME(MinvalInteger4)(const Descriptor &x,
546     const char *source, int line, int dim, const Descriptor *mask) {
547   return TotalNumericMaxOrMin<TypeCategory::Integer, 4, false>(
548       x, source, line, dim, mask, "MINVAL");
549 }
550 CppTypeFor<TypeCategory::Integer, 8> RTNAME(MinvalInteger8)(const Descriptor &x,
551     const char *source, int line, int dim, const Descriptor *mask) {
552   return TotalNumericMaxOrMin<TypeCategory::Integer, 8, false>(
553       x, source, line, dim, mask, "MINVAL");
554 }
555 #ifdef __SIZEOF_INT128__
556 CppTypeFor<TypeCategory::Integer, 16> RTNAME(MinvalInteger16)(
557     const Descriptor &x, const char *source, int line, int dim,
558     const Descriptor *mask) {
559   return TotalNumericMaxOrMin<TypeCategory::Integer, 16, false>(
560       x, source, line, dim, mask, "MINVAL");
561 }
562 #endif
563 
564 // TODO: REAL(2 & 3)
565 CppTypeFor<TypeCategory::Real, 4> RTNAME(MinvalReal4)(const Descriptor &x,
566     const char *source, int line, int dim, const Descriptor *mask) {
567   return TotalNumericMaxOrMin<TypeCategory::Real, 4, false>(
568       x, source, line, dim, mask, "MINVAL");
569 }
570 CppTypeFor<TypeCategory::Real, 8> RTNAME(MinvalReal8)(const Descriptor &x,
571     const char *source, int line, int dim, const Descriptor *mask) {
572   return TotalNumericMaxOrMin<TypeCategory::Real, 8, false>(
573       x, source, line, dim, mask, "MINVAL");
574 }
575 #if LDBL_MANT_DIG == 64
576 CppTypeFor<TypeCategory::Real, 10> RTNAME(MinvalReal10)(const Descriptor &x,
577     const char *source, int line, int dim, const Descriptor *mask) {
578   return TotalNumericMaxOrMin<TypeCategory::Real, 10, false>(
579       x, source, line, dim, mask, "MINVAL");
580 }
581 #endif
582 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
583 CppTypeFor<TypeCategory::Real, 16> RTNAME(MinvalReal16)(const Descriptor &x,
584     const char *source, int line, int dim, const Descriptor *mask) {
585   return TotalNumericMaxOrMin<TypeCategory::Real, 16, false>(
586       x, source, line, dim, mask, "MINVAL");
587 }
588 #endif
589 
590 void RTNAME(MinvalCharacter)(Descriptor &result, const Descriptor &x,
591     const char *source, int line, const Descriptor *mask) {
592   CharacterMaxOrMin<false>(result, x, 0, source, line, mask, "MINVAL");
593 }
594 
595 void RTNAME(MaxvalDim)(Descriptor &result, const Descriptor &x, int dim,
596     const char *source, int line, const Descriptor *mask) {
597   if (x.type().IsCharacter()) {
598     CharacterMaxOrMin<true>(result, x, dim, source, line, mask, "MAXVAL");
599   } else {
600     NumericMaxOrMin<true>(result, x, dim, source, line, mask, "MAXVAL");
601   }
602 }
603 void RTNAME(MinvalDim)(Descriptor &result, const Descriptor &x, int dim,
604     const char *source, int line, const Descriptor *mask) {
605   if (x.type().IsCharacter()) {
606     CharacterMaxOrMin<false>(result, x, dim, source, line, mask, "MINVAL");
607   } else {
608     NumericMaxOrMin<false>(result, x, dim, source, line, mask, "MINVAL");
609   }
610 }
611 } // extern "C"
612 
613 // NORM2
614 
615 template <int KIND> class Norm2Accumulator {
616 public:
617   using Type = CppTypeFor<TypeCategory::Real, KIND>;
618   // Use at least double precision for accumulators.
619   // Don't use __float128, it doesn't work with abs() or sqrt() yet.
620   static constexpr int largestLDKind {
621 #if LDBL_MANT_DIG == 113
622     16
623 #elif LDBL_MANT_DIG == 64
624     10
625 #else
626     8
627 #endif
628   };
629   using AccumType = CppTypeFor<TypeCategory::Real,
630       std::max(std::min(largestLDKind, KIND), 8)>;
631   explicit Norm2Accumulator(const Descriptor &array) : array_{array} {}
632   void Reinitialize() { max_ = sum_ = 0; }
633   template <typename A> void GetResult(A *p, int /*zeroBasedDim*/ = -1) const {
634     // m * sqrt(1 + sum((others(:)/m)**2))
635     *p = static_cast<Type>(max_ * std::sqrt(1 + sum_));
636   }
637   bool Accumulate(Type x) {
638     auto absX{std::abs(static_cast<AccumType>(x))};
639     if (!max_) {
640       max_ = x;
641     } else if (absX > max_) {
642       auto t{max_ / absX}; // < 1.0
643       auto tsq{t * t};
644       sum_ *= tsq; // scale sum to reflect change to the max
645       sum_ += tsq; // include a term for the previous max
646       max_ = absX;
647     } else { // absX <= max_
648       auto t{absX / max_};
649       sum_ += t * t;
650     }
651     return true;
652   }
653   template <typename A> bool AccumulateAt(const SubscriptValue at[]) {
654     return Accumulate(*array_.Element<A>(at));
655   }
656 
657 private:
658   const Descriptor &array_;
659   AccumType max_{0}; // value (m) with largest magnitude
660   AccumType sum_{0}; // sum((others(:)/m)**2)
661 };
662 
663 template <int KIND> struct Norm2Helper {
664   void operator()(Descriptor &result, const Descriptor &x, int dim,
665       const Descriptor *mask, Terminator &terminator) const {
666     DoMaxMinNorm2<TypeCategory::Real, KIND, Norm2Accumulator<KIND>>(
667         result, x, dim, mask, "NORM2", terminator);
668   }
669 };
670 
671 extern "C" {
672 // TODO: REAL(2 & 3)
673 CppTypeFor<TypeCategory::Real, 4> RTNAME(Norm2_4)(const Descriptor &x,
674     const char *source, int line, int dim, const Descriptor *mask) {
675   return GetTotalReduction<TypeCategory::Real, 4>(
676       x, source, line, dim, mask, Norm2Accumulator<4>{x}, "NORM2");
677 }
678 CppTypeFor<TypeCategory::Real, 8> RTNAME(Norm2_8)(const Descriptor &x,
679     const char *source, int line, int dim, const Descriptor *mask) {
680   return GetTotalReduction<TypeCategory::Real, 8>(
681       x, source, line, dim, mask, Norm2Accumulator<8>{x}, "NORM2");
682 }
683 #if LDBL_MANT_DIG == 64
684 CppTypeFor<TypeCategory::Real, 10> RTNAME(Norm2_10)(const Descriptor &x,
685     const char *source, int line, int dim, const Descriptor *mask) {
686   return GetTotalReduction<TypeCategory::Real, 10>(
687       x, source, line, dim, mask, Norm2Accumulator<10>{x}, "NORM2");
688 }
689 #endif
690 #if LDBL_MANT_DIG == 113
691 CppTypeFor<TypeCategory::Real, 16> RTNAME(Norm2_16)(const Descriptor &x,
692     const char *source, int line, int dim, const Descriptor *mask) {
693   return GetTotalReduction<TypeCategory::Real, 16>(
694       x, source, line, dim, mask, Norm2Accumulator<16>{x}, "NORM2");
695 }
696 #endif
697 
698 void RTNAME(Norm2Dim)(Descriptor &result, const Descriptor &x, int dim,
699     const char *source, int line, const Descriptor *mask) {
700   Terminator terminator{source, line};
701   auto type{x.type().GetCategoryAndKind()};
702   RUNTIME_CHECK(terminator, type);
703   if (type->first == TypeCategory::Real) {
704     ApplyFloatingPointKind<Norm2Helper, void>(
705         type->second, terminator, result, x, dim, mask, terminator);
706   } else {
707     terminator.Crash("NORM2: bad type code %d", x.type().raw());
708   }
709 }
710 } // extern "C"
711 } // namespace Fortran::runtime
712