I recently asked a question along the same lines for skew-Hermitian matrices. Inspired by the success of that question, and after banging my head against a wall for a couple of hours, I'm looking at the matrix exponential of real asymmetric matrices. The route to finding the eigenvalues and eigenvectors seems rather convoluted, and I'm afraid I've gotten lost.
Background: Some time ago I asked this question on the theoretical physics SE. The result allows me to phrase master equations as real asymmetric matrices. In the time-independent case, the master equation is solved by exponentiating this matrix. In the time-dependent case it will require integration. I'm only concerned with time-independence at the moment.
Upon looking at the various subroutines I think I should be calling (?gehrd, ?orghr, ?hseqr...) it is unclear if it would be simpler to cast the matrix from real*8
to complex*16
and proceed with the complex double versions of these routines, or stick with real*8
and take the hit of doubling the number of my arrays and making a complex matrix of them later.
So, which routines should I be calling (and in what order), and should I use the real double versions or the complex double versions? Below is an attempt at doing this with real double versions. I've become stuck finding the eigenvalues and eigenvectors of L*t
.
function time_indep_master(s,L,t)
! s is the length of a side of L, which is square.
! L is a real*8, asymmetric square matrix.
! t is a real*8 value corresponding to time.
! This function (will) compute expm(L*t).
integer, intent(in) :: s
real*8, intent(in) :: L(s,s), t
real*8 :: tau(s-1), work(s), wr(s), wi(s), vl
real*8, dimension(s,s) :: time_indep_master, A, H, vr
integer :: info, m, ifaill(2*s), ifailr(2*s)
logical :: sel(s)
A = L*t
sel = .true.
call dgehrd(s,1,s,A,s,tau,work,s,info)
H = A
call dorghr(s,1,s,A,s,tau,work,s,info)
call dhseqr('e','v',s,1,s,H,s,wr,wi,A,s,work,s,info)
call dhsein('r','q','n',sel,H,s,wr,wi,vl,1,vr,s,2*s,m,work,ifaill,ifailr,info)
! Confused now...
end function