# Matrix exponential of a skew-Hermitian matrix with fortran 95 and LAPACK

I'm just getting tucked into fortran 95 for some quantum mechanics simulations. Honestly, I've been spoiled by Octave so I've taken matrix exponentiation for granted. Given a (small, $n\leq 36$) skew-Hermitian matrix of size $n\times n$, what is the most efficient way of using LAPACK to solve this problem? I'm not using the LAPACK95 wrapper, just direct calls to LAPACK.

-
Do you need the matrix exponential by itself, or do you need the matrix exponential multiplied by a vector? – Paul Feb 8 '12 at 15:29
@Paul: Sorry, didn't see this before. No, I need the entire matrix. – Mark S. Everitt Feb 8 '12 at 15:52
Why would someone downvote this question? If you downvote, please leave a reason in the comments! Perhaps the question can be improved this way. – Mark S. Everitt Mar 6 '12 at 3:37
We rely on DGPADM, but at Jack Poulson says, there could be a better way. – Mike Dunlavey Sep 4 '12 at 14:39

Matrix exponentials of skew-Hermitian matrices are cheap to compute:

Suppose $A$ is your skew-Hermitian matrix, then $iA$ is Hermitian, and via zheevd and friends you can get the decomposition

$$iA = U \Lambda U^H,$$

where $U$ is the unitary eigenvector matrix and $\Lambda$ is real and diagonal. Then, trivially,

$$A = U (-i \Lambda) U^H.$$

Once you have $U$ and $\Lambda$, it is easy to compute

$$\exp(A)=\exp(U (-i\Lambda) U^H)=U \exp(-i\Lambda) U^H$$

by first exponentiating the eigenvalues, setting $B := U$ via zcopy, performing $B := B \exp(-i \Lambda)$ by running zscal on each column with an exponentiated eigenvalue, and finally setting your result to

$$\exp(A) := B U^H$$

via zgemm.

-
Thanks! I missed an obvious trick there with the $i$. You've put me onto the specific LAPACK subroutines that I need, so extra thanks for that. I won't mark this as correct yet (want to test it out first). – Mark S. Everitt Feb 8 '12 at 15:56
No rush. I've actually implemented it before, so I'm pretty confident :-) – Jack Poulson Feb 8 '12 at 16:02
This is going to be one of those magical bits of code I use all over the place. For what it's worth, I'll also put a thank you in a comment line that probably nobody else will ever see. – Mark S. Everitt Feb 8 '12 at 16:04
@JackPoulson: Well played, sir. This is what I get for picking a major that doesn't believe in imaginary numbers (other than in eigenvalues). – Geoff Oxberry Feb 8 '12 at 18:44
@JackPoulson: It works beautifully. Thanks again for this. Especially the zscal bit. I had most of the rest of the code in another subroutine, but this was something I'd overlooked. – Mark S. Everitt Feb 9 '12 at 2:54

Since I'm on my phone, I can't link things easily, and will add links later. You'll probably want to look at the paper "19 Dubious Ways to Calculate the Matrix Exponential", the Fortran library EXPOKIT, Jitse Niesen's paper on Krylov methods for calculating the Matrix exponential, and some of Nick Higham's recent papers on matrix exponentials. It's more common to need the product of a matrix exponential and a vector than the matrix exponential alone, and here, Krylov methods can be quite helpful. For smaller, dense matrices like the ones you describe, PadÃ© methods might be better, but I've had a lot of success with Krylov methods when used inside exponential methods for numerical integration of ODEs.

-
 Thanks. I'm aware of 19 dubious ways, and also expokit, but some of the people I'm working with are in industry, so I want to avoid it for copyright reasons. I'm keen on implementing it with LAPACK/BLAS since I already link to these libraries. One thing though; I do need the matrix exponential itself. I'm working on a variant of quantum process tomography, and the process in question is embodied by the matrix. Later I'll be dealing with an integrator in combination with this matrix exponential, which is when it gets really interesting! – Mark S. Everitt Feb 8 '12 at 15:50

If all you need is the matrix exponential multiplied by a vector, then this fortran subroutine may be of some use to you. It computes:

$(e^A)v$

where v is a vector, and A is a regular hermitian matrix. It is a subroutine from the EXPOKIT library

Otherwise, you may want to consider this subroutine, which works for any general complex matrix A.

-
 That does not look like a reference to Fortran libraries. – Geoff Oxberry♦ Feb 8 '12 at 15:19 @GeoffOxberry: I rewrote it to include the fortran subroutines – Paul♦ Feb 8 '12 at 15:34 @Paul: No good I'm afraid. What I'm doing is an all-matrix variation on process tomography. In addition skew-Hermitian! – Mark S. Everitt Feb 8 '12 at 15:40 I appreciate that you rewrote your answer, but based on the edit trail, it looks like you completely changed your answer, took elements of my chronologically earlier answer, and added links. – Geoff Oxberry♦ Feb 8 '12 at 15:43 @GeoffOxberry: On the contrary... My results came independently of yours, but you posted before I got a chance to re-write my answer:) – Paul♦ Feb 8 '12 at 15:46