# Extremely slow n-th root of float

**URL:** <https://discourse.julialang.org/t/extremely-slow-n-th-root-of-float/78534>\
**Category:** Performance\
**Tags:** question, profiling\
**Created:** [March 26, 2022, 11:00pm UTC](https://discourse.julialang.org/t/extremely-slow-n-th-root-of-float/78534 "2022-03-26T23:00:02Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![goerz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerz/32/3269_2.png) [@goerz](https://discourse.julialang.org/u/goerz)\
**Post date:** [March 26, 2022, 11:00pm UTC](https://discourse.julialang.org/t/extremely-slow-n-th-root-of-float/78534/1 "2022-03-26T23:00:02Z")

</div>

I’m benchmarking [an algorithm](https://github.com/JuliaQuantumControl/QuantumPropagators.jl/blob/master/src/newton.jl) translated from Fortran, doing something similar to `Expokit.expmv` – evaluate the result of exponentiation a matrix and applying it to a vector, by expanding the `exp` into a polynomial series.

It’s not performing as well as it should (by a factor of 5-10). Profiling shows that strangely, 50% of the entire runtime is spent in a [line that simply calculates the n-th root of a float](https://github.com/JuliaQuantumControl/QuantumPropagators.jl/blob/542a5a9259b691ff4a2ea40d4ccf5cc2e9784ada/src/newton.jl#L124), `Δ^ex` where `Δ` is a float and `ex=1/n` for an integer `n`.

The full benchmark script, [`profile_propagate.jl`](https://gist.github.com/goerz/48c29db842742cc872f798c4711840d9), is meant to be `include`d in the test-REPL of [`QuantumPropagators.jl`](https://github.com/JuliaQuantumControl/QuantumPropagators.jl) (`make devrepl` or just `julia --project=test` from a checkout of that repo).

The routine [`extend_leja!`](https://github.com/JuliaQuantumControl/QuantumPropagators.jl/blob/542a5a9259b691ff4a2ea40d4ccf5cc2e9784ada/src/newton.jl#L86) where this occurs in an almost 1-to-1 transcription of the original Fortran code. When I profile that Fortran code, the corresponding routine is completely negligible. The expected behavior is for this algorithm to be dominated by the matrix-vector products, not by some scalar operations for calculating the expansion coefficients. (There’s also a diagonalization of a small Hessenberg matrix that takes much longer in the Julia code than it should, but I’ll look at that after I can figure out what’s going on with the calculation of the coefficients).

Any ideas?

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [March 26, 2022, 11:03pm UTC](https://discourse.julialang.org/t/extremely-slow-n-th-root-of-float/78534/2 "2022-03-26T23:03:08Z")

</div>

Are you on Linux? `^` is known to be very slow on Linux in Julia 1.6 and 1.7, but it’s better in \<= 1.5 and \>= 1.8.

It doesn’t have this problem on Windows or Mac.  
But the problem is because it’s calling the system `^` function, which means Fortran should be experiencing slow `^` as well…

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [March 26, 2022, 11:20pm UTC](https://discourse.julialang.org/t/extremely-slow-n-th-root-of-float/78534/3 "2022-03-26T23:20:25Z")

</div>

Huh:

```julia
julia> const libdpow = "/home/chriselrod/Documents/progwork/fortran/libdpow.so"
"/home/chriselrod/Documents/progwork/fortran/libdpow.so"

julia> dpow(x,y) = @ccall libdpow.dpow(x::Float64, y::Float64)::Float64
dpow (generic function with 1 method)

julia> syspow(x,y) = ccall(:pow, Float64, (Float64,Float64), x, y)
syspow (generic function with 1 method)

julia> @btime $(Ref(2.3))[]^$(Ref(1.2))[]
  28.336 ns (0 allocations: 0 bytes)
2.716898432499149

julia> @btime dpow($(Ref(2.3))[],$(Ref(1.2))[])
  15.797 ns (0 allocations: 0 bytes)
2.716898432499149

julia> @btime syspow($(Ref(2.3))[],$(Ref(1.2))[])
  70.757 ns (0 allocations: 0 bytes)
2.716898432499149

julia> versioninfo()
Julia Version 1.9.0-DEV.229
Commit 7cde4be23d* (2022-03-22 02:44 UTC)
Platform Info:
  OS: Linux (x86_64-redhat-linux)
  CPU: 28 × Intel(R) Core(TM) i9-9940X CPU @ 3.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, skylake-avx512)
  Threads: 28 on 28 virtual cores

```

```fortran
module pow
use ISO_C_BINDING
  
implicit none

contains

real(C_double) function dpow(x, y) bind(C, name="dpow")
    real(C_double), value, intent(in) :: x, y
    dpow = x**y
    return
end function dpow

end module pow

```

I compiled with

```julia
gfortran -O3 -shared -fPIC pow.f90 -o libdpow.so

```

@Oscar_Smith  
Interestingly, I didn’t even compile the Fortran with `-march=native` (or with `-Ofast`)!

Related:  
[https://github.com/JuliaLang/julia/pull/44717](https://github.com/JuliaLang/julia/pull/44717)

If you’re curious, the Fortran assembly

```nohighlight
dpow:
.LFB0:
	.cfi_startproc
	jmp	pow@PLT
	.cfi_endproc

```

I guess it is jumping to a different `pow` than we’re calling with `ccall(:pow, ...)` above…

Also, compiling with `ifort -fast` instead, I get

```julia
julia> @btime dpow($(Ref(2.3))[],$(Ref(1.2))[])
  12.671 ns (0 allocations: 0 bytes)
2.716898432499149

```

But presumably this function is less accurate.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [March 27, 2022, 1:08am UTC](https://discourse.julialang.org/t/extremely-slow-n-th-root-of-float/78534/4 "2022-03-27T01:08:38Z")

</div>

Can you test the accuracy of this? I’m curious how good it is.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [March 27, 2022, 1:33am UTC](https://discourse.julialang.org/t/extremely-slow-n-th-root-of-float/78534/5 "2022-03-27T01:33:34Z")

</div>

I’ll get around to it, working on something else at the moment.  
But if you’re curious, you can install the Intel compilers through `apt`:  
[https://www.intel.com/content/www/us/en/develop/documentation/installation-guide-for-intel-oneapi-toolkits-linux/top/installation/install-using-package-managers/apt.html](https://www.intel.com/content/www/us/en/develop/documentation/installation-guide-for-intel-oneapi-toolkits-linux/top/installation/install-using-package-managers/apt.html)  
They’re all free as in beer now (but not as in freedom).

---

<div class="post-metadata">

**Author:** ![goerz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerz/32/3269_2.png) [@goerz](https://discourse.julialang.org/u/goerz)\
**Post date:** [March 27, 2022, 4:37am UTC](https://discourse.julialang.org/t/extremely-slow-n-th-root-of-float/78534/6 "2022-03-27T04:37:39Z")

</div>

Yes, I’m on Linux (Ubuntu), currently using Julia 1.7.2, and I’m comparing with `ifort`-compiled Fortran code (with `-O3`).

I just tried running it with the `1.8.0-beta1`, and I can confirm that `^` seems to be about twice as fast as on `1.7.2`, which matches the comment in [https://github.com/JuliaLang/julia/pull/44717#issue-1178579501](https://github.com/JuliaLang/julia/pull/44717#issue-1178579501). This is from eyeballing the `ProfileSVG` output; `Pprof` and `StatProfilerHTML` which would give more quantitative info both still crash on `1.8`.

I also tried it on my Macbook (with `1.7`), and it looked like the _relative_ amount of time spent in `^` was indeed also shorter there, although overall performance was abysmal. I’m not sure I can trust the profiler there, it’s showing the majority of time spent in top-level functions named `.L10`, `.L999`, and similar that I have no idea what they mean. I’d guess this maybe has to do with Rosetta (I haven’t had time yet to really look into the M1 support in 1.7). It doesn’t matter _that_ much, since anything serious is going to run on the Linux workstation.

The overall algorithm is still much slower than it should be, but now that Hessenberg diagonalization is the slowest part, so I’ll have a look at that next. Also, I might have to be a bit more thorough with my Fortran benchmarking, maybe try `gfortran` instead of `ifort`. For one of [the other methods in `QuantumPropagators`](https://github.com/JuliaQuantumControl/QuantumPropagators.jl/blob/master/src/cheby.jl) that’s far less involved (but limited to Hermitian matrices), I got the Julia code to be exactly as fast as the Fortran with compiled `gfortran` (which was a factor of two slower than `ifort`). I was hoping to get something similar here.

I’ll keep twiddling the code, and maybe also see how some of the existing Julia solutions like `ExpoKit` perform (I’m just implementing my old true-and-tested methods first to have something to fall back on that I know exactly what it’s doing)

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [March 27, 2022, 4:40am UTC](https://discourse.julialang.org/t/extremely-slow-n-th-root-of-float/78534/7 "2022-03-27T04:40:54Z")

</div>

[https://github.com/JuliaLang/julia/pull/44717#issue-1178579501](https://github.com/JuliaLang/julia/pull/44717#issue-1178579501) isn’t in 1.8. That will be a future speedup in 1.9. The part that did go into 1.8 is [https://github.com/JuliaLang/julia/pull/42271](https://github.com/JuliaLang/julia/pull/42271)

---

<div class="post-metadata">

**Author:** ![goerz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerz/32/3269_2.png) [@goerz](https://discourse.julialang.org/u/goerz)\
**Post date:** [March 27, 2022, 4:41am UTC](https://discourse.julialang.org/t/extremely-slow-n-th-root-of-float/78534/8 "2022-03-27T04:41:40Z")

</div>

Ah! Should I try the Nightly?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [March 27, 2022, 4:45am UTC](https://discourse.julialang.org/t/extremely-slow-n-th-root-of-float/78534/9 "2022-03-27T04:45:13Z")

</div>

Note that the PR isn’t merged yet. To see the change, you would have to build the PR from source. Also I’m planning on making a number of changes before this merges (that hopefully will have small to no performance impact for this) before it merges. Also, in the next few weeks I should have a PR that speeds up `exp` by 1-2 ns that will also help this.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [March 27, 2022, 5:13am UTC](https://discourse.julialang.org/t/extremely-slow-n-th-root-of-float/78534/10 "2022-03-27T05:13:19Z")

</div>

Also…

```julia
julia> const libm6 = "/usr/lib64/libm.so.6"
"/usr/lib64/libm.so.6"

julia> syspow(x,y) = @ccall libm6.pow(x::Float64, y::Float64)::Float64
syspow (generic function with 1 method)

julia> @btime syspow($(Ref(2.3))[],$(Ref(1.2))[])
  15.675 ns (0 allocations: 0 bytes)
2.716898432499149

```

So, actually, system libm is fast on Linux.

If it isn’t the system libm, which `pow` is Julia actually calling on 1.6 and 1.7, and why is it different on Windows, Mac, and Linux???

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [March 27, 2022, 5:18am UTC](https://discourse.julialang.org/t/extremely-slow-n-th-root-of-float/78534/11 "2022-03-27T05:18:14Z")

</div>

Wait, that’s really weird. Are 1.6 and 1.7 calling an LLVM specific pow somehow?

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [March 27, 2022, 5:23am UTC](https://discourse.julialang.org/t/extremely-slow-n-th-root-of-float/78534/12 "2022-03-27T05:23:18Z")

</div>

Apparently, but why would that be different on each OS?