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

source code of flang/runtime/extrema.cpp