-
Notifications
You must be signed in to change notification settings - Fork 407
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add math special functions: erf, erfcx, expint1, Bessel functions, Hankel functions #3920
Conversation
I have added the implementations of these functions and their unit tests. All checks passed except HIP tests in jenkin-ci check, where all 0s are returned for the Bessel functions and Hankel functions. |
Which compiler are you using on caraway and tulip? There are know bugs in the compiler used for CI. Once we merged #4025, the CI should used a better compiler. |
e1 = -infinity<RealType>::value; | ||
} | ||
else if (x == 0.0) { | ||
e1 = infinity<RealType>::value; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Shouldn't it return -infty for x=0 ?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@dalg24 I use the definition of E1(x) while std::expint uses the definition of Ei(x). So I changed the function name to expint1
to avoid confusion.
|
||
//! Compute exponential integral E1(x) (x > 0). | ||
template <class RealType> | ||
KOKKOS_INLINE_FUNCTION RealType expint(const RealType& x) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why not taking argument by value
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@dalg24 I see in Kokkos_Complex.hpp, taking argument by reference is used. Please explain why we should take argument by value.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Because this overload is for floating point numbers.
https://en.cppreference.com/w/cpp/numeric/special_functions/expint
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@dalg24 I see. Thanks.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@dalg24 I have changed it to taking argument by value.
|
||
//! Compute error function erf(z) for z=cmplx(x,y). | ||
template <class CmplxType> | ||
KOKKOS_INLINE_FUNCTION CmplxType cerf(const CmplxType& z) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Consider constraining the overload set for complex numbers.
KOKKOS_INLINE_FUNCTION CmplxType cerf(const CmplxType& z) { | |
template <class T> | |
KOKKOS_INLINE_FUNCTION Kokkos::complex<T> cerf(const Kokkos::complex<T>& z) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why did you prefix the name with "c"?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@dalg24 The names are now erf
and erfcx
.
//! Compute scaled complementary error function erfcx(z)=exp(z^2)*erfc(z) | ||
//! for z=cmplx(x,y). | ||
template <class CmplxType> | ||
KOKKOS_INLINE_FUNCTION CmplxType cerfcx(const CmplxType& z) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ok to add erfcx
but we need to provide an implementation for scalar types as well
erfcx(FloatingPoint)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@dalg24 I have added an implementation for scalar type. For now, I do not use the std::erfc() for erfcx()=exp(x^2)*erfc(x)
because it is going to overflow quickly (e.g. erfcx(35)=NaN). Instead, I use erfcx(Kokkos::complex<RealType>(x, 0.0))
for now. I will have an implementation of erfcx(real)
in the future.
//! Compute Bessel function J0(z) of the first kind of order zero | ||
//! for a complex argument | ||
template <class CmplxType, class RealType, class IntType> | ||
KOKKOS_INLINE_FUNCTION CmplxType cbesselj0(const CmplxType& z, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Implement cyl_bessel_j
and error out if the order argument is not 0 (or 1)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@dalg24 Since we decided to use cyl_bessel_j0
and cyl_bessel_j1
for order 0 and order 1, respectively, so I do not have cyl_bessel_j
for now.
//! Compute Bessel function Y1(z) of the second kind of order one | ||
//! for a complex argument | ||
template <class CmplxType, class RealType, class IntType> | ||
KOKKOS_INLINE_FUNCTION CmplxType cbessely1(const CmplxType& z, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Rename cyl_bessel_y1
for consistency with the C++ math functions
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@dalg24 I have changed all the Bessel function names accordingly.
@dalg24 I have addressed all your comments. Please let me know if you have additional comments. Thanks. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'd love for Damien to approve but not sure he got access right now.
@crtrott Thanks, Christian. Should we wait for Damien to merge it? |
I will check if he is able to review if not i will go ahead |
This PR will add mathematical special functions needed by GEMMA's deep slot code.