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

source code of flang/runtime/transformational.cpp