# High performance exponential of a skew Hermitian matrix in Fortran 95

I'm new to Fortran, and this is pretty much my first escapade. Below is a function that I wrote which relies on calls to LAPACK. The function is sat in a module with some other functions and works perfectly, but seeing as this is really the workhorse of the program I'm building I want to squeeze it hard for performance. How does it look? Am I doing anything stupid, or would another call to BLAS or LAPACK somewhere improve the performance?

The function takes a Hermitian matrix H, and returns the matrix exponential of the skew-Hermitian matrix -iHt where i is the imaginary number and t is a real number. (This is the solution to the Schrodinger equation for a time independent Hamiltonian.) s, the length of one side of the Hamiltonian, is included for automatic array initialisation rather than using allocations. Typically it'll be less than 10.

function time_indep_schrodinger(s,H,t)
! s : The length of one side of H.
! H : The Hamiltonian.
! t : The time t.
!
! The LAPACK subroutine zheev requires a tridiagonal
! subcopy of H. This array is then transformed into
! a matrix of eigenvectors. This also yields the
! eigenvalues as a 1D array, which are then
! exponentiated before using the matrices of
! eigenvectors to produce the evolution operator.
!
! Finds the eigenvalues of Ht, then uses U exp(-i*eigs) U^H.

integer,    intent(in) :: s
complex*16, intent(in) :: H(s,s)
real*8,     intent(in) :: t
complex*16             :: B(s,s), eigv(s,s), time_indep_schrodinger(s,s), work(2*s-1)
real*8                 :: rwork(3*s-2), eigs(s)
integer                :: info, n, m

! Hermitian matrix to diagonalise.
forall (n=1:s, m=1:s, m>=n) eigv(n,m) = H(n,m)*t

! Get eigenvectors and eigenvalues.
call zheev('V','U',s,eigv,s,eigs,work,2*s-1,rwork,info)

! Scale columns of copied matrix be exponentiated eigenvalues (with -i).
do n=1,s
B(:,n) = eigv(:,n)*exp((0,-1)*eigs(n))
end do

! Finally multiply scaled eigenvectors with conjugate of eigv.
time_indep_schrodinger = matmul(B,conjg(transpose(eigv)))

end function

• Is there a directive to instruct the syntax highlighter on which language to highlight for? Feb 12, 2012 at 12:53
• Turns out that there is, but fortran is not supported. Turned off syntax highlighting for the code above. Feb 13, 2012 at 4:43
• Even though it's the wrong place, you might try posting in StackOverflow. It might get more visibility by fortran followers and if you're lucky, someone will respond w/ some advice before it gets moved here by mods :P Feb 20, 2012 at 22:00
• @seand: I tried this question formatted in a different way there a while back. No attention at all. Feb 21, 2012 at 0:11
• I saw there are LAPACK implementations for a number of languages. Maybe if you could provide an implementation of your code in a more mainstream language you'd get an answer... Though I must admit I haven't got a clue what a "skew Hermitian matrix" is :p Feb 21, 2012 at 16:46

Your code looks good. When s is small, there should be no need to use BLAS instead of matmul(). I like you exploit the shorthand array notation and the forall construct.

Debatable is your use of explicit shaped arrays. They could be very slightly faster, because they are contiguous, but in many cases the compiler will have to copy the arrays when calling your routines and the overall code will be slower. The preferred way in modern code are assumed shape arrays, e.g.

 complex(complex_kind), intent(in) :: H(:,:)
complex(complex_kind)             :: B(1:ubound(H,1),1:ubound(H,1))


but there may be reasons to use your way and it may be shorter.

Generally avoid using star notation real*4 *8 *16 for sizes of your variables. It is non-standard and obsolete. Use real(some_kind_constant).

You can use selected_real_kind() and selected_int_kind() (now preferred), to get the constants.

Or if you don't mind Fortran 2008 features and need to know the size in bits you can use kind constants from the iso_fortran_env module, like integer(int32), real(real64) and so on.

Another possibility from Fortran 2003 is to use kind constants from the iso_c_binding module, like integer(c_int) if you need interoperability with C code.

The solution is OK, if you know that you have to perform a time propagation of exactly $$\t\$$ with that same Hamiltonian and you have to do so countless times.

As someone else has suggested, ditch all the matmul and conjg and use only BLAS calls.

If you apply the propagator to a small number of vectors, however, computing the full propagator is not efficient, due to the matrix-matrix multiplication it entails. Instead, you should evaluate its action on the wavefunction you are propagating by means of three matrix-vector multiplications: $$e^{-i\mathbf{H}dt}\mathbf{c} = \mathbf{U} e^{-i \mathbf{E}dt}\mathbf{U}^\dagger \mathbf{c}$$ that is $$\mathbf{c}_1 = \mathbf{U}^\dagger \mathbf{c},\qquad \mathbf{c}_2 = e^{-i \mathbf{E}dt} \mathbf{c}_1,\qquad \mathbf{c} = \mathbf{U} \mathbf{c}_2$$

