!Define a random symmetric matrix, compute eigenvalues and eigenvectors using DSYEV:
!
! Query DSYEV for optimal LWORK
!
!-----------------------------------------------------------------------------------
program     test_evs
 implicit none
 integer             :: N
 real(8),allocatable :: A(:,:),B(:,:),l(:),V(:),WORK(:)
 real(8)             :: QLWORK(1)
 integer             :: i,j
 integer             :: LDA,LWORK,INFO
 character(1)        :: JOBZ='V',UPLO='U' 

 print *,'# Enter N'
 read  *,         N
 print *,'# ',    N
 print *,'# LWORK min = 3*N-1 = ',3 * N -1
 ALLOCATE(A(N,N),B(N,N),l(N),V(N))

 !compute LWORK: a call to DSYEV with LWORK=-1 returns optimal LWORK in WORK(1)
 LWORK = -1
 call dsyev(JOBZ,UPLO,N,A,N,l,QLWORK,LWORK,INFO)
 LWORK=QLWORK(1)
 ALLOCATE  (WORK(LWORK))

 print *,'# LWORK     = ',QLWORK(1),LWORK 

 call random_number(A)
 A = A + TRANSPOSE (A)

 B=A !save A

 print *,'------------------- A(i,j) ----------------------------'
 do i = 1,N
  print '(1000F8.4)',A(i,:)
 end do

!Compute eigenvalues and eigenvectors
 call DSYEV(JOBZ,UPLO,N,A,N,l,WORK,LWORK,INFO)
 if(INFO /= 0 ) stop 'DSYEV: INFO /= 0'

 print *,'-------------------------------------------------------'
 do i = 1, N
  print *,'-------------------- i= ',i,'------------------------------'
  print *,' l_i= ',l(i)
  V = A(:,i)
  print *,' V_i= ',V
  print *,' ||A.V_i - l_i V_i||= ',SUM(ABS( MATMUL(B,V)-l(i)*V ))/N
 end do

!Check Orthonormality:
 do i = 1, N
  print *,'-------------------- i= ',i,'------------------------------'
  print *,'|1-|v_i|^2|= ', ABS(1.0D0-DOT_PRODUCT(A(:,i),A(:,i)))
  print *,'|v_i . v_j|= ',(ABS(      DOT_PRODUCT(A(:,i),A(:,j))),j=1,N)
 end do

end program test_evs