From 6c81b4a01e4956d70dac1a5733c063014e4fbe38 Mon Sep 17 00:00:00 2001 From: Jean Perier Date: Fri, 22 Apr 2022 09:37:08 +0200 Subject: [PATCH] [flang] Fold transformational bessels when host runtime has bessels Transformational bessel intrinsic functions require the same math runtime as elemental bessel intrinsics. Currently elemental bessels could be folded if f18 was linked with pgmath (cmake -DLIBPGMATH_DIR option). `j0`, `y0`, ... C libm functions were not used because they are not standard C functions: they are Posix extensions. This patch enable: - Using the Posix bessel host runtime functions when available. - folding the transformational bessel using the elemental version. Differential Revision: https://reviews.llvm.org/D124167 --- flang/lib/Evaluate/fold-real.cpp | 35 ++++++++++- flang/lib/Evaluate/intrinsics-library.cpp | 75 ++++++++++++++++++++++- flang/test/Evaluate/folding02.f90 | 16 +++++ 3 files changed, 124 insertions(+), 2 deletions(-) diff --git a/flang/lib/Evaluate/fold-real.cpp b/flang/lib/Evaluate/fold-real.cpp index a1d12c427d41..e915c5576018 100644 --- a/flang/lib/Evaluate/fold-real.cpp +++ b/flang/lib/Evaluate/fold-real.cpp @@ -11,6 +11,38 @@ namespace Fortran::evaluate { +template +static Expr FoldTransformationalBessel( + FunctionRef &&funcRef, FoldingContext &context) { + CHECK(funcRef.arguments().size() == 3); + /// Bessel runtime functions use `int` integer arguments. Convert integer + /// arguments to Int4, any overflow error will be reported during the + /// conversion folding. + using Int4 = Type; + if (auto args{ + GetConstantArguments(context, funcRef.arguments())}) { + const std::string &name{std::get(funcRef.proc().u).name}; + if (auto elementalBessel{GetHostRuntimeWrapper(name)}) { + std::vector> results; + int n1{static_cast( + std::get<0>(*args)->GetScalarValue().value().ToInt64())}; + int n2{static_cast( + std::get<1>(*args)->GetScalarValue().value().ToInt64())}; + Scalar x{std::get<2>(*args)->GetScalarValue().value()}; + for (int i{n1}; i <= n2; ++i) { + results.emplace_back((*elementalBessel)(context, Scalar{i}, x)); + } + return Expr{Constant{ + std::move(results), ConstantSubscripts{std::max(n2 - n1 + 1, 0)}}}; + } else { + context.messages().Say( + "%s(integer(kind=4), real(kind=%d)) cannot be folded on host"_warn_en_US, + name, T::kind); + } + } + return Expr{std::move(funcRef)}; +} + template Expr> FoldIntrinsicFunction( FoldingContext &context, @@ -63,6 +95,8 @@ Expr> FoldIntrinsicFunction( "%s(integer(kind=4), real(kind=%d)) cannot be folded on host"_warn_en_US, name, KIND); } + } else { + return FoldTransformationalBessel(std::move(funcRef), context); } } else if (name == "abs") { // incl. zabs & cdabs // Argument can be complex or real @@ -245,7 +279,6 @@ Expr> FoldIntrinsicFunction( // TODO: dim, dot_product, fraction, matmul, // modulo, norm2, rrspacing, // set_exponent, spacing, transfer, - // bessel_jn (transformational) and bessel_yn (transformational) return Expr{std::move(funcRef)}; } diff --git a/flang/lib/Evaluate/intrinsics-library.cpp b/flang/lib/Evaluate/intrinsics-library.cpp index 8230a595e6a6..ba3b95f3e97a 100644 --- a/flang/lib/Evaluate/intrinsics-library.cpp +++ b/flang/lib/Evaluate/intrinsics-library.cpp @@ -192,7 +192,13 @@ private: // Define host runtime libraries that can be used for folding and // fill their description if they are available. -enum class LibraryVersion { Libm, PgmathFast, PgmathRelaxed, PgmathPrecise }; +enum class LibraryVersion { + Libm, + LibmExtensions, + PgmathFast, + PgmathRelaxed, + PgmathPrecise +}; template struct HostRuntimeLibrary { // When specialized, this class holds a static constexpr table containing // all the HostRuntimeLibrary for functions of library LibraryVersion @@ -277,6 +283,64 @@ struct HostRuntimeLibrary, LibraryVersion::Libm> { static constexpr HostRuntimeMap map{table}; static_assert(map.Verify(), "map must be sorted"); }; +// Note regarding cmath: +// - cmath does not have modulo and erfc_scaled equivalent +// - C++17 defined standard Bessel math functions std::cyl_bessel_j +// and std::cyl_neumann that can be used for Fortran j and y +// bessel functions. However, they are not yet implemented in +// clang libc++ (ok in GNU libstdc++). Instead, the Posix libm +// extensions are used when available below. + +#if _POSIX_C_SOURCE >= 200112L || _XOPEN_SOURCE >= 600 +/// Define libm extensions +/// Bessel functions are defined in POSIX.1-2001. + +template <> struct HostRuntimeLibrary { + using F = FuncPointer; + using FN = FuncPointer; + static constexpr HostRuntimeFunction table[]{ + FolderFactory::Create("bessel_j0"), + FolderFactory::Create("bessel_j1"), + FolderFactory::Create("bessel_jn"), + FolderFactory::Create("bessel_y0"), + FolderFactory::Create("bessel_y1"), + FolderFactory::Create("bessel_yn"), + }; + static constexpr HostRuntimeMap map{table}; + static_assert(map.Verify(), "map must be sorted"); +}; + +template <> struct HostRuntimeLibrary { + using F = FuncPointer; + using FN = FuncPointer; + static constexpr HostRuntimeFunction table[]{ + FolderFactory::Create("bessel_j0"), + FolderFactory::Create("bessel_j1"), + FolderFactory::Create("bessel_jn"), + FolderFactory::Create("bessel_y0"), + FolderFactory::Create("bessel_y1"), + FolderFactory::Create("bessel_yn"), + }; + static constexpr HostRuntimeMap map{table}; + static_assert(map.Verify(), "map must be sorted"); +}; + +template <> +struct HostRuntimeLibrary { + using F = FuncPointer; + using FN = FuncPointer; + static constexpr HostRuntimeFunction table[]{ + FolderFactory::Create("bessel_j0"), + FolderFactory::Create("bessel_j1"), + FolderFactory::Create("bessel_jn"), + FolderFactory::Create("bessel_y0"), + FolderFactory::Create("bessel_y1"), + FolderFactory::Create("bessel_yn"), + }; + static constexpr HostRuntimeMap map{table}; + static_assert(map.Verify(), "map must be sorted"); +}; +#endif /// Define pgmath description #if LINK_WITH_LIBPGMATH @@ -409,6 +473,8 @@ static const HostRuntimeMap *GetHostRuntimeMap( switch (version) { case LibraryVersion::Libm: return GetHostRuntimeMapVersion(resultType); + case LibraryVersion::LibmExtensions: + return GetHostRuntimeMapVersion(resultType); case LibraryVersion::PgmathPrecise: return GetHostRuntimeMapVersion(resultType); case LibraryVersion::PgmathRelaxed: @@ -454,6 +520,13 @@ static const HostRuntimeFunction *SearchHostRuntime(const std::string &name, return hostFunction; } } + if (const auto *map{ + GetHostRuntimeMap(LibraryVersion::LibmExtensions, resultType)}) { + if (const auto *hostFunction{ + SearchInHostRuntimeMap(*map, name, resultType, argTypes)}) { + return hostFunction; + } + } return nullptr; } diff --git a/flang/test/Evaluate/folding02.f90 b/flang/test/Evaluate/folding02.f90 index 7ee365273484..32a465051581 100644 --- a/flang/test/Evaluate/folding02.f90 +++ b/flang/test/Evaluate/folding02.f90 @@ -249,6 +249,22 @@ module m (-0.93219375976297402797143831776338629424571990966796875_8)) TEST_R8(erfc_scaled, erfc_scaled(0.1_8), & 0.89645697996912654392787089818739332258701324462890625_8) + + real(4), parameter :: bessel_jn_transformational(*) = bessel_jn(1,3, 3.2_4) + logical, parameter :: test_bessel_jn_shape = size(bessel_jn_transformational, 1).eq.3 + logical, parameter :: test_bessel_jn_t1 = bessel_jn_transformational(1).eq.bessel_jn(1, 3.2_4) + logical, parameter :: test_bessel_jn_t2 = bessel_jn_transformational(2).eq.bessel_jn(2, 3.2_4) + logical, parameter :: test_bessel_jn_t3 = bessel_jn_transformational(3).eq.bessel_jn(3, 3.2_4) + real(4), parameter :: bessel_jn_empty(*) = bessel_jn(3,1, 3.2_4) + logical, parameter :: test_bessel_jn_empty = size(bessel_jn_empty, 1).eq.0 + + real(4), parameter :: bessel_yn_transformational(*) = bessel_yn(1,3, 1.6_4) + logical, parameter :: test_bessel_yn_shape = size(bessel_yn_transformational, 1).eq.3 + logical, parameter :: test_bessel_yn_t1 = bessel_yn_transformational(1).eq.bessel_yn(1, 1.6_4) + logical, parameter :: test_bessel_yn_t2 = bessel_yn_transformational(2).eq.bessel_yn(2, 1.6_4) + logical, parameter :: test_bessel_yn_t3 = bessel_yn_transformational(3).eq.bessel_yn(3, 1.6_4) + real(4), parameter :: bessel_yn_empty(*) = bessel_yn(3,1, 3.2_4) + logical, parameter :: test_bessel_yn_empty = size(bessel_yn_empty, 1).eq.0 #endif ! Test exponentiation by real or complex folding (it is using host runtime)