# 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, 2012 at 15:29
• @Paul: Sorry, didn't see this before. No, I need the entire matrix. Feb 8, 2012 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. Mar 6, 2012 at 3:37
• We rely on DGPADM, but at Jack Poulson says, there could be a better way. Sep 4, 2012 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). Feb 8, 2012 at 15:56
• No rush. I've actually implemented it before, so I'm pretty confident :-) Feb 8, 2012 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. Feb 8, 2012 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). Feb 8, 2012 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. Feb 9, 2012 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! Feb 8, 2012 at 15:50

The complex eigensolution approach is mathematically correct, but it does more work than necessary. Unfortunately, the improved approach I'm about to describe cannot be implemented with LAPACK calls.

Look at R. C. Ward and L. J. Gray, ACM Trans. Math. Soft. 4, 278, (1978). This describes the software that is available in TOMS algorithm 530, and which you can download from netlib. This describes how to factor the skew symmetric matrix $X$ as

$$X = U D U^T$$

where $U$ is real orthogonal and $D$ is real skew-symmetric and block diagonal. The diagonal subblocks are either $2\times 2$ or $1\times 1$. Because it is block diagonal, you can exponentiate each subblock separately. The $1\times 1$ blocks are zero, and $\exp(0)=1$, so those are trivial. The $2\times 2$ subblocks are done with

$$\exp \begin{pmatrix} 0 & -t \\ t & 0 \end{pmatrix} = \begin{pmatrix} \cos t & -\sin t \\ \sin t & \cos t \end{pmatrix}$$

The exponential matrix that you want is then given by

$$\exp(X) = U \exp(D) U^T$$

I have used this approach in my quantum chemistry codes for several decades and I have never had any problems with any of the software involved.

• Hello @Ron Shepard, and welcome to Computational Exchange SE. Can you edit your second and third equations? They are a little hard to understand. Aug 15, 2014 at 22:19

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. Feb 8, 2012 at 15:19
• @GeoffOxberry: I rewrote it to include the fortran subroutines
– Paul
Feb 8, 2012 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! Feb 8, 2012 at 15:40