Comparing Fortran and Julia's Bessel function performance

40

Thank you @rgba and @jacobwilliams, I thought I had also tried with implicit none and declaring an output variable type and didn’t see a change, but I must have accidentally kept around an old shared library or something, because doing this again, I do indeed see the full double precision output. My bad!

I updated my code with the proper output type annotations and it agrees with the nicer, more modern version @jacobwilliams showed, so I’m using their version here.

Here is what I see for the accuracy of the Fortran BESSEL_J1 relative to an arbitrary precision calculuation, and compared to Bessels.bessel_j1 for 1,000,000 random samples between 0 and 200.
acc

So the algorithm Fortran is using does indeed have a tighter error distribution near zero relative to Bessels, but also a significantly longer error tail. That makes for a rather interesting comparison where different people may value the different characteristics differently.

Click for full reproducible code 
#+begin_src fortran
module bess
  use iso_fortran_env, only: wp => real64
  implicit none
  contains
    pure function f90_bessel_j1(x) result(r)
      real(wp),intent(in) :: x
      real(wp) :: r
      r = bessel_j1(x)
    end function f90_bessel_j1
  end module bess
#+end_src

shell> gfortran bessel_j1.f90 -o bessel_j1_fort.so -shared -fPIC

#+begin_src julia
using Bessels, ArbNumerics, Libdl, Plots

setworkingprecision(ArbFloat, 500)
setextrabits(128)

dlopen("/home/mason/bessel_j1_fort.so") do lib
    fptr = dlsym(lib, :__bess_MOD_f90_bessel_j1)
    f90_besselj1(x) = ccall(fptr, Float64, (Ref{Float64},), x)
    
    xs = 200 .* rand(1_000_000)
    
    v_bessels = Bessels.besselj1.(xs)
    v_fort    = f90_besselj1.(xs)
    v_arb     = ArbNumerics.besselj1.(ArbFloat.(xs))
    
    resid_bessels = float.(v_bessels .- v_arb) ./ eps()
    resid_fort    = float.(v_fort .- v_arb) ./ eps()
    
    plot(xlabel="error / eps", ylabel="counts")
    histogram!(resid_bessels, nbins=41, yscale=:log10, ylims=(1, Inf), label="Bessels.jl besselj1")
    histogram!(resid_fort, fillalpha=0.5, nbins=41, yscale=:log10, ylims=(1, Inf), label="Fortran BESSEL_J1")
end
#+end_src
Click for system details 
shell> gfortran --version
GNU Fortran (GCC) 12.2.0
Copyright (C) 2022 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

julia> versioninfo()
Julia Version 1.8.0
Commit 5544a0fab7* (2022-08-17 13:38 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: 12 × AMD Ryzen 5 5600X 6-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, znver3)
  Threads: 6 on 12 virtual cores
Environment:
  JULIA_NUM_THREADS = 6
3 Likes