If you need to propagate multiple different times, of course, diagonalizing the Hamiltonian for every time is definitely not a good idea. Instead, you should do the diagonalization only once, and follow the triple-multiplication procedure above, whose computational cost scales as $$\N^2\$$ rather than $$\N^3\$$.

If your hamiltonian were to ever become time-dependent, $$\\mathbf{H}=\mathbf{H}(t)\$$, you must break the propagation in steps $$\mathbf{c}(t+dt) = \mathbf{U}(t+dt,t)\mathbf{c}(t),$$ where $$\\hat{U}(t+dt,t)\$$ is the time-ordered exponential propagator,

$$\mathbf{U}(t+dt,t)=\hat{T}e^{-i\int d\tau\mathbf{H}(\tau)}$$.

An excellent unitary second-order propagator is the mid-point exponential $$c(t+dt) = e^{-i \mathbf{H}(t+dt/2) dt} c(t).$$ To evaluate this step, there are many ways that are much better than a brutal direct diagonalization at each time step. First, if you can afford the diagonalization of the full Hamiltonian and if the Hamiltonian has the form $$\mathbf{H}(t) = \mathbf{H}_0 + F(t) \mathbf{H}_I$$ where both $$\\mathbf{H}_0\$$ and $$\\mathbf{H}_I\$$ are time-independent Hermitean operators and $$\F(t)\$$ is a real function, then the best approach is to split symmetrically the propagator: $$e^{-i \mathbf{H}(t+dt/2) dt} = e^{-i \mathbf{H}_0 dt/2}e^{-i \mathbf{H}_I F(t+dt/2) dt}e^{-i \mathbf{H}_0 dt/2} + o(dt^3),$$ where $$\ o(dt^3)\$$ is the difference between the two approximations to the time-step propagator. The mid-point propagator was already accurate only to second order anyway, so you are not losing any accuracy. With this arrangement, you need to perform the diagonalization only once, at the very beginning of the program, separately for $$\\mathbf{H}_0\$$ and $$\\mathbf{H}_I\$$, $$\mathbf{H}_0 = \mathbf{U}_0 \mathbf{E}^{(0)} \mathbf{U}_0^\dagger,\qquad\mathbf{H}_I = \mathbf{U}_I \boldsymbol{\Delta} \mathbf{U}_I^\dagger,$$ where $$\\mathbf{E}^{(0)}_{ij}=E^{(0)}_i\delta_{ij}\$$, $$\\boldsymbol{\Delta}_{ij}=\Delta_i\delta_{ij}\$$ are diagonal matrices. In this way, your propagator becomes $$U(t+dt,t) = \mathbf{U}_0 e^{-i \mathbf{E}^{(0)} dt/2}\mathbf{U}_0^\dagger\mathbf{U}_I e^{-i \boldsymbol{\Delta}_I F(t+dt/2) dt}\mathbf{U}_I^\dagger\mathbf{U}_0 e^{-i \mathbf{H}_0 dt/2}\mathbf{U}_0^\dagger + o(dt^3),$$ You can pre-compute the matrix $$\\mathbf{O}=\mathbf{U}_0^\dagger\mathbf{U}_I\$$, of course, in which case $$U(t+dt,t) = \mathbf{U}_0 e^{-i \mathbf{E}_0 dt/2}\mathbf{O} e^{-i \mathbf{E}_I F(t+dt/2) dt}\mathbf{O}^\dagger e^{-i \mathbf{H}_0 dt/2}\mathbf{U}_0^\dagger + o(dt^3).$$

If your time step is constant, you can neglect the first half free propagation, and merge the second with the first of the following step, thus obtaining $$U(t+dt,t) \simeq \mathbf{U}_0 e^{-i \mathbf{E}_0 dt}\mathbf{O} e^{-i \mathbf{E}_I F(t+dt/2) dt}\mathbf{U}_I^\dagger.$$ In this way, you have replaced a diagonalization (order $$\N^3\$$) with three consecutive matrix-vector multiplications (order $$\N^2\$$). If applicable to your case, therefore, this change will result in a wild speed up.

If either $$\H_0\$$ or $$\H_I\$$ are too large to be diagonalized (they are not equivalent, since $$\H_0\$$ is typically block diagonal and hence it is sufficient to diagonalize its diagonal blocks), then you can still use the symmetric splitting as a way to precondition the matrix, and use an iterative method (look up Krylov spaces) to determine the action of the operator on your target function, rather than the opearator over the whole Hilbert space, which is almost always an overkill.

PS: For folks not familiar with the simple syntax of modern Fortran, I suggest quick modern Fortran tutorial.