[SciPy-Dev] Faster expm

[Alex Leach]

Dear List,

I did some work a while back, which needed to calculate matrix exponentials, 
for which I used the EXPOKIT toolkit: 
http://www.maths.uq.edu.au/expokit/download.html

EXPOKIT is written in Fortran, and I found that f2py did a pretty good job of 
creating a Python interface. More recently, I separated out the code and did 
some speed comparisons against scipy.linalg.expm. If I could direct your 
attention to the scipy-user mailing list, you'll find the build method and 
benchmark comparisons up there:
http://mail.scipy.org/pipermail/scipy-user/2012-November/033588.html

Wrapping EXPOKIT for Python has been done before, by a colleague of the 
original author, so I can't take all the credit. A presentation he gave, with 
some usage and build instructions can be found here (Pg 16-22): 
http://sf.anu.edu.au/~mhk900/Python_Workshop/short.pdf

After some feedback from others on the SciPy-user mailing list, I obtained 
permission from the original author, Roger B. Sidje, to distribute EXPOKIT 
with SciPy under the BSD license.

So, I've done a little more work over the last couple of days, in an attempt 
to merge EXPOKIT into SciPy, but of course there's still more work to be done. 
Unfortunately, I've run out of time to do this myself, (I've got to write my 
PhD thesis, which has nothing to do with this) so I'd like to pass over what 
I've done and forget about it, at least for the next 6 months or so.

If you'd be interested in taking this off my hands, please let me know. I've 
had the small dense routines for Padé approximation working for ages, but the 
routines for Krylov projection of general sparse matrices are still giving 
some problems.

Specifically, there's some STOP statements in the Fortran code that quit the 
program without raising exceptions. I haven't been able to get past these, so 
either the signature file I've edited is incorrect,the sparse matrix I've been 
testing with is wrong or too small, or something else is wrong.

I've just made a fork of scipy on github, to which I've just pushed the files 
I've been working on:

NEW:
scipy/linalg/src/expokit.f
scipy/linalg/expokit.pyf.src
scipy/linalg/time_expm.py
scipy/linalg/expowrap.py
MODIFIED:
scipy/linalg/setup.py

expokit.f       - EXPOKIT source code, with some very minor changes (Cf2py 
comments and changed matvec function call, as per above presentation)
expokit.pyf.src - f2py signature file I've been editing from auto-generated 
version.
time_expm.py    - script I was testing everything with, before adding to 
expowrap.py.. Should not be included in distribution.
expowrap.py     - where I was going to put the python wrapper function. Was 
attempting to write a drop-in replacement for scipy.linalg.expm, but the 
sparse routines were giving problems...

If you have any questions, please don't hesitate to ask.

Kind regards,
Alex
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20121115/4333aa61/attachment.html>