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