subroutine deig(A,Eval,Evec,jobvl,jobvr)
!IN matrix - OUT eigenvectors: A.evec = eval.evec ; evec(i,j) = ith component of jth vec.
real(8),intent(in) :: A(:,:)
complex(8),intent(out) :: Eval(:)
complex(8),intent(out) :: Evec(:,:)
character(len=1),optional :: jobvl,jobvr
character(len=1) :: jobvl_,jobvr_
real(8),dimension(:,:),allocatable :: At,vl,vr
real(8),dimension(:),allocatable :: wi,wr
real(8),dimension(:),allocatable :: work
real(8),dimension(1) :: lwork_guess
integer :: info,lda,ldvl,ldvr,lwork,n,i
jobvl_='N';if(present(jobvl))jobvl_=jobvl
jobvr_='V';if(present(jobvr))jobvr_=jobvr
if(jobvl_=='V')stop "deig error: jobvl = V is not supported yet."
lda = size(A(:,1))
n = size(A(1,:))
call assert_shape(A,[n,n],"solve","A")
call assert_shape(Evec,[n,n],"solve","Evec")
ldvl = n
ldvr = n
allocate(At(n,n),wr(n),wi(n),vl(ldvl,n),vr(ldvr,n))
!Copy the input Matrix
At = A
!1st Call: Query the right size for the working array.
call dgeev(jobvl_, jobvr_, n, At, lda, wr, wi, vl, ldvl, vr, ldvr, lwork_guess, -1, info)
if(info /= 0) then
print*, "dgeev returned info = ",info
if (info < 0) then
print*, "the",-info,"-th argument had an illegal value"
else
print*, "the QR algorithm failed to compute all the"
print*, "eigenvalues, and no eigenvectors have been computed;"
print*, "elements ", info+1, ":", n, "of WR and WI contain eigenvalues which"
print*, "have converged."
end if
stop 'deig error: 1st call dgeev'
endif
lwork = int(lwork_guess(1)) ;allocate(work(lwork))
!2nd Call: Actual solution of the eigenproblem
call dgeev(jobvl_, jobvr_, n, At, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info)
if(info /= 0) then
print *, "dgeev returned info = ",info
if (info < 0) then
print*, "the",-info,"-th argument had an illegal value"
else
print*, "the QR algorithm failed to compute all the"
print*, "eigenvalues, and no eigenvectors have been computed;"
print*, "elements ", info+1, ":", n, "of WR and WI contain eigenvalues which"
print*, "have converged."
end if
stop 'deig error: 2nd call dgeev'
end if
Eval = wr + xi*wi
! as DGEEV has a rather complicated way of returning the eigenvectors,
! it is necessary to build the complex array of eigenvectors from
! two real arrays:
do i = 1,n
if(wi(i) > 0d0) then ! first of two conjugate eigenvalues
Evec(:, i) = vr(:, i) + xi*vr(:, i+1)
elseif(wi(i) < 0d0) then ! second of two conjugate eigenvalues
Evec(:, i) = vr(:, i-1) - xi*vr(:, i)
else
Evec(:, i) = vr(:, i)
end if
end do
end subroutine deig