LAPACK/ScaLAPACK Development

Hello,

When using zheev on a hermitian matrix in fortran, when some components of the input matrix are imaginary, incorrect eigenvectors are returned, however the correct real eigenvalues are found. This is highlighted by the diagonalization of a simple matrix;

0 -i 0
A = i 0 i
0 i 0

Which should yield the eigenvalues; {-sqrt(2), 0, sqrt(2)} with the corresponding eigenvectors; {(-1/2, i*sqrt(2)/2, 1/2), (sqrt(2)/2, 0, sqrt(2)/2), (-1/2, -i*sqrt(2)/2, 1/2)}.

The output of zheev however, is; {-sqrt(2), 0, sqrt(2)} with the corresponding eigenvectors; {(-1/2, sqrt(2)/2, -1/2), (i*sqrt(2)/2, 0, -i*sqrt(2)/2), (1/2, sqrt(2)/2, 1/2)}

The code is printed below - this was compiled and tested using either pgi or ifort with either the pgi lapack and blas or the intel mkl, respectively.

Really like to know what is going on here....

Cheers!
Nick




program eigenvector_test
implicit none

real(kind=8),allocatable::rwork(:),Evalues(:)
complex(kind=8),allocatable::Evectors(:,:),work(:)
integer::N,LWORK,diaginfo,i

LWORK = 10
N = 3

allocate(Evalues(N),rwork(5*N),Evectors(N,N),work(LWORK))

LWORK = -1
call zheev('V','U',N,Evectors,N,Evalues,work,LWORK,rwork,diaginfo)
LWORK = INT(work(1))
deallocate(work)
allocate(work(LWORK))
if(diaginfo /= 0) write(6,*) "DIAG FAIL"

Evectors(1,1) = (0.0,0.0)
Evectors(1,2) = (0.0,-1.0)
Evectors(1,3) = (0.0,0.0)
Evectors(2,1) = (0.0,1.0)
Evectors(2,2) = (0.0,0.0)
Evectors(2,3) = (0.0,-1.0)
Evectors(3,1) = (0.0,0.0)
Evectors(3,2) = (0.0,1.0)
Evectors(3,3) = (0.0,0.0)

do i=1,N
write(6,*) Evectors(i,:)
write(6,*) ""
end do

Evalues(:) = 0.0

call zheev('V','U',N,Evectors,N,Evalues,work,LWORK,rwork,diaginfo)
if(diaginfo /= 0) write(6,*) "DIAG FAIL"
do i=1,N
write(6,*) i, Evalues(i)
write(6,*) Evectors(i,:)
write(6,*) ""
end do

deallocate(work,Evalues,rwork,Evectors)

end program eigenvector_test

To make sure you have correct eigenvectors check || V^H * A - D * V^H || / || A ||

A couple other aspects to keep in mind when comparing results from different software routines:
(i) an eigenvector for a particular eigenvalue doesn't necessarily have to be unique, and
(ii) some routines normalize eigenvectors while others don't (this tripped me up once; I was comparing results from different programs and couldn't figure out why the eigenvalues were the same but the eigenvectors were different--until I realized one program was normalizing the eigenvectors output).

For more info, you can take a look at http://icl.cs.utk.edu/lapack-forum/viewtopic.php?t=298

Julie

My mistake, I was reading the Eigenvector matrix as row vectors - although, I'm sure I tried the output as column vectors too...

Anyhow - no errors, all working well.

Cheers!
Nick