[flang-commits] [flang] bef2bb3 - [flang] Lowering and runtime support for F08 transformational intrinsics: BESSEL_JN and BESSEL_YN

Tarun Prabhu via flang-commits flang-commits at lists.llvm.org
Mon Dec 19 07:00:36 PST 2022


Author: Tarun Prabhu
Date: 2022-12-19T07:59:38-07:00
New Revision: bef2bb34bfadd61cf399a15f096c84980121279f

URL: https://github.com/llvm/llvm-project/commit/bef2bb34bfadd61cf399a15f096c84980121279f
DIFF: https://github.com/llvm/llvm-project/commit/bef2bb34bfadd61cf399a15f096c84980121279f.diff

LOG: [flang] Lowering and runtime support for F08 transformational intrinsics: BESSEL_JN and BESSEL_YN

The runtime implementation uses the recurrence relations

`J(n-1, x) = (2.0 / x) * n * J(n, x) - J(n+1, x)`
`Y(n+1, x) = (2.0 / x) * n * Y(n, x) - Y(n-1, x)`

(see https://dlmf.nist.gov/10.74.iv and https://dlmf.nist.gov/10.6.E1).

Although the standard requires that `N1` and `N2` in `BESSEL_JN(N1, N2, x)`
and `BESSEL_YN(N1, N2, x)` be non-negative, this is not checked in the
runtime functions. This is in keeping with some other compilers which also
return some results when `N1` and/or `N2` are negative.

The special case for `x == 0` is  handled in different runtime functions
for each of `BESSEL_JN` and `BESSEL_YN`. The lowering code checks for this
case and inserts the checks and the appropriate runtime calls in FIR.

The existing tests for the two intrinsics was modified to keep the style
consistent with the additional lowering tests that were added.

Added: 
    

Modified: 
    flang/include/flang/Optimizer/Builder/Runtime/Transformational.h
    flang/include/flang/Runtime/transformational.h
    flang/lib/Lower/IntrinsicCall.cpp
    flang/lib/Optimizer/Builder/Runtime/Transformational.cpp
    flang/runtime/transformational.cpp
    flang/test/Lower/Intrinsics/bessel_jn.f90
    flang/test/Lower/Intrinsics/bessel_yn.f90
    flang/unittests/Optimizer/Builder/Runtime/TransformationalTest.cpp
    flang/unittests/Runtime/Transformational.cpp

Removed: 
    


################################################################################
diff  --git a/flang/include/flang/Optimizer/Builder/Runtime/Transformational.h b/flang/include/flang/Optimizer/Builder/Runtime/Transformational.h
index 1657dff0fd52..f084e2eb7ae3 100644
--- a/flang/include/flang/Optimizer/Builder/Runtime/Transformational.h
+++ b/flang/include/flang/Optimizer/Builder/Runtime/Transformational.h
@@ -19,6 +19,22 @@ class FirOpBuilder;
 
 namespace fir::runtime {
 
+void genBesselJn(fir::FirOpBuilder &builder, mlir::Location loc,
+                 mlir::Value resultBox, mlir::Value n1, mlir::Value n2,
+                 mlir::Value x, mlir::Value bn2, mlir::Value bn2_1);
+
+void genBesselJnX0(fir::FirOpBuilder &builder, mlir::Location loc,
+                   mlir::Type xTy, mlir::Value resultBox, mlir::Value n1,
+                   mlir::Value n2);
+
+void genBesselYn(fir::FirOpBuilder &builder, mlir::Location loc,
+                 mlir::Value resultBox, mlir::Value n1, mlir::Value n2,
+                 mlir::Value x, mlir::Value bn1, mlir::Value bn1_1);
+
+void genBesselYnX0(fir::FirOpBuilder &builder, mlir::Location loc,
+                   mlir::Type xTy, mlir::Value resultBox, mlir::Value n1,
+                   mlir::Value n2);
+
 void genCshift(fir::FirOpBuilder &builder, mlir::Location loc,
                mlir::Value resultBox, mlir::Value arrayBox,
                mlir::Value shiftBox, mlir::Value dimBox);

diff  --git a/flang/include/flang/Runtime/transformational.h b/flang/include/flang/Runtime/transformational.h
index 21a44184b622..8101c73ba32c 100644
--- a/flang/include/flang/Runtime/transformational.h
+++ b/flang/include/flang/Runtime/transformational.h
@@ -17,7 +17,9 @@
 #ifndef FORTRAN_RUNTIME_TRANSFORMATIONAL_H_
 #define FORTRAN_RUNTIME_TRANSFORMATIONAL_H_
 
+#include "flang/Runtime/cpp-type.h"
 #include "flang/Runtime/entry-names.h"
+#include "flang/Runtime/float128.h"
 #include <cinttypes>
 
 namespace Fortran::runtime {
@@ -31,6 +33,98 @@ void RTNAME(Reshape)(Descriptor &result, const Descriptor &source,
     const Descriptor *order = nullptr, const char *sourceFile = nullptr,
     int line = 0);
 
+void RTNAME(BesselJn_2)(Descriptor &result, int32_t n1, int32_t n2, float x,
+    float bn2, float bn2_1, const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselJn_3)(Descriptor &result, int32_t n1, int32_t n2, float x,
+    float bn2, float bn2_1, const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselJn_4)(Descriptor &result, int32_t n1, int32_t n2, float x,
+    float bn2, float bn2_1, const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselJn_8)(Descriptor &result, int32_t n1, int32_t n2, double x,
+    double bn2, double bn2_1, const char *sourceFile = nullptr, int line = 0);
+
+#if LDBL_MANT_DIG == 64
+void RTNAME(BesselJn_10)(Descriptor &result, int32_t n1, int32_t n2,
+    long double x, long double bn2, long double bn2_1,
+    const char *sourceFile = nullptr, int line = 0);
+#endif
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+void RTNAME(BesselJn_16)(Descriptor &result, int32_t n1, int32_t n2,
+    CppFloat128Type x, CppFloat128Type bn2, CppFloat128Type bn2_1,
+    const char *sourceFile = nullptr, int line = 0);
+#endif
+
+void RTNAME(BesselJnX0_2)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselJnX0_3)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselJnX0_4)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselJnX0_8)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+
+#if LDBL_MANT_DIG == 64
+void RTNAME(BesselJnX0_10)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+#endif
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+void RTNAME(BesselJnX0_16)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+#endif
+
+void RTNAME(BesselYn_2)(Descriptor &result, int32_t n1, int32_t n2, float x,
+    float bn1, float bn1_1, const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselYn_3)(Descriptor &result, int32_t n1, int32_t n2, float x,
+    float bn1, float bn1_1, const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselYn_4)(Descriptor &result, int32_t n1, int32_t n2, float x,
+    float bn1, float bn1_1, const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselYn_8)(Descriptor &result, int32_t n1, int32_t n2, double x,
+    double bn1, double bn1_1, const char *sourceFile = nullptr, int line = 0);
+
+#if LDBL_MANT_DIG == 64
+void RTNAME(BesselYn_10)(Descriptor &result, int32_t n1, int32_t n2,
+    long double x, long double bn1, long double bn1_1,
+    const char *sourceFile = nullptr, int line = 0);
+#endif
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+void RTNAME(BesselYn_16)(Descriptor &result, int32_t n1, int32_t n2,
+    CppFloat128Type x, CppFloat128Type bn1, CppFloat128Type bn1_1,
+    const char *sourceFile = nullptr, int line = 0);
+#endif
+
+void RTNAME(BesselYnX0_2)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselYnX0_3)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselYnX0_4)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselYnX0_8)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+
+#if LDBL_MANT_DIG == 64
+void RTNAME(BesselYnX0_10)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+#endif
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+void RTNAME(BesselYnX0_16)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+#endif
+
 void RTNAME(Cshift)(Descriptor &result, const Descriptor &source,
     const Descriptor &shift, int dim = 1, const char *sourceFile = nullptr,
     int line = 0);

