Bessel function of higher order #5897

Yaraslaut · 2023-02-21T08:48:06Z

This PR adds Bessel functions of higher order (#5637)

TODO:

update unit test so they cover all the code and ckeck precision for higher orders as well

dalg24-jenkins · 2023-02-21T08:48:08Z

Can one of the admins verify this patch?

PhilMiller · 2023-02-21T16:53:18Z

If you end up rebasing for whatever reason, there's a typo in the first commit message, 'fucntions'

PhilMiller · 2023-02-21T16:54:20Z

OK to test

Yaraslaut · 2023-02-21T17:52:39Z

If you end up rebasing for whatever reason, there's a typo in the first commit message, 'fucntions'

oh, will then squash

mhoemmen · 2023-02-21T18:18:21Z

I'm not a Kokkos dev, but would you consider calling the existing mathematical special functions where appropriate?

Kokkos requires C++17, which adds mathematical special functions
CUDA has device implementations of the mathematical special functions

Not using existing implementations where possible risks pessimizing both performance and accuracy.

Yaraslaut · 2023-02-23T07:34:21Z

I'm not a Kokkos dev, but would you consider calling the existing mathematical special functions where appropriate?
* Kokkos requires C++17, which adds [mathematical special functions](https://en.cppreference.com/w/cpp/numeric/special_functions)

* CUDA has [device implementations](https://docs.nvidia.com/cuda/cuda-math-api/group__CUDA__MATH.html#group__CUDA__MATH) of the mathematical special functions
Not using existing implementations where possible risks pessimizing both performance and accuracy.

Hi, here is small example of how this functions can be called

// adl
void t()
{
  using namespace std;
  using namespace Kokkos::Experimental;

  auto d{2.0};
  cyl_bessel_j(2,d);

  auto std_complex = std::complex<double>{2.0,2.0};
  //cyl_bessel_j(2,std_complex); // error

  auto kokkos_complex = Kokkos::complex<double>{2.0,2.0};
  //cyl_bessel_j(2,kokkos_complex); //error
  cyl_bessel_j<decltype(kokkos_complex),double,int>(2,kokkos_complex);
}


int main(int argc, char* argv[]) {
  Kokkos::ScopeGuard guard(argc, argv);
  auto d = 2.0;
  auto kokkos_complex = Kokkos::complex<double>(2.0,2.0);
  auto std_complex = std::complex<double>(2.0,2.0);
  std::cyl_bessel_i(1, d);
  //std::cyl_bessel_i(1, kokkos_complex);  error
  //std::cyl_bessel_i(1, std_complex); error

  //Kokkos::Experimental::cyl_bessel_j1<decltype(d),double,int>(d); //error
  Kokkos::Experimental::cyl_bessel_j1<decltype(kokkos_complex),double,int>(kokkos_complex);
  //Kokkos::Experimental::cyl_bessel_j1<decltype(std_complex),double,int>(std_complex); //error
}

bessel functions from kokkos require explicit template arguments since they can not deduce all types from one argument only.
functions from std have real arguments only so if user overlap namespaces and call bessel function directly without namespace specification it will end up calling std functoin ( void t() function )
kokkos functions can take as argument only Kokkos::complex<> type

I think it is a good idea to have bessel functions for real types in Kokkos namespace so users can just call Kokkos::cyl_besselj(2,2.0) but that another topic of how properly expose special functions from <cmath> in Kokkos namespace so ADL can call functions

mhoemmen · 2023-02-27T15:46:44Z

My comment relates to the implementation of these functions, not their interface. Anyway, please carry on : - )

fnrizzi · 2023-09-11T13:28:30Z

@Yaraslaut
I just started looking at this, thanks for the PR !

first comment:

auto kokkos_complex = Kokkos::complex<double>{2.0,2.0};
cyl_bessel_j<decltype(kokkos_complex),double,int>(2,kokkos_complex);

I think it would be useful to allow doing this:

auto kokkos_complex = Kokkos::complex<double>{2.0,2.0};
cyl_bessel_j(2,kokkos_complex);

where this could be handled as:

template <class RealType, class IntType>
KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_j(const IntType a,
                                              const Kokkos::complex<RealType> & z,
                                              const RealType& joint_val = 25,
                                              const IntType& bw_start   = 70) {

second comment: I agree with @mhoemmen: ideally we can expose as rich of an API as needed (for example, also to allow users to avoid having to specify all templates like in my example above), but internally we should call native libraries if available and if nothing is avaible have a fallback implementation (effectively what you have in this PR could be that)

Yaraslaut · 2023-09-11T13:59:30Z

@Yaraslaut I just started looking at this, thanks for the PR !

1. first comment:

auto kokkos_complex = Kokkos::complex<double>{2.0,2.0};
cyl_bessel_j<decltype(kokkos_complex),double,int>(2,kokkos_complex);

I think it would be useful to allow doing this:

auto kokkos_complex = Kokkos::complex<double>{2.0,2.0};
cyl_bessel_j(2,kokkos_complex);

where this could be handled as:

template <class RealType, class IntType>
KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_j(const IntType a,
                                              const Kokkos::complex<RealType> & z,
                                              const RealType& joint_val = 25,
                                              const IntType& bw_start   = 70) {

2. second comment: I agree with @mhoemmen: ideally we can expose as rich of an API as needed (for example, also to allow users to avoid having to specify all templates like in my example above), but internally we should call native libraries if available and if nothing is avaible have a fallback implementation (effectively what you have in this PR could be that)

Thanks for reviewing this PR, I will try to address this issues, and let you know once it is ready for further review.

masterleinad · 2023-09-18T15:30:57Z

Retest this please.

masterleinad · 2023-10-06T17:13:22Z

Looks like the the test fails when using std functions:

4: [ RUN      ] hip.mathspecialfunc_cbesselj2y2
4: /var/jenkins/workspace/Kokkos/core/unit_test/TestMathematicalSpecialFunctions.hpp:2132: Failure
4: Expected: (Kokkos::abs(h_cbj2(i) - h_ref_cbj2(i))) <= (Kokkos::abs(h_ref_cbj2(i)) * 1e-13), actual: 0.947984 vs 1.22779e-13
4: /var/jenkins/workspace/Kokkos/core/unit_test/TestMathematicalSpecialFunctions.hpp:2132: Failure
4: Expected: (Kokkos::abs(h_cbj2(i) - h_ref_cbj2(i))) <= (Kokkos::abs(h_ref_cbj2(i)) * 1e-13), actual: 0.947984 vs 1.22779e-13
4: /var/jenkins/workspace/Kokkos/core/unit_test/TestMathematicalSpecialFunctions.hpp:2132: Failure
4: Expected: (Kokkos::abs(h_cbj2(i) - h_ref_cbj2(i))) <= (Kokkos::abs(h_ref_cbj2(i)) * 1e-13), actual: 0.947984 vs 1.22779e-13
4: /var/jenkins/workspace/Kokkos/core/unit_test/TestMathematicalSpecialFunctions.hpp:2132: Failure
4: Expected: (Kokkos::abs(h_cbj2(i) - h_ref_cbj2(i))) <= (Kokkos::abs(h_ref_cbj2(i)) * 1e-13), actual: 0.947984 vs 1.22779e-13
4: [  FAILED  ] hip.mathspecialfunc_cbesselj2y2 (0 ms)
4: [ RUN      ] hip.mathspecialfunc_cbesseli2k2
4: [       OK ] hip.mathspecialfunc_cbesseli2k2 (0 ms)
4: [ RUN      ] hip.mathspecialfunc_cbesselh12h22
4: /var/jenkins/workspace/Kokkos/core/unit_test/TestMathematicalSpecialFunctions.hpp:2491: Failure
4: Expected: (Kokkos::abs(h_cbh12(i) - h_ref_cbh12(i))) <= (Kokkos::abs(h_ref_cbh12(i)) * 1e-13), actual: 0.432749 vs 2.41179e-13
4: /var/jenkins/workspace/Kokkos/core/unit_test/TestMathematicalSpecialFunctions.hpp:2491: Failure
4: Expected: (Kokkos::abs(h_cbh12(i) - h_ref_cbh12(i))) <= (Kokkos::abs(h_ref_cbh12(i)) * 1e-13), actual: 0.432749 vs 2.50098e-13
4: /var/jenkins/workspace/Kokkos/core/unit_test/TestMathematicalSpecialFunctions.hpp:2491: Failure
4: Expected: (Kokkos::abs(h_cbh12(i) - h_ref_cbh12(i))) <= (Kokkos::abs(h_ref_cbh12(i)) * 1e-13), actual: 0.226039 vs 5.11872e-14
4: /var/jenkins/workspace/Kokkos/core/unit_test/TestMathematicalSpecialFunctions.hpp:2491: Failure
4: Expected: (Kokkos::abs(h_cbh12(i) - h_ref_cbh12(i))) <= (Kokkos::abs(h_ref_cbh12(i)) * 1e-13), actual: 0.226039 vs 5.11872e-14
4: /var/jenkins/workspace/Kokkos/core/unit_test/TestMathematicalSpecialFunctions.hpp:2491: Failure
4: Expected: (Kokkos::abs(h_cbh12(i) - h_ref_cbh12(i))) <= (Kokkos::abs(h_ref_cbh12(i)) * 1e-13), actual: 0.00343646 vs 1.66666e-14
4: /var/jenkins/workspace/Kokkos/core/unit_test/TestMathematicalSpecialFunctions.hpp:2491: Failure
4: Expected: (Kokkos::abs(h_cbh12(i) - h_ref_cbh12(i))) <= (Kokkos::abs(h_ref_cbh12(i)) * 1e-13), actual: 0.00343646 vs 1.66666e-14
4: /var/jenkins/workspace/Kokkos/core/unit_test/TestMathematicalSpecialFunctions.hpp:2497: Failure
4: Expected: (Kokkos::abs(h_cbh22(i) - h_ref_cbh22(i))) <= (Kokkos::abs(h_ref_cbh22(i)) * 1e-13), actual: 0.182976 vs 7.73567e-15
4: /var/jenkins/workspace/Kokkos/core/unit_test/TestMathematicalSpecialFunctions.hpp:2497: Failure
4: Expected: (Kokkos::abs(h_cbh22(i) - h_ref_cbh22(i))) <= (Kokkos::abs(h_ref_cbh22(i)) * 1e-13), actual: 0.182976 vs 7.73567e-15
4: /var/jenkins/workspace/Kokkos/core/unit_test/TestMathematicalSpecialFunctions.hpp:2497: Failure
4: Expected: (Kokkos::abs(h_cbh22(i) - h_ref_cbh22(i))) <= (Kokkos::abs(h_ref_cbh22(i)) * 1e-13), actual: 1.38258e-05 vs 7.45952e-19
4: /var/jenkins/workspace/Kokkos/core/unit_test/TestMathematicalSpecialFunctions.hpp:2497: Failure
4: Expected: (Kokkos::abs(h_cbh22(i) - h_ref_cbh22(i))) <= (Kokkos::abs(h_ref_cbh22(i)) * 1e-13), actual: 1.38258e-05 vs 7.45952e-19
4: /var/jenkins/workspace/Kokkos/core/unit_test/TestMathematicalSpecialFunctions.hpp:2497: Failure
4: Expected: (Kokkos::abs(h_cbh22(i) - h_ref_cbh22(i))) <= (Kokkos::abs(h_ref_cbh22(i)) * 1e-13), actual: 1.81657e-05 vs 6.79118e-19
4: /var/jenkins/workspace/Kokkos/core/unit_test/TestMathematicalSpecialFunctions.hpp:2497: Failure
4: Expected: (Kokkos::abs(h_cbh22(i) - h_ref_cbh22(i))) <= (Kokkos::abs(h_ref_cbh22(i)) * 1e-13), actual: 1.81657e-05 vs 6.79118e-19
4: /var/jenkins/workspace/Kokkos/core/unit_test/TestMathematicalSpecialFunctions.hpp:2497: Failure
4: Expected: (Kokkos::abs(h_cbh22(i) - h_ref_cbh22(i))) <= (Kokkos::abs(h_ref_cbh22(i)) * 1e-13), actual: 1.44286e-05 vs 4.66927e-19
4: /var/jenkins/workspace/Kokkos/core/unit_test/TestMathematicalSpecialFunctions.hpp:2497: Failure
4: Expected: (Kokkos::abs(h_cbh22(i) - h_ref_cbh22(i))) <= (Kokkos::abs(h_ref_cbh22(i)) * 1e-13), actual: 1.44286e-05 vs 4.66927e-19
4: [  FAILED  ] hip.mathspecialfunc_cbesselh12h22 (1 ms)

Rombur · 2023-10-06T17:18:18Z

@masterleinad No I don't think so. I think that's the compiler bug that we discussed earlier in this PR. We just need to disable the failing tests for ROCm 5.5 and 5.6

masterleinad · 2023-10-06T17:41:29Z

@masterleinad No I don't think so. I think that's the compiler bug that we discussed earlier in this PR. We just need to disable the failing tests for ROCm 5.5 and 5.6

Even better!

Yaraslaut · 2023-10-07T05:32:19Z

I think that PR is somewhat ready for review, sorry it took so long to get it working after initial implementation it :(
not sure whom should i ping to review this changes, @fnrizzi ? @Rombur ?

masterleinad

Not a complete review. Just some obeservations.

core/src/Kokkos_MathematicalSpecialFunctions.hpp

core/unit_test/TestMathematicalSpecialFunctions.hpp

masterleinad · 2023-10-17T14:13:04Z

We still need to disable the tests failing for ROCm 5.5 and 5.6.

core/src/Kokkos_MathematicalSpecialFunctions.hpp

masterleinad

Apart from testing the public interface instead of the function in Impl, this looks good to me (pending CI).

masterleinad · 2023-10-23T12:36:37Z

core/unit_test/TestMathematicalSpecialFunctions.hpp

-    d_cby0(i) = Kokkos::Experimental::cyl_bessel_y0<Kokkos::complex<double>,
-                                                    double, int>(d_z(i));
+    d_cbj0(i) =
+        Kokkos::Experimental::Impl::cyl_bessel_j0<Kokkos::complex<double>,


Can you please change all of these to use/test the public interface?

core/unit_test/TestMathematicalSpecialFunctions.hpp

fnrizzi · 2023-10-25T09:08:50Z

core/unit_test/TestMathematicalSpecialFunctions.hpp

+    Kokkos::Experimental::cyl_bessel_k(6, 1.0);
+    Kokkos::Experimental::cyl_bessel_y(6, 1.0);
+    Kokkos::Experimental::cyl_bessel_h1(6, 1.0);
+    Kokkos::Experimental::cyl_bessel_h2(6, 1.0);


what is being tested here? just that it compiles?

Yes, only that it properly deducing templates, i can also add tests for the numerical value of those

I added tests for higher orders for both real and complex arguments, and I was lucky to find out divergene between results taken from octave and julia. So I need to rethink implementation, if anyone have any ideas i will be happy to hear them

core/src/Kokkos_MathematicalSpecialFunctions.hpp

masterleinad · 2023-10-26T18:09:06Z

You'll have to fix the indentation.

Yaraslaut · 2023-10-26T18:15:12Z

You'll have to fix the indentation.

Thanks, for pointing this our,otherwise i would discover it only after a while. I tried to addres all issues, still not sure what to do with large arguments since they require customization of default variables to get proper result

fnrizzi · 2023-11-06T08:41:23Z

making this a draft following the comment: #5897 (comment)

Yaraslaut · 2023-11-23T07:53:17Z

Search for an existing implementation led me to discover implementation of Bessel function by Amos[1]. You can see this is a Fortran code and i am not sure if kokkos wants to link/compile this library, perhaps @fnrizzi or @masterleinad can advice how we can proceed with this PR.

[1] https://netlib.org/amos/
Few other references for the history
[2] https://github.com/JuliaMath/openspecfun/tree/main/amos
[3] COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT BY D. E. AMOS, SAND83-0083, MAY, 1983.
[4] COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
[5] SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-1018, MAY, 1985
[6] PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS. MATH. SOFTWARE, 1986

masterleinad · 2023-11-27T18:23:28Z

Search for an existing implementation led me to discover implementation of Bessel function by Amos[1]. You can see this is a Fortran code and i am not sure if kokkos wants to link/compile this library, perhaps @fnrizzi or @masterleinad can advice how we can proceed with this PR.

I don't think we want to use an external library just to be able to use Bessel functions (and does that even work with the GPU backends?). I would prefer just adopting a suitable implementation and converting it to Kokkos code (adding the appropriate copyrights, of course).

fnrizzi · 2023-11-28T11:23:22Z

Search for an existing implementation led me to discover implementation of Bessel function by Amos[1]. You can see this is a Fortran code and i am not sure if kokkos wants to link/compile this library, perhaps @fnrizzi or @masterleinad can advice how we can proceed with this PR.

I don't think we want to use an external library just to be able to use Bessel functions (and does that even work with the GPU backends?). I would prefer just adopting a suitable implementation and converting it to Kokkos code (adding the appropriate copyrights, of course).

i agree with Daniel!

Yaraslaut · 2023-12-02T09:59:32Z

@masterleinad @fnrizzi existing libraries does not support complex arguments, for example see boost: boostorg/math#689

vqd8a self-requested a review February 21, 2023 16:15

Yaraslaut force-pushed the develop branch 4 times, most recently from 76e2055 to a8ecfc1 Compare February 22, 2023 13:21

vqd8a removed their request for review February 23, 2023 17:00

Yaraslaut force-pushed the develop branch from a8ecfc1 to 8eb9c31 Compare February 27, 2023 09:53

Yaraslaut marked this pull request as ready for review February 27, 2023 10:51

Yaraslaut force-pushed the develop branch 4 times, most recently from 9f95bc2 to 925b36d Compare September 14, 2023 05:47

Yaraslaut closed this Sep 14, 2023

Yaraslaut force-pushed the develop branch from 925b36d to 6a95b5f Compare September 14, 2023 06:09

Yaraslaut reopened this Sep 14, 2023

Yaraslaut force-pushed the develop branch from d233455 to 6e4e6b4 Compare September 15, 2023 05:18

Add Bessel function of higher order

279b07f

Yaraslaut force-pushed the develop branch from 6e4e6b4 to 279b07f Compare September 15, 2023 05:20

Yaraslaut force-pushed the develop branch from 7f9d5e6 to 106fb68 Compare September 22, 2023 09:34

use loop instead of recursion

c004ed7

masterleinad reviewed Oct 7, 2023

View reviewed changes

core/src/Kokkos_MathematicalSpecialFunctions.hpp Outdated Show resolved Hide resolved

core/src/Kokkos_MathematicalSpecialFunctions.hpp Outdated Show resolved Hide resolved

core/unit_test/TestMathematicalSpecialFunctions.hpp Show resolved Hide resolved

dalg24 mentioned this pull request Oct 9, 2023

Add taking the address of math constants to known issues kokkos/kokkos-core-wiki#461

Merged

Hide implementatin details and add deprecation info

63aefa6

Yaraslaut requested a review from masterleinad October 17, 2023 06:47

masterleinad reviewed Oct 17, 2023

View reviewed changes

hide default parameters inside function

33af74f

masterleinad reviewed Oct 23, 2023

View reviewed changes