[flang] Lowering and runtime support for F08 transformational intrinsics: BESSEL_JN...
authorTarun Prabhu <tarun@lanl.gov>
Mon, 19 Dec 2022 14:51:59 +0000 (07:51 -0700)
committerTarun Prabhu <tarun@lanl.gov>
Mon, 19 Dec 2022 14:59:38 +0000 (07:59 -0700)
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.

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

index 1657dff..f084e2e 100644 (file)
@@ -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);
index 21a4418..8101c73 100644 (file)
@@ -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);
index 65bb310..cdccd45 100644 (file)
@@ -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
index 0b8e145..c51f386 100644 (file)
 
 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,
index 9a01681..00bf12f 100644 (file)
@@ -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) {
index 047eac2..9e0698a 100644 (file)
 ! 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
index a87642a..cf1541d 100644 (file)
 ! 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
index 37e13cc..d5884ae 100644 (file)
 #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 =
index a467236..70ab424 100644 (file)
 #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