diff  --git a/flang/lib/Lower/IntrinsicCall.cpp b/flang/lib/Lower/IntrinsicCall.cpp
index 65bb310ba90c..cdccd451b9a6 100644
--- a/flang/lib/Lower/IntrinsicCall.cpp
+++ b/flang/lib/Lower/IntrinsicCall.cpp
@@ -462,7 +462,10 @@ struct IntrinsicLibrary {
       genCommandArgumentCount(mlir::Type, llvm::ArrayRef<fir::ExtendedValue>);
   fir::ExtendedValue genAssociated(mlir::Type,
                                    llvm::ArrayRef<fir::ExtendedValue>);
-
+  fir::ExtendedValue genBesselJn(mlir::Type,
+                                 llvm::ArrayRef<fir::ExtendedValue>);
+  fir::ExtendedValue genBesselYn(mlir::Type,
+                                 llvm::ArrayRef<fir::ExtendedValue>);
   /// Lower a bitwise comparison intrinsic using the given comparator.
   template <mlir::arith::CmpIPredicate pred>
   mlir::Value genBitwiseCompare(mlir::Type resultType,
@@ -732,6 +735,14 @@ static constexpr IntrinsicHandler handlers[]{
      &I::genAssociated,
      {{{"pointer", asInquired}, {"target", asInquired}}},
      /*isElemental=*/false},
+    {"bessel_jn",
+     &I::genBesselJn,
+     {{{"n1", asValue}, {"n2", asValue}, {"x", asValue}}},
+     /*isElemental=*/false},
+    {"bessel_yn",
+     &I::genBesselYn,
+     {{{"n1", asValue}, {"n2", asValue}, {"x", asValue}}},
+     /*isElemental=*/false},
     {"bge", &I::genBitwiseCompare<mlir::arith::CmpIPredicate::uge>},
     {"bgt", &I::genBitwiseCompare<mlir::arith::CmpIPredicate::ugt>},
     {"ble", &I::genBitwiseCompare<mlir::arith::CmpIPredicate::ule>},
@@ -2527,6 +2538,196 @@ IntrinsicLibrary::genAssociated(mlir::Type resultType,
   return Fortran::lower::genAssociated(builder, loc, pointerBox, targetBox);
 }
 
+// BESSEL_JN
+fir::ExtendedValue
+IntrinsicLibrary::genBesselJn(mlir::Type resultType,
+                              llvm::ArrayRef<fir::ExtendedValue> args) {
+  assert(args.size() == 2 || args.size() == 3);
+
+  mlir::Value x = fir::getBase(args.back());
+
+  if (args.size() == 2) {
+    mlir::Value n = fir::getBase(args[0]);
+
+    return genRuntimeCall("bessel_jn", resultType, {n, x});
+  } else {
+    mlir::Value n1 = fir::getBase(args[0]);
+    mlir::Value n2 = fir::getBase(args[1]);
+
+    mlir::Type intTy = n1.getType();
+    mlir::Type floatTy = x.getType();
+    mlir::Value zero = builder.createRealZeroConstant(loc, floatTy);
+    mlir::Value one = builder.createIntegerConstant(loc, intTy, 1);
+
+    mlir::Type resultArrayType = builder.getVarLenSeqTy(resultType, 1);
+    fir::MutableBoxValue resultMutableBox =
+        fir::factory::createTempMutableBox(builder, loc, resultArrayType);
+    mlir::Value resultBox =
+        fir::factory::getMutableIRBox(builder, loc, resultMutableBox);
+
+    mlir::Value cmpXEq0 = builder.create<mlir::arith::CmpFOp>(
+        loc, mlir::arith::CmpFPredicate::UEQ, x, zero);
+    mlir::Value cmpN1LtN2 = builder.create<mlir::arith::CmpIOp>(
+        loc, mlir::arith::CmpIPredicate::slt, n1, n2);
+    mlir::Value cmpN1EqN2 = builder.create<mlir::arith::CmpIOp>(
+        loc, mlir::arith::CmpIPredicate::eq, n1, n2);
+
+    auto genXEq0 = [&]() {
+      fir::runtime::genBesselJnX0(builder, loc, floatTy, resultBox, n1, n2);
+    };
+
+    auto genN1LtN2 = [&]() {
+      // The runtime generates the values in the range using a backward
+      // recursion from n2 to n1. (see https://dlmf.nist.gov/10.74.iv and
+      // https://dlmf.nist.gov/10.6.E1). When n1 < n2, this requires
+      // the values of BESSEL_JN(n2) and BESSEL_JN(n2 - 1) since they
+      // are the anchors of the recursion.
+      mlir::Value n2_1 = builder.create<mlir::arith::SubIOp>(loc, n2, one);
+      mlir::Value bn2 = genRuntimeCall("bessel_jn", resultType, {n2, x});
+      mlir::Value bn2_1 = genRuntimeCall("bessel_jn", resultType, {n2_1, x});
+      fir::runtime::genBesselJn(builder, loc, resultBox, n1, n2, x, bn2, bn2_1);
+    };
+
+    auto genN1EqN2 = [&]() {
+      // When n1 == n2, only BESSEL_JN(n2) is needed.
+      mlir::Value bn2 = genRuntimeCall("bessel_jn", resultType, {n2, x});
+      fir::runtime::genBesselJn(builder, loc, resultBox, n1, n2, x, bn2, zero);
+    };
+
+    auto genN1GtN2 = [&]() {
+      // The standard requires n1 <= n2. However, we still need to allocate
+      // a zero-length array and return it when n1 > n2, so we do need to call
+      // the runtime function.
+      fir::runtime::genBesselJn(builder, loc, resultBox, n1, n2, x, zero, zero);
+    };
+
+    auto genN1GeN2 = [&] {
+      builder.genIfThenElse(loc, cmpN1EqN2)
+          .genThen(genN1EqN2)
+          .genElse(genN1GtN2)
+          .end();
+    };
+
+    auto genXNeq0 = [&]() {
+      builder.genIfThenElse(loc, cmpN1LtN2)
+          .genThen(genN1LtN2)
+          .genElse(genN1GeN2)
+          .end();
+    };
+
+    builder.genIfThenElse(loc, cmpXEq0)
+        .genThen(genXEq0)
+        .genElse(genXNeq0)
+        .end();
+
+    fir::ExtendedValue res =
+        fir::factory::genMutableBoxRead(builder, loc, resultMutableBox);
+    return res.match(
+        [&](const fir::ArrayBoxValue &box) -> fir::ExtendedValue {
+          addCleanUpForTemp(loc, box.getAddr());
+          return box;
+        },
+        [&](const auto &) -> fir::ExtendedValue {
+          fir::emitFatalError(loc, "unexpected result for BESSEL_JN");
+        });
+  }
+}
+
+// BESSEL_YN
+fir::ExtendedValue
+IntrinsicLibrary::genBesselYn(mlir::Type resultType,
+                              llvm::ArrayRef<fir::ExtendedValue> args) {
+  assert(args.size() == 2 || args.size() == 3);
+
+  mlir::Value x = fir::getBase(args.back());
+
+  if (args.size() == 2) {
+    mlir::Value n = fir::getBase(args[0]);
+
+    return genRuntimeCall("bessel_yn", resultType, {n, x});
+  } else {
+    mlir::Value n1 = fir::getBase(args[0]);
+    mlir::Value n2 = fir::getBase(args[1]);
+
+    mlir::Type floatTy = x.getType();
+    mlir::Type intTy = n1.getType();
+    mlir::Value zero = builder.createRealZeroConstant(loc, floatTy);
+    mlir::Value one = builder.createIntegerConstant(loc, intTy, 1);
+
+    mlir::Type resultArrayType = builder.getVarLenSeqTy(resultType, 1);
+    fir::MutableBoxValue resultMutableBox =
+        fir::factory::createTempMutableBox(builder, loc, resultArrayType);
+    mlir::Value resultBox =
+        fir::factory::getMutableIRBox(builder, loc, resultMutableBox);
+
+    mlir::Value cmpXEq0 = builder.create<mlir::arith::CmpFOp>(
+        loc, mlir::arith::CmpFPredicate::UEQ, x, zero);
+    mlir::Value cmpN1LtN2 = builder.create<mlir::arith::CmpIOp>(
+        loc, mlir::arith::CmpIPredicate::slt, n1, n2);
+    mlir::Value cmpN1EqN2 = builder.create<mlir::arith::CmpIOp>(
+        loc, mlir::arith::CmpIPredicate::eq, n1, n2);
+
+    auto genXEq0 = [&]() {
+      fir::runtime::genBesselYnX0(builder, loc, floatTy, resultBox, n1, n2);
+    };
+
+    auto genN1LtN2 = [&]() {
+      // The runtime generates the values in the range using a forward
+      // recursion from n1 to n2. (see https://dlmf.nist.gov/10.74.iv and
+      // https://dlmf.nist.gov/10.6.E1). When n1 < n2, this requires
+      // the values of BESSEL_YN(n1) and BESSEL_YN(n1 + 1) since they
+      // are the anchors of the recursion.
+      mlir::Value n1_1 = builder.create<mlir::arith::AddIOp>(loc, n1, one);
+      mlir::Value bn1 = genRuntimeCall("bessel_yn", resultType, {n1, x});
+      mlir::Value bn1_1 = genRuntimeCall("bessel_yn", resultType, {n1_1, x});
+      fir::runtime::genBesselYn(builder, loc, resultBox, n1, n2, x, bn1, bn1_1);
+    };
+
+    auto genN1EqN2 = [&]() {
+      // When n1 == n2, only BESSEL_YN(n1) is needed.
+      mlir::Value bn1 = genRuntimeCall("bessel_yn", resultType, {n1, x});
+      fir::runtime::genBesselYn(builder, loc, resultBox, n1, n2, x, bn1, zero);
+    };
+
+    auto genN1GtN2 = [&]() {
+      // The standard requires n1 <= n2. However, we still need to allocate
+      // a zero-length array and return it when n1 > n2, so we do need to call
+      // the runtime function.
+      fir::runtime::genBesselYn(builder, loc, resultBox, n1, n2, x, zero, zero);
+    };
+
+    auto genN1GeN2 = [&] {
+      builder.genIfThenElse(loc, cmpN1EqN2)
+          .genThen(genN1EqN2)
+          .genElse(genN1GtN2)
+          .end();
+    };
+
+    auto genXNeq0 = [&]() {
+      builder.genIfThenElse(loc, cmpN1LtN2)
+          .genThen(genN1LtN2)
+          .genElse(genN1GeN2)
+          .end();
+    };
+
+    builder.genIfThenElse(loc, cmpXEq0)
+        .genThen(genXEq0)
+        .genElse(genXNeq0)
+        .end();
+
+    fir::ExtendedValue res =
+        fir::factory::genMutableBoxRead(builder, loc, resultMutableBox);
+    return res.match(
+        [&](const fir::ArrayBoxValue &box) -> fir::ExtendedValue {
+          addCleanUpForTemp(loc, box.getAddr());
+          return box;
+        },
+        [&](const auto &) -> fir::ExtendedValue {
+          fir::emitFatalError(loc, "unexpected result for BESSEL_YN");
+        });
+  }
+}
+
 // BGE, BGT, BLE, BLT
 template <mlir::arith::CmpIPredicate pred>
 mlir::Value

diff  --git a/flang/lib/Optimizer/Builder/Runtime/Transformational.cpp b/flang/lib/Optimizer/Builder/Runtime/Transformational.cpp
index 0b8e14555d1d..c51f38672d25 100644
--- a/flang/lib/Optimizer/Builder/Runtime/Transformational.cpp
+++ b/flang/lib/Optimizer/Builder/Runtime/Transformational.cpp
@@ -19,6 +19,256 @@
 
 using namespace Fortran::runtime;
 
+/// Placeholder for real*10 version of BesselJn intrinsic.
+struct ForcedBesselJn_10 {
+  static constexpr const char *name = ExpandAndQuoteKey(RTNAME(BesselJn_10));
+  static constexpr fir::runtime::FuncTypeBuilderFunc getTypeModel() {
+    return [](mlir::MLIRContext *ctx) {
+      auto ty = mlir::FloatType::getF80(ctx);
+      auto boxTy =
+          fir::runtime::getModel<Fortran::runtime::Descriptor &>()(ctx);
+      auto strTy = fir::ReferenceType::get(mlir::IntegerType::get(ctx, 8));
+      auto intTy = mlir::IntegerType::get(ctx, 32);
+      auto noneTy = mlir::NoneType::get(ctx);
+      return mlir::FunctionType::get(
+          ctx, {boxTy, intTy, intTy, ty, ty, ty, strTy, intTy}, {noneTy});
+    };
+  }
+};
+
+/// Placeholder for real*16 version of BesselJn intrinsic.
+struct ForcedBesselJn_16 {
+  static constexpr const char *name = ExpandAndQuoteKey(RTNAME(BesselJn_16));
+  static constexpr fir::runtime::FuncTypeBuilderFunc getTypeModel() {
+    return [](mlir::MLIRContext *ctx) {
+      auto ty = mlir::FloatType::getF128(ctx);
+      auto boxTy =
+          fir::runtime::getModel<Fortran::runtime::Descriptor &>()(ctx);
+      auto strTy = fir::ReferenceType::get(mlir::IntegerType::get(ctx, 8));
+      auto intTy = mlir::IntegerType::get(ctx, 32);
+      auto noneTy = mlir::NoneType::get(ctx);
+      return mlir::FunctionType::get(
+          ctx, {boxTy, intTy, intTy, ty, ty, ty, strTy, intTy}, {noneTy});
+    };
+  }
+};
+
+/// Placeholder for real*10 version of BesselJn intrinsic when `x == 0.0`.
+struct ForcedBesselJnX0_10 {
+  static constexpr const char *name = ExpandAndQuoteKey(RTNAME(BesselJnX0_10));
+  static constexpr fir::runtime::FuncTypeBuilderFunc getTypeModel() {
+    return [](mlir::MLIRContext *ctx) {
+      auto boxTy =
+          fir::runtime::getModel<Fortran::runtime::Descriptor &>()(ctx);
+      auto strTy = fir::ReferenceType::get(mlir::IntegerType::get(ctx, 8));
+      auto intTy = mlir::IntegerType::get(ctx, 32);
+      auto noneTy = mlir::NoneType::get(ctx);
+      return mlir::FunctionType::get(ctx, {boxTy, intTy, intTy, strTy, intTy},
+                                     {noneTy});
+    };
+  }
+};
+
+/// Placeholder for real*16 version of BesselJn intrinsic when `x == 0.0`.
+struct ForcedBesselJnX0_16 {
+  static constexpr const char *name = ExpandAndQuoteKey(RTNAME(BesselJnX0_16));
+  static constexpr fir::runtime::FuncTypeBuilderFunc getTypeModel() {
+    return [](mlir::MLIRContext *ctx) {
+      auto boxTy =
+          fir::runtime::getModel<Fortran::runtime::Descriptor &>()(ctx);
+      auto strTy = fir::ReferenceType::get(mlir::IntegerType::get(ctx, 8));
+      auto intTy = mlir::IntegerType::get(ctx, 32);
+      auto noneTy = mlir::NoneType::get(ctx);
+      return mlir::FunctionType::get(ctx, {boxTy, intTy, intTy, strTy, intTy},
+                                     {noneTy});
+    };
+  }
+};
+
+/// Placeholder for real*10 version of BesselYn intrinsic.
+struct ForcedBesselYn_10 {
+  static constexpr const char *name = ExpandAndQuoteKey(RTNAME(BesselYn_10));
+  static constexpr fir::runtime::FuncTypeBuilderFunc getTypeModel() {
+    return [](mlir::MLIRContext *ctx) {
+      auto ty = mlir::FloatType::getF80(ctx);
+      auto boxTy =
+          fir::runtime::getModel<Fortran::runtime::Descriptor &>()(ctx);
+      auto strTy = fir::ReferenceType::get(mlir::IntegerType::get(ctx, 8));
+      auto intTy = mlir::IntegerType::get(ctx, 32);
+      auto noneTy = mlir::NoneType::get(ctx);
+      return mlir::FunctionType::get(
+          ctx, {boxTy, intTy, intTy, ty, ty, ty, strTy, intTy}, {noneTy});
+    };
+  }
+};
+
+/// Placeholder for real*16 version of BesselYn intrinsic.
+struct ForcedBesselYn_16 {
+  static constexpr const char *name = ExpandAndQuoteKey(RTNAME(BesselYn_16));
+  static constexpr fir::runtime::FuncTypeBuilderFunc getTypeModel() {
+    return [](mlir::MLIRContext *ctx) {
+      auto ty = mlir::FloatType::getF128(ctx);
+      auto boxTy =
+          fir::runtime::getModel<Fortran::runtime::Descriptor &>()(ctx);
+      auto strTy = fir::ReferenceType::get(mlir::IntegerType::get(ctx, 8));
+      auto intTy = mlir::IntegerType::get(ctx, 32);
+      auto noneTy = mlir::NoneType::get(ctx);
+      return mlir::FunctionType::get(
+          ctx, {boxTy, intTy, intTy, ty, ty, ty, strTy, intTy}, {noneTy});
+    };
+  }
+};
+
+/// Placeholder for real*10 version of BesselYn intrinsic when `x == 0.0`.
+struct ForcedBesselYnX0_10 {
+  static constexpr const char *name = ExpandAndQuoteKey(RTNAME(BesselYnX0_10));
+  static constexpr fir::runtime::FuncTypeBuilderFunc getTypeModel() {
+    return [](mlir::MLIRContext *ctx) {
+      auto boxTy =
+          fir::runtime::getModel<Fortran::runtime::Descriptor &>()(ctx);
+      auto strTy = fir::ReferenceType::get(mlir::IntegerType::get(ctx, 8));
+      auto intTy = mlir::IntegerType::get(ctx, 32);
+      auto noneTy = mlir::NoneType::get(ctx);
+      return mlir::FunctionType::get(ctx, {boxTy, intTy, intTy, strTy, intTy},
+                                     {noneTy});
+    };
+  }
+};
+
+/// Placeholder for real*16 version of BesselYn intrinsic when `x == 0.0`.
+struct ForcedBesselYnX0_16 {
+  static constexpr const char *name = ExpandAndQuoteKey(RTNAME(BesselYnX0_16));
+  static constexpr fir::runtime::FuncTypeBuilderFunc getTypeModel() {
+    return [](mlir::MLIRContext *ctx) {
+      auto boxTy =
+          fir::runtime::getModel<Fortran::runtime::Descriptor &>()(ctx);
+      auto strTy = fir::ReferenceType::get(mlir::IntegerType::get(ctx, 8));
+      auto intTy = mlir::IntegerType::get(ctx, 32);
+      auto noneTy = mlir::NoneType::get(ctx);
+      return mlir::FunctionType::get(ctx, {boxTy, intTy, intTy, strTy, intTy},
+                                     {noneTy});
+    };
+  }
+};
+
+/// Generate call to `BesselJn` intrinsic.
+void fir::runtime::genBesselJn(fir::FirOpBuilder &builder, mlir::Location loc,
+                               mlir::Value resultBox, mlir::Value n1,
+                               mlir::Value n2, mlir::Value x, mlir::Value bn2,
+                               mlir::Value bn2_1) {
+  mlir::func::FuncOp func;
+  auto xTy = x.getType();
+
+  if (xTy.isF16() || xTy.isBF16())
+    TODO(loc, "half-precision BESSEL_JN");
+  else if (xTy.isF32())
+    func = fir::runtime::getRuntimeFunc<mkRTKey(BesselJn_4)>(loc, builder);
+  else if (xTy.isF64())
+    func = fir::runtime::getRuntimeFunc<mkRTKey(BesselJn_8)>(loc, builder);
+  else if (xTy.isF80())
+    func = fir::runtime::getRuntimeFunc<ForcedBesselJn_10>(loc, builder);
+  else if (xTy.isF128())
+    func = fir::runtime::getRuntimeFunc<ForcedBesselJn_16>(loc, builder);
+  else
+    fir::emitFatalError(loc, "invalid type in BESSEL_JN");
+
+  auto fTy = func.getFunctionType();
+  auto sourceFile = fir::factory::locationToFilename(builder, loc);
+  auto sourceLine =
+      fir::factory::locationToLineNo(builder, loc, fTy.getInput(7));
+  auto args =
+      fir::runtime::createArguments(builder, loc, fTy, resultBox, n1, n2, x,
+                                    bn2, bn2_1, sourceFile, sourceLine);
+  builder.create<fir::CallOp>(loc, func, args);
+}
+
+/// Generate call to `BesselJn` intrinsic. This is used when `x == 0.0`.
+void fir::runtime::genBesselJnX0(fir::FirOpBuilder &builder, mlir::Location loc,
+                                 mlir::Type xTy, mlir::Value resultBox,
+                                 mlir::Value n1, mlir::Value n2) {
+  mlir::func::FuncOp func;
+
+  if (xTy.isF16() || xTy.isBF16())
+    TODO(loc, "half-precision BESSEL_JN");
+  else if (xTy.isF32())
+    func = fir::runtime::getRuntimeFunc<mkRTKey(BesselJnX0_4)>(loc, builder);
+  else if (xTy.isF64())
+    func = fir::runtime::getRuntimeFunc<mkRTKey(BesselJnX0_8)>(loc, builder);
+  else if (xTy.isF80())
+    func = fir::runtime::getRuntimeFunc<ForcedBesselJnX0_10>(loc, builder);
+  else if (xTy.isF128())
+    func = fir::runtime::getRuntimeFunc<ForcedBesselJnX0_16>(loc, builder);
+  else
+    fir::emitFatalError(loc, "invalid type in BESSEL_JN");
+
+  auto fTy = func.getFunctionType();
+  auto sourceFile = fir::factory::locationToFilename(builder, loc);
+  auto sourceLine =
+      fir::factory::locationToLineNo(builder, loc, fTy.getInput(4));
+  auto args = fir::runtime::createArguments(builder, loc, fTy, resultBox, n1,
+                                            n2, sourceFile, sourceLine);
+  builder.create<fir::CallOp>(loc, func, args);
+}
+
+/// Generate call to `BesselYn` intrinsic.
+void fir::runtime::genBesselYn(fir::FirOpBuilder &builder, mlir::Location loc,
+                               mlir::Value resultBox, mlir::Value n1,
+                               mlir::Value n2, mlir::Value x, mlir::Value bn1,
+                               mlir::Value bn1_1) {
+  mlir::func::FuncOp func;
+  auto xTy = x.getType();
+
+  if (xTy.isF16() || xTy.isBF16())
+    TODO(loc, "half-precision BESSEL_YN");
+  else if (xTy.isF32())
+    func = fir::runtime::getRuntimeFunc<mkRTKey(BesselYn_4)>(loc, builder);
+  else if (xTy.isF64())
+    func = fir::runtime::getRuntimeFunc<mkRTKey(BesselYn_8)>(loc, builder);
+  else if (xTy.isF80())
+    func = fir::runtime::getRuntimeFunc<ForcedBesselYn_10>(loc, builder);
+  else if (xTy.isF128())
+    func = fir::runtime::getRuntimeFunc<ForcedBesselYn_16>(loc, builder);
+  else
+    fir::emitFatalError(loc, "invalid type in BESSEL_YN");
+
+  auto fTy = func.getFunctionType();
+  auto sourceFile = fir::factory::locationToFilename(builder, loc);
+  auto sourceLine =
+      fir::factory::locationToLineNo(builder, loc, fTy.getInput(7));
+  auto args =
+      fir::runtime::createArguments(builder, loc, fTy, resultBox, n1, n2, x,
+                                    bn1, bn1_1, sourceFile, sourceLine);
+  builder.create<fir::CallOp>(loc, func, args);
+}
+
+/// Generate call to `BesselYn` intrinsic. This is used when `x == 0.0`.
+void fir::runtime::genBesselYnX0(fir::FirOpBuilder &builder, mlir::Location loc,
+                                 mlir::Type xTy, mlir::Value resultBox,
+                                 mlir::Value n1, mlir::Value n2) {
+  mlir::func::FuncOp func;
+
+  if (xTy.isF16() || xTy.isBF16())
+    TODO(loc, "half-precision BESSEL_YN");
+  else if (xTy.isF32())
+    func = fir::runtime::getRuntimeFunc<mkRTKey(BesselYnX0_4)>(loc, builder);
+  else if (xTy.isF64())
+    func = fir::runtime::getRuntimeFunc<mkRTKey(BesselYnX0_8)>(loc, builder);
+  else if (xTy.isF80())
+    func = fir::runtime::getRuntimeFunc<ForcedBesselYnX0_10>(loc, builder);
+  else if (xTy.isF128())
+    func = fir::runtime::getRuntimeFunc<ForcedBesselYnX0_16>(loc, builder);
+  else
+    fir::emitFatalError(loc, "invalid type in BESSEL_YN");
+
+  auto fTy = func.getFunctionType();
+  auto sourceFile = fir::factory::locationToFilename(builder, loc);
+  auto sourceLine =
+      fir::factory::locationToLineNo(builder, loc, fTy.getInput(4));
+  auto args = fir::runtime::createArguments(builder, loc, fTy, resultBox, n1,
+                                            n2, sourceFile, sourceLine);
+  builder.create<fir::CallOp>(loc, func, args);
+}
+
 /// Generate call to Cshift intrinsic
 void fir::runtime::genCshift(fir::FirOpBuilder &builder, mlir::Location loc,
                              mlir::Value resultBox, mlir::Value arrayBox,

diff  --git a/flang/runtime/transformational.cpp b/flang/runtime/transformational.cpp
index 9a0168167886..00bf12f7a7be 100644
--- a/flang/runtime/transformational.cpp
+++ b/flang/runtime/transformational.cpp
@@ -133,8 +133,319 @@ static inline std::size_t AllocateResult(Descriptor &result,
   return elementLen;
 }
 
+template <TypeCategory CAT, int KIND>
+static inline std::size_t AllocateBesselResult(Descriptor &result, int32_t n1,
+    int32_t n2, Terminator &terminator, const char *function) {
+  int rank{1};
+  SubscriptValue extent[maxRank];
+  for (int j{0}; j < maxRank; j++) {
+    extent[j] = 0;
+  }
+  if (n1 <= n2) {
+    extent[0] = n2 - n1 + 1;
+  }
+
+  std::size_t elementLen{Descriptor::BytesFor(CAT, KIND)};
+  result.Establish(TypeCode{CAT, KIND}, elementLen, nullptr, rank, extent,
+      CFI_attribute_allocatable, false);
+  for (int j{0}; j < rank; ++j) {
+    result.GetDimension(j).SetBounds(1, extent[j]);
+  }
+  if (int stat{result.Allocate()}) {
+    terminator.Crash(
+        "%s: Could not allocate memory for result (stat=%d)", function, stat);
+  }
+  return elementLen;
+}
+
+template <TypeCategory CAT, int KIND>
+static inline void DoBesselJn(Descriptor &result, int32_t n1, int32_t n2,
+    CppTypeFor<CAT, KIND> x, CppTypeFor<CAT, KIND> bn2,
+    CppTypeFor<CAT, KIND> bn2_1, const char *sourceFile, int line) {
+  Terminator terminator{sourceFile, line};
+  AllocateBesselResult<CAT, KIND>(result, n1, n2, terminator, "BESSEL_JN");
+
+  // The standard requires that n1 and n2 be non-negative. However, some other
+  // compilers generate results even when n1 and/or n2 are negative. For now,
+  // we also do not enforce the non-negativity constraint.
+  if (n2 < n1) {
+    return;
+  }
+
+  SubscriptValue at[maxRank];
+  for (int j{0}; j < maxRank; ++j) {
+    at[j] = 0;
+  }
+
+  // if n2 >= n1, there will be at least one element in the result.
+  at[0] = n2 - n1 + 1;
+  *result.Element<CppTypeFor<CAT, KIND>>(at) = bn2;
+
+  if (n2 == n1) {
+    return;
+  }
+
+  at[0] = n2 - n1;
+  *result.Element<CppTypeFor<CAT, KIND>>(at) = bn2_1;
+
+  // Bessel functions of the first kind are stable for a backward recursion
+  // (see https://dlmf.nist.gov/10.74.iv and https://dlmf.nist.gov/10.6.E1).
+  //
+  //     J(n-1, x) = (2.0 / x) * n * J(n, x) - J(n+1, x)
+  //
+  // which is equivalent to
+  //
+  //     J(n, x) = (2.0 / x) * (n + 1) * J(n+1, x) - J(n+2, x)
+  //
+  CppTypeFor<CAT, KIND> bn_2 = bn2;
+  CppTypeFor<CAT, KIND> bn_1 = bn2_1;
+  CppTypeFor<CAT, KIND> twoOverX = 2.0 / x;
+  for (int n{n2 - 2}; n >= n1; --n) {
+    auto bn = twoOverX * (n + 1) * bn_1 - bn_2;
+
+    at[0] = n - n1 + 1;
+    *result.Element<CppTypeFor<CAT, KIND>>(at) = bn;
+
+    bn_2 = bn_1;
+    bn_1 = bn;
+  }
+}
+
+template <TypeCategory CAT, int KIND>
+static inline void DoBesselJnX0(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile, int line) {
+  Terminator terminator{sourceFile, line};
+  AllocateBesselResult<CAT, KIND>(result, n1, n2, terminator, "BESSEL_JN");
+
+  // The standard requires that n1 and n2 be non-negative. However, some other
+  // compilers generate results even when n1 and/or n2 are negative. For now,
+  // we also do not enforce the non-negativity constraint.
+  if (n2 < n1) {
+    return;
+  }
+
+  SubscriptValue at[maxRank];
+  for (int j{0}; j < maxRank; ++j) {
+    at[j] = 0;
+  }
+
+  // J(0, 0.0) = 1.0, when n == 0.
+  // J(n, 0.0) = 0.0, when n > 0.
+  at[0] = 1;
+  *result.Element<CppTypeFor<CAT, KIND>>(at) = (n1 == 0) ? 1.0 : 0.0;
+  for (int j{2}; j <= n2 - n1 + 1; ++j) {
+    at[0] = j;
+    *result.Element<CppTypeFor<CAT, KIND>>(at) = 0.0;
+  }
+}
+
+template <TypeCategory CAT, int KIND>
+static inline void DoBesselYn(Descriptor &result, int32_t n1, int32_t n2,
+    CppTypeFor<CAT, KIND> x, CppTypeFor<CAT, KIND> bn1,
+    CppTypeFor<CAT, KIND> bn1_1, const char *sourceFile, int line) {
+  Terminator terminator{sourceFile, line};
+  AllocateBesselResult<CAT, KIND>(result, n1, n2, terminator, "BESSEL_YN");
+
+  // The standard requires that n1 and n2 be non-negative. However, some other
+  // compilers generate results even when n1 and/or n2 are negative. For now,
+  // we also do not enforce the non-negativity constraint.
+  if (n2 < n1) {
+    return;
+  }
+
+  SubscriptValue at[maxRank];
+  for (int j{0}; j < maxRank; ++j) {
+    at[j] = 0;
+  }
+
+  // if n2 >= n1, there will be at least one element in the result.
+  at[0] = 1;
+  *result.Element<CppTypeFor<CAT, KIND>>(at) = bn1;
+
+  if (n2 == n1) {
+    return;
+  }
+
+  at[0] = 2;
+  *result.Element<CppTypeFor<CAT, KIND>>(at) = bn1_1;
+
+  // Bessel functions of the second kind are stable for a forward recursion
+  // (see https://dlmf.nist.gov/10.74.iv and https://dlmf.nist.gov/10.6.E1).
+  //
+  //     Y(n+1, x) = (2.0 / x) * n * Y(n, x) - Y(n-1, x)
+  //
+  // which is equivalent to
+  //
+  //     Y(n, x) = (2.0 / x) * (n - 1) * Y(n-1, x) - Y(n-2, x)
+  //
+  CppTypeFor<CAT, KIND> bn_2 = bn1;
+  CppTypeFor<CAT, KIND> bn_1 = bn1_1;
+  CppTypeFor<CAT, KIND> twoOverX = 2.0 / x;
+  for (int n{n1 + 2}; n <= n2; ++n) {
+    auto bn = twoOverX * (n - 1) * bn_1 - bn_2;
+
+    at[0] = n - n1 + 1;
+    *result.Element<CppTypeFor<CAT, KIND>>(at) = bn;
+
+    bn_2 = bn_1;
+    bn_1 = bn;
+  }
+}
+
+template <TypeCategory CAT, int KIND>
+static inline void DoBesselYnX0(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile, int line) {
+  Terminator terminator{sourceFile, line};
+  AllocateBesselResult<CAT, KIND>(result, n1, n2, terminator, "BESSEL_YN");
+
+  // The standard requires that n1 and n2 be non-negative. However, some other
+  // compilers generate results even when n1 and/or n2 are negative. For now,
+  // we also do not enforce the non-negativity constraint.
+  if (n2 < n1) {
+    return;
+  }
+
+  SubscriptValue at[maxRank];
+  for (int j{0}; j < maxRank; ++j) {
+    at[j] = 0;
+  }
+
+  // Y(n, 0.0) = -Inf, when n >= 0
+  for (int j{1}; j <= n2 - n1 + 1; ++j) {
+    at[0] = j;
+    *result.Element<CppTypeFor<CAT, KIND>>(at) =
+        -std::numeric_limits<CppTypeFor<CAT, KIND>>::infinity();
+  }
+}
+
 extern "C" {
 
+// BESSEL_JN
+// TODO: REAL(2 & 3)
+void RTNAME(BesselJn_4)(Descriptor &result, int32_t n1, int32_t n2,
+    CppTypeFor<TypeCategory::Real, 4> x, CppTypeFor<TypeCategory::Real, 4> bn2,
+    CppTypeFor<TypeCategory::Real, 4> bn2_1, const char *sourceFile, int line) {
+  DoBesselJn<TypeCategory::Real, 4>(
+      result, n1, n2, x, bn2, bn2_1, sourceFile, line);
+}
+
+void RTNAME(BesselJn_8)(Descriptor &result, int32_t n1, int32_t n2,
+    CppTypeFor<TypeCategory::Real, 8> x, CppTypeFor<TypeCategory::Real, 8> bn2,
+    CppTypeFor<TypeCategory::Real, 8> bn2_1, const char *sourceFile, int line) {
+  DoBesselJn<TypeCategory::Real, 8>(
+      result, n1, n2, x, bn2, bn2_1, sourceFile, line);
+}
+
+#if LDBL_MANT_DIG == 64
+void RTNAME(BesselJn_10)(Descriptor &result, int32_t n1, int32_t n2,
+    CppTypeFor<TypeCategory::Real, 10> x,
+    CppTypeFor<TypeCategory::Real, 10> bn2,
+    CppTypeFor<TypeCategory::Real, 10> bn2_1, const char *sourceFile,
+    int line) {
+  DoBesselJn<TypeCategory::Real, 10>(
+      result, n1, n2, x, bn2, bn2_1, sourceFile, line);
+}
+#endif
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+void RTNAME(BesselJn_16)(Descriptor &result, int32_t n1, int32_t n2,
+    CppTypeFor<TypeCategory::Real, 16> x,
+    CppTypeFor<TypeCategory::Real, 16> bn2,
+    CppTypeFor<TypeCategory::Real, 16> bn2_1, const char *sourceFile,
+    int line) {
+  DoBesselJn<TypeCategory::Real, 16>(
+      result, n1, n2, x, bn2, bn2_1, sourceFile, line);
+}
+#endif
+
+// TODO: REAL(2 & 3)
+void RTNAME(BesselJnX0_4)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile, int line) {
+  DoBesselJnX0<TypeCategory::Real, 4>(result, n1, n2, sourceFile, line);
+}
+
+void RTNAME(BesselJnX0_8)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile, int line) {
+  DoBesselJnX0<TypeCategory::Real, 8>(result, n1, n2, sourceFile, line);
+}
+
+#if LDBL_MANT_DIG == 64
+void RTNAME(BesselJnX0_10)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile, int line) {
+  DoBesselJnX0<TypeCategory::Real, 10>(result, n1, n2, sourceFile, line);
+}
+#endif
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+void RTNAME(BesselJnX0_16)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile, int line) {
+  DoBesselJnX0<TypeCategory::Real, 16>(result, n1, n2, sourceFile, line);
+}
+#endif
+
+// BESSEL_YN
+// TODO: REAL(2 & 3)
+void RTNAME(BesselYn_4)(Descriptor &result, int32_t n1, int32_t n2,
+    CppTypeFor<TypeCategory::Real, 4> x, CppTypeFor<TypeCategory::Real, 4> bn1,
+    CppTypeFor<TypeCategory::Real, 4> bn1_1, const char *sourceFile, int line) {
+  DoBesselYn<TypeCategory::Real, 4>(
+      result, n1, n2, x, bn1, bn1_1, sourceFile, line);
+}
+
+void RTNAME(BesselYn_8)(Descriptor &result, int32_t n1, int32_t n2,
+    CppTypeFor<TypeCategory::Real, 8> x, CppTypeFor<TypeCategory::Real, 8> bn1,
+    CppTypeFor<TypeCategory::Real, 8> bn1_1, const char *sourceFile, int line) {
+  DoBesselYn<TypeCategory::Real, 8>(
+      result, n1, n2, x, bn1, bn1_1, sourceFile, line);
+}
+
+#if LDBL_MANT_DIG == 64
+void RTNAME(BesselYn_10)(Descriptor &result, int32_t n1, int32_t n2,
+    CppTypeFor<TypeCategory::Real, 10> x,
+    CppTypeFor<TypeCategory::Real, 10> bn1,
+    CppTypeFor<TypeCategory::Real, 10> bn1_1, const char *sourceFile,
+    int line) {
+  DoBesselYn<TypeCategory::Real, 10>(
+      result, n1, n2, x, bn1, bn1_1, sourceFile, line);
+}
+#endif
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+void RTNAME(BesselYn_16)(Descriptor &result, int32_t n1, int32_t n2,
+    CppTypeFor<TypeCategory::Real, 16> x,
+    CppTypeFor<TypeCategory::Real, 16> bn1,
+    CppTypeFor<TypeCategory::Real, 16> bn1_1, const char *sourceFile,
+    int line) {
+  DoBesselYn<TypeCategory::Real, 16>(
+      result, n1, n2, x, bn1, bn1_1, sourceFile, line);
+}
+#endif
+
+// TODO: REAL(2 & 3)
+void RTNAME(BesselYnX0_4)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile, int line) {
+  DoBesselYnX0<TypeCategory::Real, 4>(result, n1, n2, sourceFile, line);
+}
+
+void RTNAME(BesselYnX0_8)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile, int line) {
+  DoBesselYnX0<TypeCategory::Real, 8>(result, n1, n2, sourceFile, line);
+}
+
+#if LDBL_MANT_DIG == 64
+void RTNAME(BesselYnX0_10)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile, int line) {
+  DoBesselYnX0<TypeCategory::Real, 10>(result, n1, n2, sourceFile, line);
+}
+#endif
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+void RTNAME(BesselYnX0_16)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile, int line) {
+  DoBesselYnX0<TypeCategory::Real, 16>(result, n1, n2, sourceFile, line);
+}
+#endif
+
 // CSHIFT where rank of ARRAY argument > 1
 void RTNAME(Cshift)(Descriptor &result, const Descriptor &source,
     const Descriptor &shift, int dim, const char *sourceFile, int line) {

diff  --git a/flang/test/Lower/Intrinsics/bessel_jn.f90 b/flang/test/Lower/Intrinsics/bessel_jn.f90
index 047eac274c43..9e0698a5986f 100644
--- a/flang/test/Lower/Intrinsics/bessel_jn.f90
+++ b/flang/test/Lower/Intrinsics/bessel_jn.f90
@@ -5,20 +5,112 @@
 ! RUN: bbc -emit-fir %s -o - --math-runtime=precise | FileCheck --check-prefixes=ALL %s
 ! RUN: %flang_fc1 -emit-fir -mllvm -math-runtime=precise %s -o - | FileCheck --check-prefixes=ALL %s
 
+! ALL-LABEL: @_QPtest_real4
+! ALL-SAME: (%[[argx:.*]]: !fir.ref<f32>{{.*}}, %[[argn:.*]]: !fir.ref<i32>{{.*}}) -> f32
 function test_real4(x, n)
   real :: x, test_real4
   integer :: n
+
+  ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref<f32>
+  ! ALL-DAG: %[[n:.*]] = fir.load %[[argn]] : !fir.ref<i32>
+  ! ALL: fir.call @jnf(%[[n]], %[[x]]) {{.*}} : (i32, f32) -> f32
   test_real4 = bessel_jn(n, x)
 end function
 
-! ALL-LABEL: @_QPtest_real4
-! ALL: {{%[A-Za-z0-9._]+}} = fir.call @jnf({{%[A-Za-z0-9._]+}}, {{%[A-Za-z0-9._]+}}) {{.*}}: (i32, f32) -> f32
-
+! ALL-LABEL: @_QPtest_real8
+! ALL-SAME: (%[[argx:.*]]: !fir.ref<f64>{{.*}}, %[[argn:.*]]: !fir.ref<i32>{{.*}}) -> f64
 function test_real8(x, n)
   real(8) :: x, test_real8
   integer :: n
+
+  ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref<f64>
+  ! ALL-DAG: %[[n:.*]] = fir.load %[[argn]] : !fir.ref<i32>
+  ! ALL: fir.call @jn(%[[n]], %[[x]]) {{.*}} : (i32, f64) -> f64
   test_real8 = bessel_jn(n, x)
 end function
 
-! ALL-LABEL: @_QPtest_real8
-! ALL: {{%[A-Za-z0-9._]+}} = fir.call @jn({{%[A-Za-z0-9._]+}}, {{%[A-Za-z0-9._]+}}) {{.*}}: (i32, f64) -> f64
+! ALL-LABEL: @_QPtest_transformational_real4
+! ALL-SAME: %[[argx:.*]]: !fir.ref<f32>{{.*}}, %[[argn1:.*]]: !fir.ref<i32>{{.*}}, %[[argn2:.*]]: !fir.ref<i32>{{.*}}
+subroutine test_transformational_real4(x, n1, n2, r)
+  real(4) :: x
+  integer :: n1, n2
+  real(4) :: r(:)
+
+  ! ALL-DAG: %[[zero:.*]] = arith.constant 0{{.*}} : f32
+  ! ALL-DAG: %[[one:.*]] = arith.constant 1 : i32
+  ! ALL-DAG: %[[r:.*]] = fir.alloca !fir.box<!fir.heap<!fir.array<?xf32>>>
+  ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref<f32>
+  ! ALL-DAG: %[[n1:.*]] = fir.load %[[argn1]] : !fir.ref<i32>
+  ! ALL-DAG: %[[n2:.*]] = fir.load %[[argn2]] : !fir.ref<i32>
+  ! ALL-DAG: %[[xeq0:.*]] = arith.cmpf ueq, %[[x]], %[[zero]] : f32
+  ! ALL-DAG: %[[n1ltn2:.*]] = arith.cmpi slt, %[[n1]], %[[n2]] : i32
+  ! ALL-DAG: %[[n1eqn2:.*]] = arith.cmpi eq, %[[n1]], %[[n2]] : i32
+  ! ALL: fir.if %[[xeq0]] {
+  ! ALL: %[[resxeq0:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselJnX0_4(%[[resxeq0]], %[[n1]], %[[n2]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-NEXT: fir.if %[[n1ltn2]] {
+  ! ALL-DAG: %[[n2_1:.*]] = arith.subi %[[n2]], %[[one]] : i32
+  ! ALL-DAG: %[[bn2:.*]] = fir.call @jnf(%[[n2]], %[[x]]) {{.*}} : (i32, f32) -> f32
+  ! ALL-DAG: %[[bn2_1:.*]] = fir.call @jnf(%[[n2_1]], %[[x]]) {{.*}} : (i32, f32) -> f32
+  ! ALL-DAG: %[[resn1ltn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselJn_4(%[[resn1ltn2]], %[[n1]], %[[n2]], %[[x]], %[[bn2]], %[[bn2_1]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f32, f32, f32, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-NEXT: fir.if %[[n1eqn2]] {
+  ! ALL-DAG: %[[bn2:.*]] = fir.call @jnf(%[[n2]], %[[x]]) {{.*}} : (i32, f32) -> f32
+  ! ALL-DAG: %[[resn1eqn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselJn_4(%[[resn1eqn2]], %[[n1]], %[[n2]], %[[x]], %[[bn2]], %[[zero]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f32, f32, f32, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-DAG: %[[resn1gtn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselJn_4(%[[resn1gtn2]], %[[n1]], %[[n2]], %[[x]], %[[zero]], %[[zero]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f32, f32, f32, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: }
+  ! ALL-NEXT: }
+  ! ALL-NEXT: }
+  r = bessel_jn(n1, n2, x)
+  ! ALL: %[[box:.*]] = fir.load %[[r]] : !fir.ref<!fir.box<!fir.heap<!fir.array<?xf32>>>>
+  ! ALL: %[[addr:.*]] = fir.box_addr %[[box]] : (!fir.box<!fir.heap<!fir.array<?xf32>>>) -> !fir.heap<!fir.array<?xf32>>
+  ! ALL: fir.freemem %[[addr]] : !fir.heap<!fir.array<?xf32>>
+end subroutine test_transformational_real4
+
+! ALL-LABEL: @_QPtest_transformational_real8
+! ALL-SAME: %[[argx:.*]]: !fir.ref<f64>{{.*}}, %[[argn1:.*]]: !fir.ref<i32>{{.*}}, %[[argn2:.*]]: !fir.ref<i32>{{.*}}
+subroutine test_transformational_real8(x, n1, n2, r)
+  real(8) :: x
+  integer :: n1, n2
+  real(8) :: r(:)
+
+  ! ALL-DAG: %[[zero:.*]] = arith.constant 0{{.*}} : f64
+  ! ALL-DAG: %[[one:.*]] = arith.constant 1 : i32
+  ! ALL-DAG: %[[r:.*]] = fir.alloca !fir.box<!fir.heap<!fir.array<?xf64>>>
+  ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref<f64>
+  ! ALL-DAG: %[[n1:.*]] = fir.load %[[argn1]] : !fir.ref<i32>
+  ! ALL-DAG: %[[n2:.*]] = fir.load %[[argn2]] : !fir.ref<i32>
+  ! ALL-DAG: %[[xeq0:.*]] = arith.cmpf ueq, %[[x]], %[[zero]] : f64
+  ! ALL-DAG: %[[n1ltn2:.*]] = arith.cmpi slt, %[[n1]], %[[n2]] : i32
+  ! ALL-DAG: %[[n1eqn2:.*]] = arith.cmpi eq, %[[n1]], %[[n2]] : i32
+  ! ALL: fir.if %[[xeq0]] {
+  ! ALL: %[[resxeq0:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselJnX0_8(%[[resxeq0]], %[[n1]], %[[n2]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-NEXT: fir.if %[[n1ltn2]] {
+  ! ALL-DAG: %[[n2_1:.*]] = arith.subi %[[n2]], %[[one]] : i32
+  ! ALL-DAG: %[[bn2:.*]] = fir.call @jn(%[[n2]], %[[x]]) {{.*}} : (i32, f64) -> f64
+  ! ALL-DAG: %[[bn2_1:.*]] = fir.call @jn(%[[n2_1]], %[[x]]) {{.*}} : (i32, f64) -> f64
+  ! ALL-DAG: %[[resn1ltn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselJn_8(%[[resn1ltn2]], %[[n1]], %[[n2]], %[[x]], %[[bn2]], %[[bn2_1]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f64, f64, f64, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-NEXT: fir.if %[[n1eqn2]] {
+  ! ALL-DAG: %[[bn2:.*]] = fir.call @jn(%[[n2]], %[[x]]) {{.*}} : (i32, f64) -> f64
+  ! ALL-DAG: %[[resn1eqn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselJn_8(%[[resn1eqn2]], %[[n1]], %[[n2]], %[[x]], %[[bn2]], %[[zero]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f64, f64, f64, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-DAG: %[[resn1gtn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselJn_8(%[[resn1gtn2]], %[[n1]], %[[n2]], %[[x]], %[[zero]], %[[zero]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f64, f64, f64, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: }
+  ! ALL-NEXT: }
+  ! ALL-NEXT: }
+  r = bessel_jn(n1, n2, x)
+  ! ALL: %[[box:.*]] = fir.load %[[r]] : !fir.ref<!fir.box<!fir.heap<!fir.array<?xf64>>>>
+  ! ALL: %[[addr:.*]] = fir.box_addr %[[box]] : (!fir.box<!fir.heap<!fir.array<?xf64>>>) -> !fir.heap<!fir.array<?xf64>>
+  ! ALL: fir.freemem %[[addr]] : !fir.heap<!fir.array<?xf64>>
+end subroutine test_transformational_real8

diff  --git a/flang/test/Lower/Intrinsics/bessel_yn.f90 b/flang/test/Lower/Intrinsics/bessel_yn.f90
index a87642a2d3ca..cf1541dac0ec 100644
--- a/flang/test/Lower/Intrinsics/bessel_yn.f90
+++ b/flang/test/Lower/Intrinsics/bessel_yn.f90
@@ -5,20 +5,112 @@
 ! RUN: bbc -emit-fir %s -o - --math-runtime=precise | FileCheck --check-prefixes=ALL %s
 ! RUN: %flang_fc1 -emit-fir -mllvm -math-runtime=precise %s -o - | FileCheck --check-prefixes=ALL %s
 
+! ALL-LABEL: @_QPtest_real4
+! ALL-SAME: (%[[argx:.*]]: !fir.ref<f32>{{.*}}, %[[argn:.*]]: !fir.ref<i32>{{.*}}) -> f32
 function test_real4(x, n)
   real :: x, test_real4
   integer :: n
+
+  ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref<f32>
+  ! ALL-DAG: %[[n:.*]] = fir.load %[[argn]] : !fir.ref<i32>
+  ! ALL: fir.call @ynf(%[[n]], %[[x]]) {{.*}} : (i32, f32) -> f32
   test_real4 = bessel_yn(n, x)
 end function
 
-! ALL-LABEL: @_QPtest_real4
-! ALL: {{%[A-Za-z0-9._]+}} = fir.call @ynf({{%[A-Za-z0-9._]+}}, {{%[A-Za-z0-9._]+}}) {{.*}}: (i32, f32) -> f32
-
+! ALL-LABEL: @_QPtest_real8
+! ALL-SAME: (%[[argx:.*]]: !fir.ref<f64>{{.*}}, %[[argn:.*]]: !fir.ref<i32>{{.*}}) -> f64
 function test_real8(x, n)
   real(8) :: x, test_real8
   integer :: n
+
+  ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref<f64>
+  ! ALL-DAG: %[[n:.*]] = fir.load %[[argn]] : !fir.ref<i32>
+  ! ALL: fir.call @yn(%[[n]], %[[x]]) {{.*}} : (i32, f64) -> f64
   test_real8 = bessel_yn(n, x)
 end function
 
-! ALL-LABEL: @_QPtest_real8
-! ALL: {{%[A-Za-z0-9._]+}} = fir.call @yn({{%[A-Za-z0-9._]+}}, {{%[A-Za-z0-9._]+}}) {{.*}}: (i32, f64) -> f64
+! ALL-LABEL: @_QPtest_transformational_real4
+! ALL-SAME: %[[argx:.*]]: !fir.ref<f32>{{.*}}, %[[argn1:.*]]: !fir.ref<i32>{{.*}}, %[[argn2:.*]]: !fir.ref<i32>{{.*}}
+subroutine test_transformational_real4(x, n1, n2, r)
+  real(4) :: x
+  integer :: n1, n2
+  real(4) :: r(:)
+
+  ! ALL-DAG: %[[zero:.*]] = arith.constant 0{{.*}} : f32
+  ! ALL-DAG: %[[one:.*]] = arith.constant 1 : i32
+  ! ALL-DAG: %[[r:.*]] = fir.alloca !fir.box<!fir.heap<!fir.array<?xf32>>>
+  ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref<f32>
+  ! ALL-DAG: %[[n1:.*]] = fir.load %[[argn1]] : !fir.ref<i32>
+  ! ALL-DAG: %[[n2:.*]] = fir.load %[[argn2]] : !fir.ref<i32>
+  ! ALL-DAG: %[[xeq0:.*]] = arith.cmpf ueq, %[[x]], %[[zero]] : f32
+  ! ALL-DAG: %[[n1ltn2:.*]] = arith.cmpi slt, %[[n1]], %[[n2]] : i32
+  ! ALL-DAG: %[[n1eqn2:.*]] = arith.cmpi eq, %[[n1]], %[[n2]] : i32
+  ! ALL: fir.if %[[xeq0]] {
+  ! ALL: %[[resxeq0:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselYnX0_4(%[[resxeq0]], %[[n1]], %[[n2]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-NEXT: fir.if %[[n1ltn2]] {
+  ! ALL-DAG: %[[n1_1:.*]] = arith.addi %[[n1]], %[[one]] : i32
+  ! ALL-DAG: %[[bn1:.*]] = fir.call @ynf(%[[n1]], %[[x]]) {{.*}} : (i32, f32) -> f32
+  ! ALL-DAG: %[[bn1_1:.*]] = fir.call @ynf(%[[n1_1]], %[[x]]) {{.*}} : (i32, f32) -> f32
+  ! ALL-DAG: %[[resn1ltn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselYn_4(%[[resn1ltn2]], %[[n1]], %[[n2]], %[[x]], %[[bn1]], %[[bn1_1]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f32, f32, f32, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-NEXT: fir.if %[[n1eqn2]] {
+  ! ALL-DAG: %[[bn1:.*]] = fir.call @ynf(%[[n1]], %[[x]]) {{.*}} : (i32, f32) -> f32
+  ! ALL-DAG: %[[resn1eqn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselYn_4(%[[resn1eqn2]], %[[n1]], %[[n2]], %[[x]], %[[bn1]], %[[zero]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f32, f32, f32, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-DAG: %[[resn1gtn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselYn_4(%[[resn1gtn2]], %[[n1]], %[[n2]], %[[x]], %[[zero]], %[[zero]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f32, f32, f32, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: }
+  ! ALL-NEXT: }
+  ! ALL-NEXT: }
+  r = bessel_yn(n1, n2, x)
+  ! ALL: %[[box:.*]] = fir.load %[[r]] : !fir.ref<!fir.box<!fir.heap<!fir.array<?xf32>>>>
+  ! ALL: %[[addr:.*]] = fir.box_addr %[[box]] : (!fir.box<!fir.heap<!fir.array<?xf32>>>) -> !fir.heap<!fir.array<?xf32>>
+  ! ALL: fir.freemem %[[addr]] : !fir.heap<!fir.array<?xf32>>
+end subroutine test_transformational_real4
+
+! ALL-LABEL: @_QPtest_transformational_real8
+! ALL-SAME: %[[argx:.*]]: !fir.ref<f64>{{.*}}, %[[argn1:.*]]: !fir.ref<i32>{{.*}}, %[[argn2:.*]]: !fir.ref<i32>{{.*}}
+subroutine test_transformational_real8(x, n1, n2, r)
+  real(8) :: x
+  integer :: n1, n2
+  real(8) :: r(:)
+
+  ! ALL-DAG: %[[zero:.*]] = arith.constant 0{{.*}} : f64
+  ! ALL-DAG: %[[one:.*]] = arith.constant 1 : i32
+  ! ALL-DAG: %[[r:.*]] = fir.alloca !fir.box<!fir.heap<!fir.array<?xf64>>>
+  ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref<f64>
+  ! ALL-DAG: %[[n1:.*]] = fir.load %[[argn1]] : !fir.ref<i32>
+  ! ALL-DAG: %[[n2:.*]] = fir.load %[[argn2]] : !fir.ref<i32>
+  ! ALL-DAG: %[[xeq0:.*]] = arith.cmpf ueq, %[[x]], %[[zero]] : f64
+  ! ALL-DAG: %[[n1ltn2:.*]] = arith.cmpi slt, %[[n1]], %[[n2]] : i32
+  ! ALL-DAG: %[[n1eqn2:.*]] = arith.cmpi eq, %[[n1]], %[[n2]] : i32
+  ! ALL: fir.if %[[xeq0]] {
+  ! ALL: %[[resxeq0:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselYnX0_8(%[[resxeq0]], %[[n1]], %[[n2]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-NEXT: fir.if %[[n1ltn2]] {
+  ! ALL-DAG: %[[n1_1:.*]] = arith.addi %[[n1]], %[[one]] : i32
+  ! ALL-DAG: %[[bn1:.*]] = fir.call @yn(%[[n1]], %[[x]]) {{.*}} : (i32, f64) -> f64
+  ! ALL-DAG: %[[bn1_1:.*]] = fir.call @yn(%[[n1_1]], %[[x]]) {{.*}} : (i32, f64) -> f64
+  ! ALL-DAG: %[[resn1ltn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselYn_8(%[[resn1ltn2]], %[[n1]], %[[n2]], %[[x]], %[[bn1]], %[[bn1_1]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f64, f64, f64, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-NEXT: fir.if %[[n1eqn2]] {
+  ! ALL-DAG: %[[bn1:.*]] = fir.call @yn(%[[n1]], %[[x]]) {{.*}} : (i32, f64) -> f64
+  ! ALL-DAG: %[[resn1eqn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselYn_8(%[[resn1eqn2]], %[[n1]], %[[n2]], %[[x]], %[[bn1]], %[[zero]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f64, f64, f64, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-DAG: %[[resn1gtn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselYn_8(%[[resn1gtn2]], %[[n1]], %[[n2]], %[[x]], %[[zero]], %[[zero]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f64, f64, f64, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: }
+  ! ALL-NEXT: }
+  ! ALL-NEXT: }
+  r = bessel_yn(n1, n2, x)
+  ! ALL: %[[box:.*]] = fir.load %[[r]] : !fir.ref<!fir.box<!fir.heap<!fir.array<?xf64>>>>
+  ! ALL: %[[addr:.*]] = fir.box_addr %[[box]] : (!fir.box<!fir.heap<!fir.array<?xf64>>>) -> !fir.heap<!fir.array<?xf64>>
+  ! ALL: fir.freemem %[[addr]] : !fir.heap<!fir.array<?xf64>>
+end subroutine test_transformational_real8

diff  --git a/flang/unittests/Optimizer/Builder/Runtime/TransformationalTest.cpp b/flang/unittests/Optimizer/Builder/Runtime/TransformationalTest.cpp
index 37e13ccb06b5..d5884ae3febb 100644
--- a/flang/unittests/Optimizer/Builder/Runtime/TransformationalTest.cpp
+++ b/flang/unittests/Optimizer/Builder/Runtime/TransformationalTest.cpp
@@ -10,6 +10,92 @@
 #include "RuntimeCallTestBase.h"
 #include "gtest/gtest.h"
 
+void testGenBesselJn(
+    fir::FirOpBuilder &builder, mlir::Type realTy, llvm::StringRef fctName) {
+  mlir::Location loc = builder.getUnknownLoc();
+  mlir::Type i32Ty = builder.getIntegerType(32);
+  mlir::Type seqTy =
+      fir::SequenceType::get(fir::SequenceType::Shape(1, 10), realTy);
+  mlir::Value result = builder.create<fir::UndefOp>(loc, seqTy);
+  mlir::Value n1 = builder.create<fir::UndefOp>(loc, i32Ty);
+  mlir::Value n2 = builder.create<fir::UndefOp>(loc, i32Ty);
+  mlir::Value x = builder.create<fir::UndefOp>(loc, realTy);
+  mlir::Value bn1 = builder.create<fir::UndefOp>(loc, realTy);
+  mlir::Value bn2 = builder.create<fir::UndefOp>(loc, realTy);
+  fir::runtime::genBesselJn(builder, loc, result, n1, n2, x, bn1, bn2);
+  checkCallOpFromResultBox(result, fctName, 6);
+}
+
+TEST_F(RuntimeCallTest, genBesselJnTest) {
+  testGenBesselJn(*firBuilder, f32Ty, "_FortranABesselJn_4");
+  testGenBesselJn(*firBuilder, f64Ty, "_FortranABesselJn_8");
+  testGenBesselJn(*firBuilder, f80Ty, "_FortranABesselJn_10");
+  testGenBesselJn(*firBuilder, f128Ty, "_FortranABesselJn_16");
+}
+
+void testGenBesselJnX0(
+    fir::FirOpBuilder &builder, mlir::Type realTy, llvm::StringRef fctName) {
+  mlir::Location loc = builder.getUnknownLoc();
+  mlir::Type i32Ty = builder.getIntegerType(32);
+  mlir::Type seqTy =
+      fir::SequenceType::get(fir::SequenceType::Shape(1, 10), realTy);
+  mlir::Value result = builder.create<fir::UndefOp>(loc, seqTy);
+  mlir::Value n1 = builder.create<fir::UndefOp>(loc, i32Ty);
+  mlir::Value n2 = builder.create<fir::UndefOp>(loc, i32Ty);
+  fir::runtime::genBesselJnX0(builder, loc, realTy, result, n1, n2);
+  checkCallOpFromResultBox(result, fctName, 3);
+}
+
+TEST_F(RuntimeCallTest, genBesselJnX0Test) {
+  testGenBesselJnX0(*firBuilder, f32Ty, "_FortranABesselJnX0_4");
+  testGenBesselJnX0(*firBuilder, f64Ty, "_FortranABesselJnX0_8");
+  testGenBesselJnX0(*firBuilder, f80Ty, "_FortranABesselJnX0_10");
+  testGenBesselJnX0(*firBuilder, f128Ty, "_FortranABesselJnX0_16");
+}
+
+void testGenBesselYn(
+    fir::FirOpBuilder &builder, mlir::Type realTy, llvm::StringRef fctName) {
+  mlir::Location loc = builder.getUnknownLoc();
+  mlir::Type i32Ty = builder.getIntegerType(32);
+  mlir::Type seqTy =
+      fir::SequenceType::get(fir::SequenceType::Shape(1, 10), realTy);
+  mlir::Value result = builder.create<fir::UndefOp>(loc, seqTy);
+  mlir::Value n1 = builder.create<fir::UndefOp>(loc, i32Ty);
+  mlir::Value n2 = builder.create<fir::UndefOp>(loc, i32Ty);
+  mlir::Value x = builder.create<fir::UndefOp>(loc, realTy);
+  mlir::Value bn1 = builder.create<fir::UndefOp>(loc, realTy);
+  mlir::Value bn2 = builder.create<fir::UndefOp>(loc, realTy);
+  fir::runtime::genBesselYn(builder, loc, result, n1, n2, x, bn1, bn2);
+  checkCallOpFromResultBox(result, fctName, 6);
+}
+
+TEST_F(RuntimeCallTest, genBesselYnTest) {
+  testGenBesselYn(*firBuilder, f32Ty, "_FortranABesselYn_4");
+  testGenBesselYn(*firBuilder, f64Ty, "_FortranABesselYn_8");
+  testGenBesselYn(*firBuilder, f80Ty, "_FortranABesselYn_10");
+  testGenBesselYn(*firBuilder, f128Ty, "_FortranABesselYn_16");
+}
+
+void testGenBesselYnX0(
+    fir::FirOpBuilder &builder, mlir::Type realTy, llvm::StringRef fctName) {
+  mlir::Location loc = builder.getUnknownLoc();
+  mlir::Type i32Ty = builder.getIntegerType(32);
+  mlir::Type seqTy =
+      fir::SequenceType::get(fir::SequenceType::Shape(1, 10), realTy);
+  mlir::Value result = builder.create<fir::UndefOp>(loc, seqTy);
+  mlir::Value n1 = builder.create<fir::UndefOp>(loc, i32Ty);
+  mlir::Value n2 = builder.create<fir::UndefOp>(loc, i32Ty);
+  fir::runtime::genBesselYnX0(builder, loc, realTy, result, n1, n2);
+  checkCallOpFromResultBox(result, fctName, 3);
+}
+
+TEST_F(RuntimeCallTest, genBesselYnX0Test) {
+  testGenBesselYnX0(*firBuilder, f32Ty, "_FortranABesselYnX0_4");
+  testGenBesselYnX0(*firBuilder, f64Ty, "_FortranABesselYnX0_8");
+  testGenBesselYnX0(*firBuilder, f80Ty, "_FortranABesselYnX0_10");
+  testGenBesselYnX0(*firBuilder, f128Ty, "_FortranABesselYnX0_16");
+}
+
 TEST_F(RuntimeCallTest, genCshiftTest) {
   auto loc = firBuilder->getUnknownLoc();
   mlir::Type seqTy =

diff  --git a/flang/unittests/Runtime/Transformational.cpp b/flang/unittests/Runtime/Transformational.cpp
index a4672369f782..70ab42416564 100644
--- a/flang/unittests/Runtime/Transformational.cpp
+++ b/flang/unittests/Runtime/Transformational.cpp
@@ -10,10 +10,215 @@
 #include "gtest/gtest.h"
 #include "tools.h"
 #include "flang/Runtime/type-code.h"
+#include <vector>
 
 using namespace Fortran::runtime;
 using Fortran::common::TypeCategory;
 
+template <int KIND>
+using BesselFuncType = std::function<void(Descriptor &, int32_t, int32_t,
+    CppTypeFor<TypeCategory::Real, KIND>, CppTypeFor<TypeCategory::Real, KIND>,
+    CppTypeFor<TypeCategory::Real, KIND>, const char *, int)>;
+
+template <int KIND>
+using BesselX0FuncType =
+    std::function<void(Descriptor &, int32_t, int32_t, const char *, int)>;
+
+template <int KIND>
+constexpr CppTypeFor<TypeCategory::Real, KIND>
+    besselEpsilon = CppTypeFor<TypeCategory::Real, KIND>(1e-4);
+
+template <int KIND>
+static void testBesselJn(BesselFuncType<KIND> rtFunc, int32_t n1, int32_t n2,
+    CppTypeFor<TypeCategory::Real, KIND> x,
+    const std::vector<CppTypeFor<TypeCategory::Real, KIND>> &expected) {
+  StaticDescriptor<1> desc;
+  Descriptor &result{desc.descriptor()};
+  unsigned len = expected.size();
+
+  CppTypeFor<TypeCategory::Real, KIND> anc0 = len > 0 ? expected[len - 1] : 0.0;
+  CppTypeFor<TypeCategory::Real, KIND> anc1 = len > 1 ? expected[len - 2] : 0.0;
+
+  rtFunc(result, n1, n2, x, anc0, anc1, __FILE__, __LINE__);
+
+  EXPECT_EQ(result.type(), (TypeCode{TypeCategory::Real, KIND}));
+  EXPECT_EQ(result.rank(), 1);
+  EXPECT_EQ(
+      result.ElementBytes(), sizeof(CppTypeFor<TypeCategory::Real, KIND>));
+  EXPECT_EQ(result.GetDimension(0).LowerBound(), 1);
+  EXPECT_EQ(result.GetDimension(0).Extent(), len);
+
+  for (size_t j{0}; j < len; ++j) {
+    EXPECT_NEAR(
+        (*result.ZeroBasedIndexedElement<CppTypeFor<TypeCategory::Real, KIND>>(
+            j)),
+        expected[j], besselEpsilon<KIND>);
+  }
+}
+
+template <int KIND>
+static void testBesselJnX0(
+    BesselX0FuncType<KIND> rtFunc, int32_t n1, int32_t n2) {
+  StaticDescriptor<1> desc;
+  Descriptor &result{desc.descriptor()};
+
+  rtFunc(result, n1, n2, __FILE__, __LINE__);
+
+  EXPECT_EQ(result.type(), (TypeCode{TypeCategory::Real, KIND}));
+  EXPECT_EQ(result.rank(), 1);
+  EXPECT_EQ(result.GetDimension(0).LowerBound(), 1);
+  EXPECT_EQ(result.GetDimension(0).Extent(), n2 >= n1 ? n2 - n1 + 1 : 0);
+
+  if (n2 < n1) {
+    return;
+  }
+
+  EXPECT_NEAR(
+      (*result.ZeroBasedIndexedElement<CppTypeFor<TypeCategory::Real, KIND>>(
+          0)),
+      (n1 == 0) ? 1.0 : 0.0, 1e-5);
+
+  for (int j{1}; j < (n2 - n1 + 1); ++j) {
+    EXPECT_NEAR(
+        (*result.ZeroBasedIndexedElement<CppTypeFor<TypeCategory::Real, KIND>>(
+            j)),
+        0.0, besselEpsilon<KIND>);
+  }
+}
+
+template <int KIND> static void testBesselJn(BesselFuncType<KIND> rtFunc) {
+  testBesselJn<KIND>(rtFunc, 1, 0, 1.0, {});
+  testBesselJn<KIND>(rtFunc, 0, 0, 1.0, {0.765197694});
+  testBesselJn<KIND>(rtFunc, 0, 1, 1.0, {0.765197694, 0.440050572});
+  testBesselJn<KIND>(
+      rtFunc, 0, 2, 1.0, {0.765197694, 0.440050572, 0.114903487});
+  testBesselJn<KIND>(rtFunc, 1, 5, 5.0,
+      {-0.327579111, 0.046565145, 0.364831239, 0.391232371, 0.261140555});
+}
+
+template <int KIND> static void testBesselJnX0(BesselX0FuncType<KIND> rtFunc) {
+  testBesselJnX0<KIND>(rtFunc, 1, 0);
+  testBesselJnX0<KIND>(rtFunc, 0, 0);
+  testBesselJnX0<KIND>(rtFunc, 1, 1);
+  testBesselJnX0<KIND>(rtFunc, 0, 3);
+  testBesselJnX0<KIND>(rtFunc, 1, 4);
+}
+
+static void testBesselJn() {
+  testBesselJn<4>(RTNAME(BesselJn_4));
+  testBesselJn<8>(RTNAME(BesselJn_8));
+#if LDBL_MANT_DIG == 64
+  testBesselJn<10>(RTNAME(BesselJn_10));
+#endif
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+  testBesselJn<16>(RTNAME(BesselJn_16));
+#endif
+
+  testBesselJnX0<4>(RTNAME(BesselJnX0_4));
+  testBesselJnX0<8>(RTNAME(BesselJnX0_8));
+#if LDBL_MANT_DIG == 64
+  testBesselJnX0<10>(RTNAME(BesselJnX0_10));
+#endif
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+  testBesselJnX0<16>(RTNAME(BesselJnX0_16));
+#endif
+}
+
+TEST(Transformational, BesselJn) { testBesselJn(); }
+
+template <int KIND>
+static void testBesselYn(BesselFuncType<KIND> rtFunc, int32_t n1, int32_t n2,
+    CppTypeFor<TypeCategory::Real, KIND> x,
+    const std::vector<CppTypeFor<TypeCategory::Real, KIND>> &expected) {
+  StaticDescriptor<1> desc;
+  Descriptor &result{desc.descriptor()};
+  unsigned len = expected.size();
+
+  CppTypeFor<TypeCategory::Real, KIND> anc0 = len > 0 ? expected[0] : 0.0;
+  CppTypeFor<TypeCategory::Real, KIND> anc1 = len > 1 ? expected[1] : 0.0;
+
+  rtFunc(result, n1, n2, x, anc0, anc1, __FILE__, __LINE__);
+
+  EXPECT_EQ(result.type(), (TypeCode{TypeCategory::Real, KIND}));
+  EXPECT_EQ(result.rank(), 1);
+  EXPECT_EQ(
+      result.ElementBytes(), sizeof(CppTypeFor<TypeCategory::Real, KIND>));
+  EXPECT_EQ(result.GetDimension(0).LowerBound(), 1);
+  EXPECT_EQ(result.GetDimension(0).Extent(), len);
+
+  for (size_t j{0}; j < len; ++j) {
+    EXPECT_NEAR(
+        (*result.ZeroBasedIndexedElement<CppTypeFor<TypeCategory::Real, KIND>>(
+            j)),
+        expected[j], besselEpsilon<KIND>);
+  }
+}
+
+template <int KIND>
+static void testBesselYnX0(
+    BesselX0FuncType<KIND> rtFunc, int32_t n1, int32_t n2) {
+  StaticDescriptor<1> desc;
+  Descriptor &result{desc.descriptor()};
+
+  rtFunc(result, n1, n2, __FILE__, __LINE__);
+
+  EXPECT_EQ(result.type(), (TypeCode{TypeCategory::Real, KIND}));
+  EXPECT_EQ(result.rank(), 1);
+  EXPECT_EQ(result.GetDimension(0).LowerBound(), 1);
+  EXPECT_EQ(result.GetDimension(0).Extent(), n2 >= n1 ? n2 - n1 + 1 : 0);
+
+  if (n2 < n1) {
+    return;
+  }
+
+  for (int j{0}; j < (n2 - n1 + 1); ++j) {
+    EXPECT_EQ(
+        (*result.ZeroBasedIndexedElement<CppTypeFor<TypeCategory::Real, KIND>>(
+            j)),
+        (-std::numeric_limits<
+            CppTypeFor<TypeCategory::Real, KIND>>::infinity()));
+  }
+}
+
+template <int KIND> static void testBesselYn(BesselFuncType<KIND> rtFunc) {
+  testBesselYn<KIND>(rtFunc, 1, 0, 1.0, {});
+  testBesselYn<KIND>(rtFunc, 0, 0, 1.0, {0.08825695});
+  testBesselYn<KIND>(rtFunc, 0, 1, 1.0, {0.08825695, -0.7812128});
+  testBesselYn<KIND>(rtFunc, 0, 2, 1.0, {0.0882569555, -0.7812128, -1.6506826});
+  testBesselYn<KIND>(rtFunc, 1, 5, 1.0,
+      {-0.7812128, -1.6506826, -5.8215175, -33.278423, -260.40588});
+}
+
+template <int KIND> static void testBesselYnX0(BesselX0FuncType<KIND> rtFunc) {
+  testBesselYnX0<KIND>(rtFunc, 1, 0);
+  testBesselYnX0<KIND>(rtFunc, 0, 0);
+  testBesselYnX0<KIND>(rtFunc, 1, 1);
+  testBesselYnX0<KIND>(rtFunc, 0, 3);
+  testBesselYnX0<KIND>(rtFunc, 1, 4);
+}
+
+static void testBesselYn() {
+  testBesselYn<4>(RTNAME(BesselYn_4));
+  testBesselYn<8>(RTNAME(BesselYn_8));
+#if LDBL_MANT_DIG == 64
+  testBesselYn<10>(RTNAME(BesselYn_10));
+#endif
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+  testBesselYn<16>(RTNAME(BesselYn_16));
+#endif
+
+  testBesselYnX0<4>(RTNAME(BesselYnX0_4));
+  testBesselYnX0<8>(RTNAME(BesselYnX0_8));
+#if LDBL_MANT_DIG == 64
+  testBesselYnX0<10>(RTNAME(BesselYnX0_10));
+#endif
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+  testBesselYnX0<16>(RTNAME(BesselYnX0_16));
+#endif
+}
+
+TEST(Transformational, BesselYn) { testBesselYn(); }
+
 TEST(Transformational, Shifts) {
   // ARRAY  1 3 5
   //        2 4 6


        


More information about the flang-commits mailing list