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