flang/runtime/transformational.cpp

   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
  26 namespace 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.
  30 class ShiftControl {
  31 public:
  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
  81 private:
  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
  92 static 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(result.raw().base_addr, 0, bytes);
 118   }
 119 }
 120
 121 static 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
 141 template <TypeCategory CAT, int KIND>
 142 static 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
 166 template <TypeCategory CAT, int KIND>
 167 static 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
 219 template <TypeCategory CAT, int KIND>
 220 static 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
 247 template <TypeCategory CAT, int KIND>
 248 static 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
 300 template <TypeCategory CAT, int KIND>
 301 static 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
 326 extern "C" {
 327 RT_EXT_API_GROUP_BEGIN
 328
 329 // BESSEL_JN
 330 // TODO: REAL(2 & 3)
 331 void 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
 338 void 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
 346 void 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
 357 void 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)
 368 void 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
 373 void 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
 379 void 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
 386 void 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)
 394 void 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
 401 void 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
 409 void 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
 420 void 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)
 431 void 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
 436 void 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
 442 void 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
 449 void 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
 456 void 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
 501 void 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
 520 void 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
 600 void 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
 632 void 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
 706 void 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(resultElements, 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
 818 void 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
 852 void 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
 871 void 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
 917 RT_EXT_API_GROUP_END
 918 } // extern "C"
 919 } // namespace Fortran::runtime