[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.
This commit is contained in:
Tarun Prabhu 2022-12-19 07:51:59 -07:00
parent 05b060b0b0
commit bef2bb34bf
9 changed files with 1358 additions and 11 deletions

View 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);

View 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);

View 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

View File

@ -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,

View 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) {

View File

@ -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

View File

@ -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

View File

@ -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 =

View File

@ -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