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