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? – qubyte Feb 12 '12 at 12:53
• Turns out that there is, but fortran is not supported. Turned off syntax highlighting for the code above. – qubyte Feb 13 '12 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 – seand Feb 20 '12 at 22:00
• @seand: I tried this question formatted in a different way there a while back. No attention at all. – qubyte Feb 21 '12 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 – dvdvorle Feb 21 '12 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.