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