I’m benchmarking an algorithm 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, `Δ^ex`

where `Δ`

is a float and `ex=1/n`

for an integer `n`

.

The full benchmark script, `profile_propagate.jl`

, is meant to be `include`

d in the test-REPL of `QuantumPropagators.jl`

(`make devrepl`

or just `julia --project=test`

from a checkout of that repo).

The routine `extend_leja!`

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?