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,
Δ is a float and
ex=1/n for an integer
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).