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
complex16 :: B(s,s), eigv(s,s), time_indep_schrodinger(s,s), work(2s-1)
real8 :: rwork(3s-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)

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 real4 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.

StackExchange Code Review Q#8917, answer score: 6