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.
Tell me more
×
Computational Science Stack Exchange is a question and answer site for
scientists using computers to solve scientific problems. It's 100% free, no registration required.
|
|
|||||||||||
|
|
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. |
|||||||||||||||||
|
|
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. |
|||
|
|
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. |
|||||||||||
|
