1//===-- lib/runtime/transformational.cpp ------------------------*- C++ -*-===//
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 the transformational intrinsic functions of Fortran 2018 that
10// rearrange or duplicate data without (much) regard to type. These are
11// CSHIFT, EOSHIFT, PACK, RESHAPE, SPREAD, TRANSPOSE, and UNPACK.
12//
13// Many of these are defined in the 2018 standard with text that makes sense
14// only if argument arrays have lower bounds of one. Rather than interpret
15// these cases as implying a hidden constraint, these implementations
16// work with arbitrary lower bounds. This may be technically an extension
17// of the standard but it more likely to conform with its intent.
18
19#include "flang/Runtime/transformational.h"
20#include "copy.h"
21#include "flang-rt/runtime/descriptor.h"
22#include "flang-rt/runtime/terminator.h"
23#include "flang-rt/runtime/tools.h"
24#include "flang-rt/runtime/type-info.h"
25#include "flang/Common/float128.h"
26
27namespace Fortran::runtime {
28
29// Utility for CSHIFT & EOSHIFT rank > 1 cases that determines the shift count
30// for each of the vector sections of the result.
31class ShiftControl {
32public:
33 RT_API_ATTRS ShiftControl(const Descriptor &s, Terminator &t, int dim)
34 : shift_{s}, terminator_{t}, shiftRank_{s.rank()}, dim_{dim} {}
35 RT_API_ATTRS void Init(const Descriptor &source, const char *which) {
36 int rank{source.rank()};
37 RUNTIME_CHECK(terminator_, shiftRank_ == 0 || shiftRank_ == rank - 1);
38 auto catAndKind{shift_.type().GetCategoryAndKind()};
39 RUNTIME_CHECK(
40 terminator_, catAndKind && catAndKind->first == TypeCategory::Integer);
41 shiftElemLen_ = catAndKind->second;
42 if (shiftRank_ > 0) {
43 int k{0};
44 for (int j{0}; j < rank; ++j) {
45 if (j + 1 != dim_) {
46 const Dimension &shiftDim{shift_.GetDimension(k)};
47 lb_[k++] = shiftDim.LowerBound();
48 if (shiftDim.Extent() != source.GetDimension(j).Extent()) {
49 terminator_.Crash("%s: on dimension %d, SHIFT= has extent %jd but "
50 "ARRAY= has extent %jd",
51 which, k, static_cast<std::intmax_t>(shiftDim.Extent()),
52 static_cast<std::intmax_t>(source.GetDimension(j).Extent()));
53 }
54 }
55 }
56 } else if (auto count{GetInt64Safe(
57 shift_.OffsetElement<char>(), shiftElemLen_, terminator_)}) {
58 shiftCount_ = *count;
59 } else {
60 terminator_.Crash("%s: SHIFT= value exceeds 64 bits", which);
61 }
62 }
63 RT_API_ATTRS SubscriptValue GetShift(const SubscriptValue resultAt[]) const {
64 if (shiftRank_ > 0) {
65 SubscriptValue shiftAt[maxRank];
66 int k{0};
67 for (int j{0}; j < shiftRank_ + 1; ++j) {
68 if (j + 1 != dim_) {
69 shiftAt[k] = lb_[k] + resultAt[j] - 1;
70 ++k;
71 }
72 }
73 auto count{GetInt64Safe(
74 shift_.Element<char>(shiftAt), shiftElemLen_, terminator_)};
75 RUNTIME_CHECK(terminator_, count.has_value());
76 return *count;
77 } else {
78 return shiftCount_; // invariant count extracted in Init()
79 }
80 }
81
82private:
83 const Descriptor &shift_;
84 Terminator &terminator_;
85 int shiftRank_;
86 int dim_;
87 SubscriptValue lb_[maxRank];
88 std::size_t shiftElemLen_;
89 SubscriptValue shiftCount_{};
90};
91
92// Fill an EOSHIFT result with default boundary values
93static RT_API_ATTRS void DefaultInitialize(
94 const Descriptor &result, Terminator &terminator) {
95 auto catAndKind{result.type().GetCategoryAndKind()};
96 RUNTIME_CHECK(
97 terminator, catAndKind && catAndKind->first != TypeCategory::Derived);
98 std::size_t elementLen{result.ElementBytes()};
99 std::size_t bytes{result.Elements() * elementLen};
100 if (catAndKind->first == TypeCategory::Character) {
101 switch (int kind{catAndKind->second}) {
102 case 1:
103 Fortran::runtime::fill_n(result.OffsetElement<char>(), bytes, ' ');
104 break;
105 case 2:
106 Fortran::runtime::fill_n(result.OffsetElement<char16_t>(), bytes / 2,
107 static_cast<char16_t>(' '));
108 break;
109 case 4:
110 Fortran::runtime::fill_n(result.OffsetElement<char32_t>(), bytes / 4,
111 static_cast<char32_t>(' '));
112 break;
113 default:
114 terminator.Crash(
115 "not yet implemented: CHARACTER(KIND=%d) in EOSHIFT intrinsic", kind);
116 }
117 } else {
118 std::memset(result.raw().base_addr, 0, bytes);
119 }
120}
121
122static inline RT_API_ATTRS std::size_t AllocateResult(Descriptor &result,
123 const Descriptor &source, int rank, const SubscriptValue extent[],
124 Terminator &terminator, const char *function) {
125 std::size_t elementLen{source.ElementBytes()};
126 const DescriptorAddendum *sourceAddendum{source.Addendum()};
127 result.Establish(source.type(), elementLen, nullptr, rank, extent,
128 CFI_attribute_allocatable, sourceAddendum != nullptr);
129 if (sourceAddendum) {
130 *result.Addendum() = *sourceAddendum;
131 }
132 for (int j{0}; j < rank; ++j) {
133 result.GetDimension(j).SetBounds(1, extent[j]);
134 }
135 if (int stat{result.Allocate(kNoAsyncObject)}) {
136 terminator.Crash(
137 "%s: Could not allocate memory for result (stat=%d)", function, stat);
138 }
139 return elementLen;
140}
141
142template <TypeCategory CAT, int KIND>
143static inline RT_API_ATTRS std::size_t AllocateBesselResult(Descriptor &result,
144 int32_t n1, int32_t n2, Terminator &terminator, const char *function) {
145 int rank{1};
146 SubscriptValue extent[maxRank];
147 for (int j{0}; j < maxRank; j++) {
148 extent[j] = 0;
149 }
150 if (n1 <= n2) {
151 extent[0] = n2 - n1 + 1;
152 }
153
154 std::size_t elementLen{Descriptor::BytesFor(CAT, KIND)};
155 result.Establish(TypeCode{CAT, KIND}, elementLen, nullptr, rank, extent,
156 CFI_attribute_allocatable, false);
157 for (int j{0}; j < rank; ++j) {
158 result.GetDimension(j).SetBounds(1, extent[j]);
159 }
160 if (int stat{result.Allocate(kNoAsyncObject)}) {
161 terminator.Crash(
162 "%s: Could not allocate memory for result (stat=%d)", function, stat);
163 }
164 return elementLen;
165}
166
167template <TypeCategory CAT, int KIND>
168static inline RT_API_ATTRS void DoBesselJn(Descriptor &result, int32_t n1,
169 int32_t n2, CppTypeFor<CAT, KIND> x, CppTypeFor<CAT, KIND> bn2,
170 CppTypeFor<CAT, KIND> bn2_1, const char *sourceFile, int line) {
171 Terminator terminator{sourceFile, line};
172 AllocateBesselResult<CAT, KIND>(result, n1, n2, terminator, "BESSEL_JN");
173
174 // The standard requires that n1 and n2 be non-negative. However, some other
175 // compilers generate results even when n1 and/or n2 are negative. For now,
176 // we also do not enforce the non-negativity constraint.
177 if (n2 < n1) {
178 return;
179 }
180
181 SubscriptValue at[maxRank];
182 for (int j{0}; j < maxRank; ++j) {
183 at[j] = 0;
184 }
185
186 // if n2 >= n1, there will be at least one element in the result.
187 at[0] = n2 - n1 + 1;
188 *result.Element<CppTypeFor<CAT, KIND>>(at) = bn2;
189
190 if (n2 == n1) {
191 return;
192 }
193
194 at[0] = n2 - n1;
195 *result.Element<CppTypeFor<CAT, KIND>>(at) = bn2_1;
196
197 // Bessel functions of the first kind are stable for a backward recursion
198 // (see https://dlmf.nist.gov/10.74.iv and https://dlmf.nist.gov/10.6.E1).
199 //
200 // J(n-1, x) = (2.0 / x) * n * J(n, x) - J(n+1, x)
201 //
202 // which is equivalent to
203 //
204 // J(n, x) = (2.0 / x) * (n + 1) * J(n+1, x) - J(n+2, x)
205 //
206 CppTypeFor<CAT, KIND> bn_2 = bn2;
207 CppTypeFor<CAT, KIND> bn_1 = bn2_1;
208 CppTypeFor<CAT, KIND> twoOverX = 2.0 / x;
209 for (int n{n2 - 2}; n >= n1; --n) {
210 auto bn = twoOverX * (n + 1) * bn_1 - bn_2;
211
212 at[0] = n - n1 + 1;
213 *result.Element<CppTypeFor<CAT, KIND>>(at) = bn;
214
215 bn_2 = bn_1;
216 bn_1 = bn;
217 }
218}
219
220template <TypeCategory CAT, int KIND>
221static inline RT_API_ATTRS void DoBesselJnX0(Descriptor &result, int32_t n1,
222 int32_t n2, const char *sourceFile, int line) {
223 Terminator terminator{sourceFile, line};
224 AllocateBesselResult<CAT, KIND>(result, n1, n2, terminator, "BESSEL_JN");
225
226 // The standard requires that n1 and n2 be non-negative. However, some other
227 // compilers generate results even when n1 and/or n2 are negative. For now,
228 // we also do not enforce the non-negativity constraint.
229 if (n2 < n1) {
230 return;
231 }
232
233 SubscriptValue at[maxRank];
234 for (int j{0}; j < maxRank; ++j) {
235 at[j] = 0;
236 }
237
238 // J(0, 0.0) = 1.0, when n == 0.
239 // J(n, 0.0) = 0.0, when n > 0.
240 at[0] = 1;
241 *result.Element<CppTypeFor<CAT, KIND>>(at) = (n1 == 0) ? 1.0 : 0.0;
242 for (int j{2}; j <= n2 - n1 + 1; ++j) {
243 at[0] = j;
244 *result.Element<CppTypeFor<CAT, KIND>>(at) = 0.0;
245 }
246}
247
248template <TypeCategory CAT, int KIND>
249static inline RT_API_ATTRS void DoBesselYn(Descriptor &result, int32_t n1,
250 int32_t n2, CppTypeFor<CAT, KIND> x, CppTypeFor<CAT, KIND> bn1,
251 CppTypeFor<CAT, KIND> bn1_1, const char *sourceFile, int line) {
252 Terminator terminator{sourceFile, line};
253 AllocateBesselResult<CAT, KIND>(result, n1, n2, terminator, "BESSEL_YN");
254
255 // The standard requires that n1 and n2 be non-negative. However, some other
256 // compilers generate results even when n1 and/or n2 are negative. For now,
257 // we also do not enforce the non-negativity constraint.
258 if (n2 < n1) {
259 return;
260 }
261
262 SubscriptValue at[maxRank];
263 for (int j{0}; j < maxRank; ++j) {
264 at[j] = 0;
265 }
266
267 // if n2 >= n1, there will be at least one element in the result.
268 at[0] = 1;
269 *result.Element<CppTypeFor<CAT, KIND>>(at) = bn1;
270
271 if (n2 == n1) {
272 return;
273 }
274
275 at[0] = 2;
276 *result.Element<CppTypeFor<CAT, KIND>>(at) = bn1_1;
277
278 // Bessel functions of the second kind are stable for a forward recursion
279 // (see https://dlmf.nist.gov/10.74.iv and https://dlmf.nist.gov/10.6.E1).
280 //
281 // Y(n+1, x) = (2.0 / x) * n * Y(n, x) - Y(n-1, x)
282 //
283 // which is equivalent to
284 //
285 // Y(n, x) = (2.0 / x) * (n - 1) * Y(n-1, x) - Y(n-2, x)
286 //
287 CppTypeFor<CAT, KIND> bn_2 = bn1;
288 CppTypeFor<CAT, KIND> bn_1 = bn1_1;
289 CppTypeFor<CAT, KIND> twoOverX = 2.0 / x;
290 for (int n{n1 + 2}; n <= n2; ++n) {
291 auto bn = twoOverX * (n - 1) * bn_1 - bn_2;
292
293 at[0] = n - n1 + 1;
294 *result.Element<CppTypeFor<CAT, KIND>>(at) = bn;
295
296 bn_2 = bn_1;
297 bn_1 = bn;
298 }
299}
300
301template <TypeCategory CAT, int KIND>
302static inline RT_API_ATTRS void DoBesselYnX0(Descriptor &result, int32_t n1,
303 int32_t n2, const char *sourceFile, int line) {
304 Terminator terminator{sourceFile, line};
305 AllocateBesselResult<CAT, KIND>(result, n1, n2, terminator, "BESSEL_YN");
306
307 // The standard requires that n1 and n2 be non-negative. However, some other
308 // compilers generate results even when n1 and/or n2 are negative. For now,
309 // we also do not enforce the non-negativity constraint.
310 if (n2 < n1) {
311 return;
312 }
313
314 SubscriptValue at[maxRank];
315 for (int j{0}; j < maxRank; ++j) {
316 at[j] = 0;
317 }
318
319 // Y(n, 0.0) = -Inf, when n >= 0
320 for (int j{1}; j <= n2 - n1 + 1; ++j) {
321 at[0] = j;
322 *result.Element<CppTypeFor<CAT, KIND>>(at) =
323 -std::numeric_limits<CppTypeFor<CAT, KIND>>::infinity();
324 }
325}
326
327static inline RT_API_ATTRS void CheckConformabilityForShallowCopy(
328 const Descriptor &d1, const Descriptor &d2, Terminator &terminator,
329 const char *funcName, const char *d1Name, const char *d2Name) {
330 if (d1.rank() != d2.rank()) {
331 terminator.Crash(
332 "Incompatible arguments to %s: %s has rank %d, %s has rank %d",
333 funcName, d1Name, d1.rank(), d1Name, d2.rank());
334 }
335
336 // Check that the shapes conform.
337 CheckConformability(d1, d2, terminator, funcName, d1Name, d2Name);
338
339 if (d1.ElementBytes() != d2.ElementBytes()) {
340 terminator.Crash("Incompatible arguments to %s: %s has element byte length "
341 "%zd, %s has length %zd",
342 funcName, d1Name, d1.ElementBytes(), d2Name, d2.ElementBytes());
343 }
344 if (d1.type() != d2.type()) {
345 terminator.Crash("Incompatible arguments to %s: %s has type code %d, %s "
346 "has type code %d",
347 funcName, d1Name, d1.type().raw(), d2Name, d2.type().raw());
348 }
349 const DescriptorAddendum *d1Addendum{d1.Addendum()};
350 const typeInfo::DerivedType *d1Derived{
351 d1Addendum ? d1Addendum->derivedType() : nullptr};
352 const DescriptorAddendum *d2Addendum{d2.Addendum()};
353 const typeInfo::DerivedType *d2Derived{
354 d2Addendum ? d2Addendum->derivedType() : nullptr};
355 if (d1Derived != d2Derived) {
356 terminator.Crash(
357 "Incompatible arguments to %s: %s and %s have different derived types",
358 funcName, d1Name, d2Name);
359 }
360 if (d2Derived) {
361 // Compare LEN parameters.
362 std::size_t lenParms{d2Derived->LenParameters()};
363 for (std::size_t j{0}; j < lenParms; ++j) {
364 if (d1Addendum->LenParameterValue(j) !=
365 d2Addendum->LenParameterValue(j)) {
366 terminator.Crash("Incompatible arguments to %s: type length parameter "
367 "%zd for %s is %zd, for %s is %zd",
368 funcName, j, d1Name,
369 static_cast<std::size_t>(d1Addendum->LenParameterValue(j)), d2Name,
370 static_cast<std::size_t>(d2Addendum->LenParameterValue(j)));
371 }
372 }
373 }
374}
375
376template <bool IS_ALLOCATING>
377static inline RT_API_ATTRS void DoShallowCopy(
378 std::conditional_t<IS_ALLOCATING, Descriptor, const Descriptor> &result,
379 const Descriptor &source, Terminator &terminator, const char *funcName) {
380 if constexpr (IS_ALLOCATING) {
381 SubscriptValue extent[maxRank];
382 source.GetShape(extent);
383 AllocateResult(result, source, source.rank(), extent, terminator, funcName);
384 } else {
385 CheckConformabilityForShallowCopy(
386 result, source, terminator, funcName, "RESULT=", "SOURCE=");
387 }
388
389 ShallowCopy(result, source);
390}
391
392extern "C" {
393RT_EXT_API_GROUP_BEGIN
394
395// BESSEL_JN
396// TODO: REAL(2 & 3)
397void RTDEF(BesselJn_4)(Descriptor &result, int32_t n1, int32_t n2,
398 CppTypeFor<TypeCategory::Real, 4> x, CppTypeFor<TypeCategory::Real, 4> bn2,
399 CppTypeFor<TypeCategory::Real, 4> bn2_1, const char *sourceFile, int line) {
400 DoBesselJn<TypeCategory::Real, 4>(
401 result, n1, n2, x, bn2, bn2_1, sourceFile, line);
402}
403
404void RTDEF(BesselJn_8)(Descriptor &result, int32_t n1, int32_t n2,
405 CppTypeFor<TypeCategory::Real, 8> x, CppTypeFor<TypeCategory::Real, 8> bn2,
406 CppTypeFor<TypeCategory::Real, 8> bn2_1, const char *sourceFile, int line) {
407 DoBesselJn<TypeCategory::Real, 8>(
408 result, n1, n2, x, bn2, bn2_1, sourceFile, line);
409}
410
411#if HAS_FLOAT80
412void RTDEF(BesselJn_10)(Descriptor &result, int32_t n1, int32_t n2,
413 CppTypeFor<TypeCategory::Real, 10> x,
414 CppTypeFor<TypeCategory::Real, 10> bn2,
415 CppTypeFor<TypeCategory::Real, 10> bn2_1, const char *sourceFile,
416 int line) {
417 DoBesselJn<TypeCategory::Real, 10>(
418 result, n1, n2, x, bn2, bn2_1, sourceFile, line);
419}
420#endif
421
422#if HAS_LDBL128 || HAS_FLOAT128
423void RTDEF(BesselJn_16)(Descriptor &result, int32_t n1, int32_t n2,
424 CppTypeFor<TypeCategory::Real, 16> x,
425 CppTypeFor<TypeCategory::Real, 16> bn2,
426 CppTypeFor<TypeCategory::Real, 16> bn2_1, const char *sourceFile,
427 int line) {
428 DoBesselJn<TypeCategory::Real, 16>(
429 result, n1, n2, x, bn2, bn2_1, sourceFile, line);
430}
431#endif
432
433// TODO: REAL(2 & 3)
434void RTDEF(BesselJnX0_4)(Descriptor &result, int32_t n1, int32_t n2,
435 const char *sourceFile, int line) {
436 DoBesselJnX0<TypeCategory::Real, 4>(result, n1, n2, sourceFile, line);
437}
438
439void RTDEF(BesselJnX0_8)(Descriptor &result, int32_t n1, int32_t n2,
440 const char *sourceFile, int line) {
441 DoBesselJnX0<TypeCategory::Real, 8>(result, n1, n2, sourceFile, line);
442}
443
444#if HAS_FLOAT80
445void RTDEF(BesselJnX0_10)(Descriptor &result, int32_t n1, int32_t n2,
446 const char *sourceFile, int line) {
447 DoBesselJnX0<TypeCategory::Real, 10>(result, n1, n2, sourceFile, line);
448}
449#endif
450
451#if HAS_LDBL128 || HAS_FLOAT128
452void RTDEF(BesselJnX0_16)(Descriptor &result, int32_t n1, int32_t n2,
453 const char *sourceFile, int line) {
454 DoBesselJnX0<TypeCategory::Real, 16>(result, n1, n2, sourceFile, line);
455}
456#endif
457
458// BESSEL_YN
459// TODO: REAL(2 & 3)
460void RTDEF(BesselYn_4)(Descriptor &result, int32_t n1, int32_t n2,
461 CppTypeFor<TypeCategory::Real, 4> x, CppTypeFor<TypeCategory::Real, 4> bn1,
462 CppTypeFor<TypeCategory::Real, 4> bn1_1, const char *sourceFile, int line) {
463 DoBesselYn<TypeCategory::Real, 4>(
464 result, n1, n2, x, bn1, bn1_1, sourceFile, line);
465}
466
467void RTDEF(BesselYn_8)(Descriptor &result, int32_t n1, int32_t n2,
468 CppTypeFor<TypeCategory::Real, 8> x, CppTypeFor<TypeCategory::Real, 8> bn1,
469 CppTypeFor<TypeCategory::Real, 8> bn1_1, const char *sourceFile, int line) {
470 DoBesselYn<TypeCategory::Real, 8>(
471 result, n1, n2, x, bn1, bn1_1, sourceFile, line);
472}
473
474#if HAS_FLOAT80
475void RTDEF(BesselYn_10)(Descriptor &result, int32_t n1, int32_t n2,
476 CppTypeFor<TypeCategory::Real, 10> x,
477 CppTypeFor<TypeCategory::Real, 10> bn1,
478 CppTypeFor<TypeCategory::Real, 10> bn1_1, const char *sourceFile,
479 int line) {
480 DoBesselYn<TypeCategory::Real, 10>(
481 result, n1, n2, x, bn1, bn1_1, sourceFile, line);
482}
483#endif
484
485#if HAS_LDBL128 || HAS_FLOAT128
486void RTDEF(BesselYn_16)(Descriptor &result, int32_t n1, int32_t n2,
487 CppTypeFor<TypeCategory::Real, 16> x,
488 CppTypeFor<TypeCategory::Real, 16> bn1,
489 CppTypeFor<TypeCategory::Real, 16> bn1_1, const char *sourceFile,
490 int line) {
491 DoBesselYn<TypeCategory::Real, 16>(
492 result, n1, n2, x, bn1, bn1_1, sourceFile, line);
493}
494#endif
495
496// TODO: REAL(2 & 3)
497void RTDEF(BesselYnX0_4)(Descriptor &result, int32_t n1, int32_t n2,
498 const char *sourceFile, int line) {
499 DoBesselYnX0<TypeCategory::Real, 4>(result, n1, n2, sourceFile, line);
500}
501
502void RTDEF(BesselYnX0_8)(Descriptor &result, int32_t n1, int32_t n2,
503 const char *sourceFile, int line) {
504 DoBesselYnX0<TypeCategory::Real, 8>(result, n1, n2, sourceFile, line);
505}
506
507#if HAS_FLOAT80
508void RTDEF(BesselYnX0_10)(Descriptor &result, int32_t n1, int32_t n2,
509 const char *sourceFile, int line) {
510 DoBesselYnX0<TypeCategory::Real, 10>(result, n1, n2, sourceFile, line);
511}
512#endif
513
514#if HAS_LDBL128 || HAS_FLOAT128
515void RTDEF(BesselYnX0_16)(Descriptor &result, int32_t n1, int32_t n2,
516 const char *sourceFile, int line) {
517 DoBesselYnX0<TypeCategory::Real, 16>(result, n1, n2, sourceFile, line);
518}
519#endif
520
521// CSHIFT where rank of ARRAY argument > 1
522void RTDEF(Cshift)(Descriptor &result, const Descriptor &source,
523 const Descriptor &shift, int dim, const char *sourceFile, int line) {
524 Terminator terminator{sourceFile, line};
525 int rank{source.rank()};
526 RUNTIME_CHECK(terminator, rank > 1);
527 if (dim < 1 || dim > rank) {
528 terminator.Crash(
529 "CSHIFT: DIM=%d must be >= 1 and <= ARRAY= rank %d", dim, rank);
530 }
531 ShiftControl shiftControl{shift, terminator, dim};
532 shiftControl.Init(source, "CSHIFT");
533 SubscriptValue extent[maxRank];
534 source.GetShape(extent);
535 AllocateResult(result, source, rank, extent, terminator, "CSHIFT");
536 SubscriptValue resultAt[maxRank];
537 for (int j{0}; j < rank; ++j) {
538 resultAt[j] = 1;
539 }
540 SubscriptValue sourceLB[maxRank];
541 source.GetLowerBounds(sourceLB);
542 SubscriptValue dimExtent{extent[dim - 1]};
543 SubscriptValue dimLB{sourceLB[dim - 1]};
544 SubscriptValue &resDim{resultAt[dim - 1]};
545 for (std::size_t n{result.Elements()}; n > 0; n -= dimExtent) {
546 SubscriptValue shiftCount{shiftControl.GetShift(resultAt)};
547 SubscriptValue sourceAt[maxRank];
548 for (int j{0}; j < rank; ++j) {
549 sourceAt[j] = sourceLB[j] + resultAt[j] - 1;
550 }
551 SubscriptValue &sourceDim{sourceAt[dim - 1]};
552 sourceDim = dimLB + shiftCount % dimExtent;
553 if (sourceDim < dimLB) {
554 sourceDim += dimExtent;
555 }
556 for (resDim = 1; resDim <= dimExtent; ++resDim) {
557 CopyElement(result, resultAt, source, sourceAt, terminator);
558 if (++sourceDim == dimLB + dimExtent) {
559 sourceDim = dimLB;
560 }
561 }
562 result.IncrementSubscripts(resultAt);
563 }
564}
565
566// CSHIFT where rank of ARRAY argument == 1
567void RTDEF(CshiftVector)(Descriptor &result, const Descriptor &source,
568 std::int64_t shift, const char *sourceFile, int line) {
569 Terminator terminator{sourceFile, line};
570 RUNTIME_CHECK(terminator, source.rank() == 1);
571 const Dimension &sourceDim{source.GetDimension(0)};
572 SubscriptValue extent{sourceDim.Extent()};
573 AllocateResult(result, source, 1, &extent, terminator, "CSHIFT");
574 SubscriptValue lb{sourceDim.LowerBound()};
575 for (SubscriptValue j{0}; j < extent; ++j) {
576 SubscriptValue resultAt{1 + j};
577 SubscriptValue sourceAt{
578 lb + static_cast<SubscriptValue>(j + shift) % extent};
579 if (sourceAt < lb) {
580 sourceAt += extent;
581 }
582 CopyElement(result, &resultAt, source, &sourceAt, terminator);
583 }
584}
585
586// EOSHIFT of rank > 1
587void RTDEF(Eoshift)(Descriptor &result, const Descriptor &source,
588 const Descriptor &shift, const Descriptor *boundary, int dim,
589 const char *sourceFile, int line) {
590 Terminator terminator{sourceFile, line};
591 SubscriptValue extent[maxRank];
592 int rank{source.GetShape(extent)};
593 RUNTIME_CHECK(terminator, rank > 1);
594 if (dim < 1 || dim > rank) {
595 terminator.Crash(
596 "EOSHIFT: DIM=%d must be >= 1 and <= ARRAY= rank %d", dim, rank);
597 }
598 std::size_t elementLen{
599 AllocateResult(result, source, rank, extent, terminator, "EOSHIFT")};
600 int boundaryRank{-1};
601 if (boundary) {
602 boundaryRank = boundary->rank();
603 RUNTIME_CHECK(terminator, boundaryRank == 0 || boundaryRank == rank - 1);
604 RUNTIME_CHECK(terminator, boundary->type() == source.type());
605 if (boundary->ElementBytes() != elementLen) {
606 terminator.Crash("EOSHIFT: BOUNDARY= has element byte length %zd, but "
607 "ARRAY= has length %zd",
608 boundary->ElementBytes(), elementLen);
609 }
610 if (boundaryRank > 0) {
611 int k{0};
612 for (int j{0}; j < rank; ++j) {
613 if (j != dim - 1) {
614 if (boundary->GetDimension(k).Extent() != extent[j]) {
615 terminator.Crash("EOSHIFT: BOUNDARY= has extent %jd on dimension "
616 "%d but must conform with extent %jd of ARRAY=",
617 static_cast<std::intmax_t>(boundary->GetDimension(k).Extent()),
618 k + 1, static_cast<std::intmax_t>(extent[j]));
619 }
620 ++k;
621 }
622 }
623 }
624 }
625 ShiftControl shiftControl{shift, terminator, dim};
626 shiftControl.Init(source, "EOSHIFT");
627 SubscriptValue resultAt[maxRank];
628 for (int j{0}; j < rank; ++j) {
629 resultAt[j] = 1;
630 }
631 if (!boundary) {
632 DefaultInitialize(result, terminator);
633 }
634 SubscriptValue sourceLB[maxRank];
635 source.GetLowerBounds(sourceLB);
636 SubscriptValue boundaryAt[maxRank];
637 if (boundaryRank > 0) {
638 boundary->GetLowerBounds(boundaryAt);
639 }
640 SubscriptValue dimExtent{extent[dim - 1]};
641 SubscriptValue dimLB{sourceLB[dim - 1]};
642 SubscriptValue &resDim{resultAt[dim - 1]};
643 for (std::size_t n{result.Elements()}; n > 0; n -= dimExtent) {
644 SubscriptValue shiftCount{shiftControl.GetShift(resultAt)};
645 SubscriptValue sourceAt[maxRank];
646 for (int j{0}; j < rank; ++j) {
647 sourceAt[j] = sourceLB[j] + resultAt[j] - 1;
648 }
649 SubscriptValue &sourceDim{sourceAt[dim - 1]};
650 sourceDim = dimLB + shiftCount;
651 for (resDim = 1; resDim <= dimExtent; ++resDim) {
652 if (sourceDim >= dimLB && sourceDim < dimLB + dimExtent) {
653 CopyElement(result, resultAt, source, sourceAt, terminator);
654 } else if (boundary) {
655 CopyElement(result, resultAt, *boundary, boundaryAt, terminator);
656 }
657 ++sourceDim;
658 }
659 result.IncrementSubscripts(resultAt);
660 if (boundaryRank > 0) {
661 boundary->IncrementSubscripts(boundaryAt);
662 }
663 }
664}
665
666// EOSHIFT of vector
667void RTDEF(EoshiftVector)(Descriptor &result, const Descriptor &source,
668 std::int64_t shift, const Descriptor *boundary, const char *sourceFile,
669 int line) {
670 Terminator terminator{sourceFile, line};
671 RUNTIME_CHECK(terminator, source.rank() == 1);
672 SubscriptValue extent{source.GetDimension(0).Extent()};
673 std::size_t elementLen{
674 AllocateResult(result, source, 1, &extent, terminator, "EOSHIFT")};
675 if (boundary) {
676 RUNTIME_CHECK(terminator, boundary->rank() == 0);
677 RUNTIME_CHECK(terminator, boundary->type() == source.type());
678 if (boundary->ElementBytes() != elementLen) {
679 terminator.Crash("EOSHIFT: BOUNDARY= has element byte length %zd but "
680 "ARRAY= has length %zd",
681 boundary->ElementBytes(), elementLen);
682 }
683 }
684 if (!boundary) {
685 DefaultInitialize(result, terminator);
686 }
687 SubscriptValue lb{source.GetDimension(0).LowerBound()};
688 for (SubscriptValue j{1}; j <= extent; ++j) {
689 SubscriptValue sourceAt{lb + j - 1 + static_cast<SubscriptValue>(shift)};
690 if (sourceAt >= lb && sourceAt < lb + extent) {
691 CopyElement(result, &j, source, &sourceAt, terminator);
692 } else if (boundary) {
693 CopyElement(result, &j, *boundary, 0, terminator);
694 }
695 }
696}
697
698// PACK
699void RTDEF(Pack)(Descriptor &result, const Descriptor &source,
700 const Descriptor &mask, const Descriptor *vector, const char *sourceFile,
701 int line) {
702 Terminator terminator{sourceFile, line};
703 CheckConformability(source, mask, terminator, "PACK", "ARRAY=", "MASK=");
704 auto maskType{mask.type().GetCategoryAndKind()};
705 RUNTIME_CHECK(
706 terminator, maskType && maskType->first == TypeCategory::Logical);
707 SubscriptValue trues{0};
708 if (mask.rank() == 0) {
709 if (IsLogicalElementTrue(mask, nullptr)) {
710 trues = source.Elements();
711 }
712 } else {
713 SubscriptValue maskAt[maxRank];
714 mask.GetLowerBounds(maskAt);
715 for (std::size_t n{mask.Elements()}; n > 0; --n) {
716 if (IsLogicalElementTrue(mask, maskAt)) {
717 ++trues;
718 }
719 mask.IncrementSubscripts(maskAt);
720 }
721 }
722 SubscriptValue extent{trues};
723 if (vector) {
724 RUNTIME_CHECK(terminator, vector->rank() == 1);
725 RUNTIME_CHECK(terminator, source.type() == vector->type());
726 if (source.ElementBytes() != vector->ElementBytes()) {
727 terminator.Crash("PACK: ARRAY= has element byte length %zd, but VECTOR= "
728 "has length %zd",
729 source.ElementBytes(), vector->ElementBytes());
730 }
731 extent = vector->GetDimension(0).Extent();
732 if (extent < trues) {
733 terminator.Crash("PACK: VECTOR= has extent %jd but there are %jd MASK= "
734 "elements that are .TRUE.",
735 static_cast<std::intmax_t>(extent),
736 static_cast<std::intmax_t>(trues));
737 }
738 }
739 AllocateResult(result, source, 1, &extent, terminator, "PACK");
740 SubscriptValue sourceAt[maxRank], resultAt{1};
741 source.GetLowerBounds(sourceAt);
742 if (mask.rank() == 0) {
743 if (IsLogicalElementTrue(mask, nullptr)) {
744 for (SubscriptValue n{trues}; n > 0; --n) {
745 CopyElement(result, &resultAt, source, sourceAt, terminator);
746 ++resultAt;
747 source.IncrementSubscripts(sourceAt);
748 }
749 }
750 } else {
751 SubscriptValue maskAt[maxRank];
752 mask.GetLowerBounds(maskAt);
753 for (std::size_t n{source.Elements()}; n > 0; --n) {
754 if (IsLogicalElementTrue(mask, maskAt)) {
755 CopyElement(result, &resultAt, source, sourceAt, terminator);
756 ++resultAt;
757 }
758 source.IncrementSubscripts(sourceAt);
759 mask.IncrementSubscripts(maskAt);
760 }
761 }
762 if (vector) {
763 SubscriptValue vectorAt{
764 vector->GetDimension(0).LowerBound() + resultAt - 1};
765 for (; resultAt <= extent; ++resultAt, ++vectorAt) {
766 CopyElement(result, &resultAt, *vector, &vectorAt, terminator);
767 }
768 }
769}
770
771// RESHAPE
772// F2018 16.9.163
773void RTDEF(Reshape)(Descriptor &result, const Descriptor &source,
774 const Descriptor &shape, const Descriptor *pad, const Descriptor *order,
775 const char *sourceFile, int line) {
776 // Compute and check the rank of the result.
777 Terminator terminator{sourceFile, line};
778 RUNTIME_CHECK(terminator, shape.rank() == 1);
779 RUNTIME_CHECK(terminator, shape.type().IsInteger());
780 SubscriptValue resultRank{shape.GetDimension(0).Extent()};
781 if (resultRank < 0 || resultRank > static_cast<SubscriptValue>(maxRank)) {
782 terminator.Crash(
783 "RESHAPE: SHAPE= vector length %jd implies a bad result rank",
784 static_cast<std::intmax_t>(resultRank));
785 }
786
787 // Extract and check the shape of the result; compute its element count.
788 SubscriptValue resultExtent[maxRank];
789 std::size_t shapeElementBytes{shape.ElementBytes()};
790 std::size_t resultElements{1};
791 SubscriptValue shapeSubscript{shape.GetDimension(0).LowerBound()};
792 for (int j{0}; j < resultRank; ++j, ++shapeSubscript) {
793 auto extent{GetInt64Safe(
794 shape.Element<char>(&shapeSubscript), shapeElementBytes, terminator)};
795 if (!extent) {
796 terminator.Crash("RESHAPE: value of SHAPE(%d) exceeds 64 bits", j + 1);
797 } else if (*extent < 0) {
798 terminator.Crash("RESHAPE: bad value for SHAPE(%d)=%jd", j + 1,
799 static_cast<std::intmax_t>(*extent));
800 }
801 resultExtent[j] = *extent;
802 resultElements *= resultExtent[j];
803 }
804
805 // Check that there are sufficient elements in the SOURCE=, or that
806 // the optional PAD= argument is present and nonempty.
807 std::size_t elementBytes{source.ElementBytes()};
808 std::size_t sourceElements{source.Elements()};
809 std::size_t padElements{pad ? pad->Elements() : 0};
810 if (resultElements > sourceElements) {
811 if (padElements <= 0) {
812 terminator.Crash(
813 "RESHAPE: not enough elements, need %zd but only have %zd",
814 resultElements, sourceElements);
815 }
816 if (pad->ElementBytes() != elementBytes) {
817 terminator.Crash("RESHAPE: PAD= has element byte length %zd but SOURCE= "
818 "has length %zd",
819 pad->ElementBytes(), elementBytes);
820 }
821 }
822
823 // Extract and check the optional ORDER= argument, which must be a
824 // permutation of [1..resultRank].
825 int dimOrder[maxRank];
826 if (order) {
827 RUNTIME_CHECK(terminator, order->rank() == 1);
828 RUNTIME_CHECK(terminator, order->type().IsInteger());
829 if (order->GetDimension(0).Extent() != resultRank) {
830 terminator.Crash("RESHAPE: the extent of ORDER (%jd) must match the rank"
831 " of the SHAPE (%d)",
832 static_cast<std::intmax_t>(order->GetDimension(0).Extent()),
833 resultRank);
834 }
835 std::uint64_t values{0};
836 SubscriptValue orderSubscript{order->GetDimension(0).LowerBound()};
837 std::size_t orderElementBytes{order->ElementBytes()};
838 for (SubscriptValue j{0}; j < resultRank; ++j, ++orderSubscript) {
839 auto k{GetInt64Safe(order->Element<char>(&orderSubscript),
840 orderElementBytes, terminator)};
841 if (!k) {
842 terminator.Crash("RESHAPE: ORDER element value exceeds 64 bits");
843 } else if (*k < 1 || *k > resultRank || ((values >> *k) & 1)) {
844 terminator.Crash("RESHAPE: bad value for ORDER element (%jd)",
845 static_cast<std::intmax_t>(*k));
846 }
847 values |= std::uint64_t{1} << *k;
848 dimOrder[j] = *k - 1;
849 }
850 } else {
851 for (int j{0}; j < resultRank; ++j) {
852 dimOrder[j] = j;
853 }
854 }
855
856 // Allocate result descriptor
857 AllocateResult(
858 result, source, resultRank, resultExtent, terminator, "RESHAPE");
859
860 // Populate the result's elements.
861 SubscriptValue resultSubscript[maxRank];
862 result.GetLowerBounds(resultSubscript);
863 SubscriptValue sourceSubscript[maxRank];
864 source.GetLowerBounds(sourceSubscript);
865 std::size_t resultElement{0};
866 std::size_t elementsFromSource{std::min(resultElements, sourceElements)};
867 for (; resultElement < elementsFromSource; ++resultElement) {
868 CopyElement(result, resultSubscript, source, sourceSubscript, terminator);
869 source.IncrementSubscripts(sourceSubscript);
870 result.IncrementSubscripts(resultSubscript, dimOrder);
871 }
872 if (resultElement < resultElements) {
873 // Remaining elements come from the optional PAD= argument.
874 SubscriptValue padSubscript[maxRank];
875 pad->GetLowerBounds(padSubscript);
876 for (; resultElement < resultElements; ++resultElement) {
877 CopyElement(result, resultSubscript, *pad, padSubscript, terminator);
878 pad->IncrementSubscripts(padSubscript);
879 result.IncrementSubscripts(resultSubscript, dimOrder);
880 }
881 }
882}
883
884// ShallowCopy
885void RTDEF(ShallowCopy)(Descriptor &result, const Descriptor &source,
886 const char *sourceFile, int line) {
887 Terminator terminator{sourceFile, line};
888 DoShallowCopy<true>(result, source, terminator, "ShallowCopy");
889}
890
891void RTDEF(ShallowCopyDirect)(const Descriptor &result,
892 const Descriptor &source, const char *sourceFile, int line) {
893 Terminator terminator{sourceFile, line};
894 DoShallowCopy<false>(result, source, terminator, "ShallowCopyDirect");
895}
896
897// SPREAD
898void RTDEF(Spread)(Descriptor &result, const Descriptor &source, int dim,
899 std::int64_t ncopies, const char *sourceFile, int line) {
900 Terminator terminator{sourceFile, line};
901 int rank{source.rank() + 1};
902 RUNTIME_CHECK(terminator, rank <= maxRank);
903 if (dim < 1 || dim > rank) {
904 terminator.Crash("SPREAD: DIM=%d argument for rank-%d source array "
905 "must be greater than 1 and less than or equal to %d",
906 dim, rank - 1, rank);
907 }
908 ncopies = std::max<std::int64_t>(ncopies, 0);
909 SubscriptValue extent[maxRank];
910 int k{0};
911 for (int j{0}; j < rank; ++j) {
912 extent[j] = j == dim - 1 ? ncopies : source.GetDimension(k++).Extent();
913 }
914 AllocateResult(result, source, rank, extent, terminator, "SPREAD");
915 SubscriptValue resultAt[maxRank];
916 for (int j{0}; j < rank; ++j) {
917 resultAt[j] = 1;
918 }
919 SubscriptValue &resultDim{resultAt[dim - 1]};
920 SubscriptValue sourceAt[maxRank];
921 source.GetLowerBounds(sourceAt);
922 for (std::size_t n{result.Elements()}; n > 0; n -= ncopies) {
923 for (resultDim = 1; resultDim <= ncopies; ++resultDim) {
924 CopyElement(result, resultAt, source, sourceAt, terminator);
925 }
926 result.IncrementSubscripts(resultAt);
927 source.IncrementSubscripts(sourceAt);
928 }
929}
930
931// TRANSPOSE
932void RTDEF(Transpose)(Descriptor &result, const Descriptor &matrix,
933 const char *sourceFile, int line) {
934 Terminator terminator{sourceFile, line};
935 RUNTIME_CHECK(terminator, matrix.rank() == 2);
936 SubscriptValue extent[2]{
937 matrix.GetDimension(1).Extent(), matrix.GetDimension(0).Extent()};
938 AllocateResult(result, matrix, 2, extent, terminator, "TRANSPOSE");
939 SubscriptValue resultAt[2]{1, 1};
940 SubscriptValue matrixLB[2];
941 matrix.GetLowerBounds(matrixLB);
942 for (std::size_t n{result.Elements()}; n-- > 0;
943 result.IncrementSubscripts(resultAt)) {
944 SubscriptValue matrixAt[2]{
945 matrixLB[0] + resultAt[1] - 1, matrixLB[1] + resultAt[0] - 1};
946 CopyElement(result, resultAt, matrix, matrixAt, terminator);
947 }
948}
949
950// UNPACK
951void RTDEF(Unpack)(Descriptor &result, const Descriptor &vector,
952 const Descriptor &mask, const Descriptor &field, const char *sourceFile,
953 int line) {
954 Terminator terminator{sourceFile, line};
955 RUNTIME_CHECK(terminator, vector.rank() == 1);
956 int rank{mask.rank()};
957 RUNTIME_CHECK(terminator, rank > 0);
958 SubscriptValue extent[maxRank];
959 mask.GetShape(extent);
960 CheckConformability(mask, field, terminator, "UNPACK", "MASK=", "FIELD=");
961 std::size_t elementLen{
962 AllocateResult(result, field, rank, extent, terminator, "UNPACK")};
963 RUNTIME_CHECK(terminator, vector.type() == field.type());
964 if (vector.ElementBytes() != elementLen) {
965 terminator.Crash(
966 "UNPACK: VECTOR= has element byte length %zd but FIELD= has length %zd",
967 vector.ElementBytes(), elementLen);
968 }
969 SubscriptValue resultAt[maxRank], maskAt[maxRank], fieldAt[maxRank],
970 vectorAt{vector.GetDimension(0).LowerBound()};
971 for (int j{0}; j < rank; ++j) {
972 resultAt[j] = 1;
973 }
974 mask.GetLowerBounds(maskAt);
975 field.GetLowerBounds(fieldAt);
976 SubscriptValue vectorElements{vector.GetDimension(0).Extent()};
977 SubscriptValue vectorLeft{vectorElements};
978 for (std::size_t n{result.Elements()}; n-- > 0;) {
979 if (IsLogicalElementTrue(mask, maskAt)) {
980 if (vectorLeft-- == 0) {
981 terminator.Crash(
982 "UNPACK: VECTOR= argument has fewer elements (%d) than "
983 "MASK= has .TRUE. entries",
984 vectorElements);
985 }
986 CopyElement(result, resultAt, vector, &vectorAt, terminator);
987 ++vectorAt;
988 } else {
989 CopyElement(result, resultAt, field, fieldAt, terminator);
990 }
991 result.IncrementSubscripts(resultAt);
992 mask.IncrementSubscripts(maskAt);
993 field.IncrementSubscripts(fieldAt);
994 }
995}
996
997RT_EXT_API_GROUP_END
998} // extern "C"
999} // namespace Fortran::runtime
1000

source code of flang-rt/lib/runtime/transformational.cpp