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 | |
24 | namespace Fortran::runtime { |
25 | |
26 | // MAXLOC & MINLOC |
27 | |
28 | template <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 | |
44 | template <typename T, bool IS_MAX, bool BACK> class CharacterCompare { |
45 | public: |
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 | |
60 | private: |
61 | std::size_t chars_; |
62 | }; |
63 | |
64 | template <typename COMPARE> class ExtremumLocAccumulator { |
65 | public: |
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 | |
101 | private: |
102 | const Descriptor &array_; |
103 | int argRank_; |
104 | SubscriptValue extremumLoc_[maxRank]; |
105 | const Type *previous_{nullptr}; |
106 | COMPARE compare_; |
107 | }; |
108 | |
109 | template <typename ACCUMULATOR, typename CPPTYPE> |
110 | static 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 | |
119 | template <TypeCategory CAT, int KIND, bool IS_MAX, |
120 | template <typename, bool, bool> class COMPARE> |
121 | inline 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 | |
135 | template <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 | |
146 | template <bool IS_MAX> |
147 | inline 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 | |
175 | template <TypeCategory CAT, int KIND, bool IS_MAXVAL> |
176 | inline 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 | |
195 | extern "C" { |
196 | RT_EXT_API_GROUP_BEGIN |
197 | |
198 | void 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 | } |
203 | void 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 | } |
208 | void 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 | } |
213 | void 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 | } |
218 | void 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__ |
224 | void 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 |
230 | void 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 | } |
235 | void 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 |
241 | void 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 |
248 | void 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 |
254 | void 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 | } |
259 | void 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 | } |
264 | void 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 | } |
269 | void 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 | } |
274 | void 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__ |
280 | void 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 |
286 | void 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 | } |
291 | void 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 |
297 | void 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 |
304 | void 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 | |
311 | RT_EXT_API_GROUP_END |
312 | } // extern "C" |
313 | |
314 | // MAXLOC/MINLOC with DIM= |
315 | |
316 | template <TypeCategory CAT, int KIND, bool IS_MAX, |
317 | template <typename, bool, bool> class COMPARE, bool BACK> |
318 | static 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 | |
329 | template <TypeCategory CAT, int KIND, bool IS_MAX, |
330 | template <typename, bool, bool> class COMPARE> |
331 | inline 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 | |
343 | template <TypeCategory CAT, bool IS_MAX, |
344 | template <typename, bool, bool> class COMPARE> |
345 | struct 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 | |
356 | template <bool IS_MAX> |
357 | inline 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 | |
408 | extern "C" { |
409 | RT_EXT_API_GROUP_BEGIN |
410 | |
411 | void 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 | } |
416 | void 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 | |
422 | RT_EXT_API_GROUP_END |
423 | } // extern "C" |
424 | |
425 | // MAXVAL and MINVAL |
426 | |
427 | template <TypeCategory CAT, int KIND, bool IS_MAXVAL, typename Enable = void> |
428 | struct 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 |
437 | template <bool IS_MAXVAL> |
438 | struct 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. |
454 | template <bool IS_MAXVAL> |
455 | struct 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 | |
483 | template <TypeCategory CAT, int KIND, bool IS_MAXVAL> |
484 | class NumericExtremumAccumulator { |
485 | public: |
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 | |
517 | private: |
518 | const Descriptor &array_; |
519 | bool any_{false}; |
520 | Type extremum_{MaxOrMinIdentity<CAT, KIND, IS_MAXVAL>::Value()}; |
521 | }; |
522 | |
523 | template <TypeCategory CAT, int KIND, bool IS_MAXVAL> |
524 | inline 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 | |
531 | template <TypeCategory CAT, int KIND, typename ACCUMULATOR> |
532 | static 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 | |
560 | template <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 | |
572 | template <bool IS_MAXVAL> |
573 | inline 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 | |
596 | template <int KIND, bool IS_MAXVAL> class CharacterExtremumAccumulator { |
597 | public: |
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 | |
630 | private: |
631 | const Descriptor &array_; |
632 | std::size_t charLen_; |
633 | const Type *extremum_{nullptr}; |
634 | }; |
635 | |
636 | template <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 | |
648 | template <bool IS_MAXVAL> |
649 | inline 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 | |
660 | extern "C" { |
661 | RT_EXT_API_GROUP_BEGIN |
662 | |
663 | CppTypeFor<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 | } |
668 | CppTypeFor<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 | } |
673 | CppTypeFor<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 | } |
678 | CppTypeFor<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__ |
684 | CppTypeFor<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) |
693 | CppTypeFor<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 | } |
698 | CppTypeFor<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 |
704 | CppTypeFor<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 |
711 | CppTypeFor<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 | |
718 | void 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 | |
723 | CppTypeFor<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 | } |
728 | CppTypeFor<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 | } |
733 | CppTypeFor<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 | } |
738 | CppTypeFor<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__ |
744 | CppTypeFor<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) |
753 | CppTypeFor<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 | } |
758 | CppTypeFor<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 |
764 | CppTypeFor<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 |
771 | CppTypeFor<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 | |
778 | void 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 | |
783 | void 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 | } |
791 | void 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 | |
800 | RT_EXT_API_GROUP_END |
801 | } // extern "C" |
802 | |
803 | // NORM2 |
804 | |
805 | RT_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. |
809 | static 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 | |
819 | RT_VAR_GROUP_END |
820 | |
821 | template <int KIND> class Norm2Accumulator { |
822 | public: |
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 | |
855 | private: |
856 | const Descriptor &array_; |
857 | AccumType max_{0}; // value (m) with largest magnitude |
858 | AccumType sum_{0}; // sum((others(:)/m)**2) |
859 | }; |
860 | |
861 | template <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 | |
869 | extern "C" { |
870 | RT_EXT_API_GROUP_BEGIN |
871 | |
872 | // TODO: REAL(2 & 3) |
873 | CppTypeFor<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 | } |
878 | CppTypeFor<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 |
884 | CppTypeFor<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 |
891 | CppTypeFor<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 | |
898 | void 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 | |
911 | RT_EXT_API_GROUP_END |
912 | } // extern "C" |
913 | } // namespace Fortran::runtime |
914 | |