LAPACK/ScaLAPACK Development

Hi there,

I've noted some anormal behavior in my code when my matrix size turns to be between 41 and 47 regarding its diagonalisation. Due to the symmetric nature of my matrix I expect the spectrum of the eigenstate amplitudes (squared of each eigenstate component) to be symmetric as well. This has been working fine until now that I had to use a 41x41 matrix, giving as a result a non symmetric spectrum. I have tried to diagonalise the same matrix with MATLAB and it gives the expected symmetric result. After reviewing my code and read the corresponding lapack subroutines documentation I can't find any error, that's why I thought it could be a bug (?). I would really appreciate if you could either confirm or discard it.

I paste my code (due to I'm only allowed to attach 3 files):

program diagonal

implicit none

INTEGER :: i,j

INTEGER N
PARAMETER ( N = 41 )
INTEGER LDA
PARAMETER ( LDA = N )
INTEGER LWMAX
PARAMETER ( LWMAX = 1394 )

INTEGER INFO, LWORK

DOUBLE PRECISION A( LDA, N ), W( N ), WORK(LWMAX)


!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!open hamiltonian matrix!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!
open(unit=1,file="hamiltoniano.matrix")
read(1,*) ((A(i,j),i=1,N),j=1,N)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!optimize workspace size (lwork)!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
LWORK=-1

call dsyev('V','U',N,A,LDA,W,WORK,LWORK,INFO)
if(INFO/=0) stop 'Error in dsyev'
LWORK=MIN(LWMAX,INT(WORK(1)))
print*, LWORK

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!diagonalisation at the opt lwork!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
call dsyev('V','U',N,A,LDA,W,WORK,LWORK,INFO)
if(INFO/=0) stop 'Error in dsyev'

!!!!!!!!!!!!!!!!!!!!!!
!!store eigenvectors!!
!!!!!!!!!!!!!!!!!!!!!!
open (unit=60,file='eigen.matrix',status='unknown')
do i=1,N
write(60,*) i,(abs(A(i,j)*A(i,j)),j=1,N)
enddo
close(60)

end program

Thank you very much,
Marta P. Estarellas


p.d. :I also attach my matrix, the results obtained with MATLAB and the results obtained with dsyev (I also tried the complex analogous zheev and gives the same problem)

Try to use dsyevd. I understand that your matrix has some structure ( eg Hamiltonian) and that, based on this structure, you expect the eigenvector matrix to have some properties (eg symmetry). If eigenvalues are multiple, two general purpose symmetric eigensolvers are likely to return different eigenvecto matrices. One that you may like. Maybe you will not be happy with any. Matlab is using dsyevd. So maybe try using dsyevd and this might be better for you. In general we will not be able to guarantee you to find back the property you expect in particular with multiplicities. To get properties in eigenvector matrix, either used a specialized eigensolver ( see Daniel Kressner PhD thesis for example ) or try to post process LAPACK output.