Skip to content

Commit

Permalink
Fix implementation for cyl_bessel_i0
Browse files Browse the repository at this point in the history
  • Loading branch information
@masterleinad
masterleinad committed Oct 4, 2023
1 parent 4ce289b commit 6ff0deb
Show file tree
Hide file tree
Showing 2 changed files with 36 additions and 50 deletions.
81 changes: 32 additions & 49 deletions core/src/Kokkos_MathematicalSpecialFunctions.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -846,69 +846,52 @@ KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_y1(const CmplxType& z,
//! for a complex argument
template <class CmplxType, class RealType, class IntType>
KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_i0(const CmplxType& z,
const RealType& joint_val = 25,
const IntType& bw_start = 70) {
const RealType& joint_val = 18,
const IntType& n_terms = 50) {
// This function is converted and modified from the corresponding Fortran
// programs CIKNB and CIK01 in S. Zhang & J. Jin "Computation of Special
// programs CIK01 in S. Zhang & J. Jin "Computation of Special
// Functions" (Wiley, 1996).
// Input : z --- Complex argument
// joint_val --- Joint point of abs(z) separating small and large
// argument regions
// bw_start --- Starting point for backward recurrence
// n_terms --- Numbers of terms used in the power series
// Output: cbi0 --- I0(z)
using Kokkos::numbers::pi_v;

CmplxType cbi0;
constexpr auto pi = pi_v<RealType>;
const RealType a[12] = {0.125,
7.03125e-2,
7.32421875e-2,
1.1215209960938e-1,
2.2710800170898e-1,
5.7250142097473e-1,
1.7277275025845e0,
6.0740420012735e0,
2.4380529699556e1,
1.1001714026925e2,
5.5133589612202e2,
3.0380905109224e3};

CmplxType cbi0(1.0, 0.0);
RealType a0 = Kokkos::abs(z);
CmplxType z1 = z;

if (a0 < 1e-100) { // Treat z=0 as a special case
cbi0 = CmplxType(1.0, 0.0);
} else {
if (a0 > 1e-100) {
if (z.real() < 0.0) z1 = -z;
if (a0 <= joint_val) { // Using backward recurrence for |z|<=joint_val
// (default:25)
CmplxType cbs = CmplxType(0.0, 0.0);
// CmplxType csk0 = CmplxType(0.0,0.0);
CmplxType cf0 = CmplxType(0.0, 0.0);
CmplxType cf1 = CmplxType(1e-100, 0.0);
CmplxType cf, cs0;
for (int k = bw_start; k >= 0; k--) { // Backward recurrence (default:
// 70)
cf = 2.0 * (k + 1.0) * cf1 / z1 + cf0;
if (k == 0) cbi0 = cf;
// if ((k == 2*(k/2)) && (k != 0)) {
// csk0 = csk0+4.0*cf/static_cast<RealType>(k);
//}
cbs = cbs + 2.0 * cf;
cf0 = cf1;
cf1 = cf;
if (a0 <= joint_val) {
// Using power series definition for |z|<=joint_val (default:18)
CmplxType cr = CmplxType(1.0e+00, 0.0e+00);
CmplxType z2 = z * z;
for (int k = 1; k < n_terms; ++k) {
cr = RealType(.25) * cr * z2 / CmplxType(k * k);
cbi0 += cr;
if (Kokkos::abs(cr / cbi0) < RealType(1.e-15)) continue;
}
cs0 = Kokkos::exp(z1) / (cbs - cf);
cbi0 = cbi0 * cs0;
} else { // Using asymptotic expansion (6.2.1) for |z|>joint_val
// (default:25)
CmplxType ca = Kokkos::exp(z1) / Kokkos::sqrt(2.0 * pi * z1);
cbi0 = CmplxType(1.0, 0.0);
CmplxType zr = 1.0 / z1;
} else {
// Using asymptotic expansion (6.2.1) for |z|>joint_val (default:18)
const RealType a[12] = {0.125,
7.03125e-2,
7.32421875e-2,
1.1215209960938e-1,
2.2710800170898e-1,
5.7250142097473e-1,
1.7277275025845e0,
6.0740420012735e0,
2.4380529699556e1,
1.1001714026925e2,
5.5133589612202e2,
3.0380905109224e3};

for (int k = 1; k <= 12; k++) {
cbi0 = cbi0 + a[k - 1] * Kokkos::pow(zr, 1.0 * k);
cbi0 += a[k - 1] * Kokkos::pow(z1, -k);
}
cbi0 = ca * cbi0;
cbi0 *= Kokkos::exp(z1) /
Kokkos::sqrt(2.0 * Kokkos::numbers::pi_v<RealType> * z1);
}
}
return cbi0;
Expand Down
5 changes: 4 additions & 1 deletion core/unit_test/TestMathematicalSpecialFunctions.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1058,7 +1058,7 @@ struct TestComplexBesselI0K0Function {
void testit() {
using Kokkos::Experimental::infinity;

int N = 25;
int N = 26;
d_z = ViewType("d_z", N);
d_cbi0 = ViewType("d_cbi0", N);
d_cbk0 = ViewType("d_cbk0", N);
Expand Down Expand Up @@ -1094,6 +1094,7 @@ struct TestComplexBesselI0K0Function {
h_z(22) = Kokkos::complex<double>(-28.0, 0.0);
h_z(23) = Kokkos::complex<double>(60.0, 0.0);
h_z(24) = Kokkos::complex<double>(-60.0, 0.0);
h_z(25) = Kokkos::complex<double>(7.998015e-5, 0.0);

Kokkos::deep_copy(d_z, h_z);

Expand Down Expand Up @@ -1152,6 +1153,7 @@ struct TestComplexBesselI0K0Function {
h_ref_cbi0(22) = Kokkos::complex<double>(1.095346047317573e+11, 0);
h_ref_cbi0(23) = Kokkos::complex<double>(5.894077055609803e+24, 0);
h_ref_cbi0(24) = Kokkos::complex<double>(5.894077055609803e+24, 0);
h_ref_cbi0(25) = Kokkos::complex<double>(1.0000000015992061009, 0);

h_ref_cbk0(0) = Kokkos::complex<double>(infinity<double>::value, 0);
h_ref_cbk0(1) =
Expand Down Expand Up @@ -1198,6 +1200,7 @@ struct TestComplexBesselI0K0Function {
h_ref_cbk0(23) = Kokkos::complex<double>(1.413897840559108e-27, 0);
h_ref_cbk0(24) =
Kokkos::complex<double>(1.413897840559108e-27, -1.851678917759592e+25);
h_ref_cbk0(25) = Kokkos::complex<double>(9.5496636116079915979, 0.);

// FIXME_HIP Disable the test when using ROCm 5.5 and 5.6 due to a known
// compiler bug
Expand Down

0 comments on commit 6ff0deb

Please sign in to comment.