zgeev giving eigenvectors which are not orthogonal

I try to diagonalize a matrix using zgeev and it giving correct eigenvalues but the eigenvectors are not orthogonal.

program complex_diagonalization 
implicit none
integer,parameter :: N=3
integer::i,j
integer,parameter :: LDA=N,LDVL=N,LDVR=N
real(kind=8),parameter::q=dsqrt(2.0d0),q1=1.0d0/q
integer,parameter :: LWMAX=1000
integer :: INFO,LWORK
real(kind=8) :: RWORK(2*N)
complex(kind=8) :: B(LDA,N),VL(LDVL,N),VR(LDVR,N),W(N),WORK(LWMAX)
external::zgeev
!matrix defining
B(1,1)=0.0d0;B(1,2)=-q1;B(1,3)=-q1
B(2,1)=-q1;B(2,2)=0.50d0;B(2,3)=-0.50d0
B(3,1)=-q1;B(3,2)=-0.5d0;B(3,3)=0.50d0  
LWORK=-1
 CALL ZGEEV('Vectors','Vectors',N,B,LDA,W,VL,LDVL,VR,LDVR,WORK,LWORK,RWORK,INFO)
LWORK=MIN(LWMAX,INT(WORK(1)))
CALL ZGEEV('Vectors','Vectors',N,B,LDA,W,VL,LDVL,VR,LDVR,WORK,LWORK,RWORK,INFO)

IF( INFO.GT.0 ) THEN
 WRITE(*,*)'The algorithm failed to compute eigenvalues.'
 STOP
END IF
!eigenvalues
do i=1,N
WRITE(*,*)W(i)
enddo

!eigenvectors
do i=1,N
WRITE(*,*)(VR(i,j),j=1,N)
ENDDO

end

and the result I am getting are this: eigenvalues:

( 0.99999999999999978,0.0000000000000000)
(-0.99999999999999978,0.0000000000000000)
( 0.99999999999999978,0.0000000000000000)

eigenvectors

 (0.70710678118654746,0.0000000000000000)
 (-0.50000000000000000,0.0000000000000000)
 (-0.50000000000000000,0.0000000000000000)


 (0.70710678118654746,0.0000000000000000)
(0.50000000000000000,0.0000000000000000)
 (0.50000000000000000,0.0000000000000000)

(-0.11982367636731203,0.0000000000000000)
( 0.78160853028734012,0.0000000000000000)
(-0.61215226207528295,0.0000000000000000)

you can see that the third eigenvector is not orthogonal with one of the two eigenvectors. What I am expecting is that in the third eigenvector first entry should be zero and second entry will be minus of third entry and because it's a unit vector it will be 0.707.

Solution

A real symmetric matrix has three orthogonal eigenvectors if the three eigenvalues are unique. Only the eigenvectors corresponding to distinct eigenvalues have tobe orthogonal. https://math.stackexchange.com/a/1368948/134138

The Hermitian specialized routine ZHEEV should guarantee orthogonality of the eigenvectors as suggested by Ian Bush. Or in your case you can also consider DSYEV (because your matrix is real).

The situation is well described in this post from LAPACK Forum http://icl.cs.utk.edu/lapack-forum/archives/lapack/msg01352.html

From the documentation:

DSYEV:  
*          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
*          orthonormal eigenvectors of the matrix A.

ZHEEV:
*          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
*          orthonormal eigenvectors of the matrix A.