This notebook compares the wrapped Fortran SLATEC implementation contained in py3nj
with the native Julia WignerFamilies.jl
. I use very high quantum numbers for this test, in order to amortize the cost of calling SLATEC with some Python overhead.
using WignerFamilies
using PyCall
using BenchmarkTools
py3nj = pyimport("py3nj")
j₂, j₃, m₂, m₃ = 40000, 40000, -2, 2
w = WignerF(Float64, j₂, j₃, m₂, m₃)
w3j = get_wigner_array(w)
print("WignerFamilies.jl: ")
@btime WignerFamilies.wigner3j_f!(w, w3j)
print("SLATEC wrapper: ")
@btime two_l1, three_j = py3nj.wigner.drc3jj(2j₂, 2j₃, 2m₂, 2m₃);
WignerFamilies.jl: 619.900 μs (0 allocations: 0 bytes) SLATEC wrapper: 1.637 ms (83 allocations: 2.44 MiB)
This is perhaps a more fair test, since the wrapped SLATEC is allocating each time.
print("WignerFamilies.jl: ")
@btime WignerFamilies.wigner3j_f(Float64, j₂, j₃, m₂, m₃);
print("SLATEC wrapper: ")
@btime two_l1, three_j = py3nj.wigner.drc3jj(2j₂, 2j₃, 2m₂, 2m₃);
WignerFamilies.jl: 652.100 μs (3 allocations: 625.17 KiB) SLATEC wrapper: 1.671 ms (83 allocations: 2.44 MiB)
WignerFamilies.jl uses a special case for m_i = 0
that gives it an extra factor of 2.
j₂, j₃, m₂, m₃ = 40000, 40000, 0, 0
print("WignerFamilies.jl: ")
@btime WignerFamilies.wigner3j_f(Float64, j₂, j₃, m₂, m₃);
print("SLATEC wrapper: ")
@btime two_l1, three_j = py3nj.wigner.drc3jj(2j₂, 2j₃, 2m₂, 2m₃);
WignerFamilies.jl: 387.400 μs (3 allocations: 625.17 KiB) SLATEC wrapper: 1.673 ms (83 allocations: 2.44 MiB)