Open i2000s opened 10 years ago
I did a line by line translation of Al-Mohy and Higham's MATLAB implementation of their algorithm, but I am missing their fast 1-norm estimation routine. They also describe that algorithm in a paper, but I have not had the time to try to implement it. If I recall correctly, I just use Julia's 1-norm, and I still see speedups in matrix exponentiation of sparse matrices.
The code can be found here.
@marcusps : That would very be interesting. We should certainly do some direct benchmarking. In acroy/QuDOS.jl#1 we speculated that the dramatic speedup of Al-Mohy & Higham compared to Expokit was partly due to memory allocations required for the Arnoldi procedure. In Julia we can mostly avoid those, so the comparison would be fairer.
BTW, do you know what kind of license the MATLAB implementation has?
It's a BSD license.
Thanks! I will try to have a look.
As started in another thread, the expmv() function can also be written in the Al-Mohy's algorithm (A. H. Al-Mohy and N. J. Higham. Computing the action of the matrix exponential, with an application to exponential integrators. MIMS EPrint 2010.30, The University of Manchester, 2010; SIAM J. Sci. Comput.), beyond the Krylov, Chebyshev and other usual methods. The MATLAB implementation of the Al-Mohy algorithm (with intensive for-loops) and a C rewritten can be found here and here. It would be good to compare those different methods in Julia and provide parametric interface to competitive algorithms. See also the roadmap of IterativeSolvers.jl.