Dear Python-users,
I am currently very confused about the Scipy routine to obtain the eigenvectors
of a complex matrix.In attached you find two files to diagonalize a 2X2 complex
Hermitian matrix, however, on my computer,
If I run python, I got:
[[ 0.80322132+0.j 0.59500941+0.02827207j] [-0.59500941+0.02827207j
0.80322132+0.j ]]
If I compile the fortran code, I got:
( -0.595009410289, -0.028272068905) ( 0.802316135182, 0.038122316497) (
-0.803221321796, 0.000000000000) ( -0.595680709955, 0.000000000000)
>From the scipy webpage, it is said that numpy.linalg.eig() provides nothing
>butan interface to lapack zheevd subroutine, which is used in my fortran code.
Would somebody be kind to tell me how to get consistent results?
Many thanks in advance.
Best wishes,
Hongbin
Ad hoc, ad loc and quid
pro quo
--- Jeremy Hilary Boob
!!*********************************************************************************************
!!
!!
!!*********************************************************************************************
!! ifort -CB *.f90 -L/usr/local/intel/lib/intel64 -lifcore -limf
-Wl,-rpath,/usr/local/intel/lib/intel64 -llapack_ifort -lblas_ifort -o 2X2.x
Program twobytwo
!!
!!
Implicit None
!!>>>>>>>>>>>>>>>>>> HAMILTONINA <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
Complex*16 :: Hamiltonian(2,2)
Complex*16 :: Hamilt(2,2)
Real*8,ALLOCATABLE :: eigenvals(:)
Real*8 :: eigentmp
Complex*16, ALLOCATABLE :: eigenvects(:,:),ceigenvals(:)
Complex*16 :: dummy(1,1)
Integer :: lwork, liwork, lrwork, info
Integer,PARAMETER :: lwmax=180000
Complex*16 :: work(lwmax)
Real*8,ALLOCATABLE :: rwork(:)
real*8 :: vl, vu
integer :: ne
Real*8 :: abstol
integer,allocatable :: ifail(:),iwork(:)
integer,allocatable :: order(:)
!!LOCAL
Integer :: i, j, k
Real*8, PARAMETER :: pi=3.1415926535897932384626d0
Complex*16,PARAMETER :: im_i=cmplx(0.0d0,1.0d0)
abstol=2.0d0*tiny(abstol)
allocate(eigenvals(2))
allocate(ceigenvals(2))
allocate(eigenvects(2,2))
lwork=12.0*2
lrwork=1+2*5+2*2**2
liwork=3+5*2
allocate(rwork(17*2))
allocate(iwork(15*2))
allocate(ifail(15*2))
allocate(order(2))
Hamiltonian=cmplx(0.0d0,0.0d0)
Hamiltonian(1,1)=cmplx(0.6d0,0.d0)
Hamiltonian(1,2)=cmplx(-1.97537668d0,-0.09386068d0)
Hamiltonian(2,1)=cmplx(-1.97537668d0,+0.09386068d0)
Hamiltonian(2,2)=cmplx(-0.6d0,0.0d0)
Hamilt=Hamiltonian
Call print_matrix("Hamiltonian",2, 2, Hamiltonian)
!!@_@
! lwork=-1
! work(:)=cmplx(0.0d0,0.0d0)
! rwork=0.0d0
!@_@ digonalizing with ZGEEV ---- omitting the overlap matrix
! call ZGEEV('No leftvectors','Vectors right', 2, HH1, 2,
ceigenvals,dummy,1,eigenvects, 2, work,lwork,rwork,info)
! lwork=min(lwmax,INT(work(1)))
! call ZGEEV('No leftvectors','Vectors right', 2, HH1, 2,
ceigenvals,dummy,1,eigenvects, 2, work,lwork,rwork,info)
! eigenvals(:)=REAL(ceigenvals)
!!@_@
! call
zheevx('V','A','U',2,Hamiltonian,2,vl,vu,1,2,abstol,ne,eigenvals,eigenvects,2,work,lwork,rwork,iwork,
ifail,info)
!!@_@
lwork=-1
liwork=-1
lrwork=-1
call zheevd('V','U',2,Hamiltonian,2,eigenvals, work, lwork, rwork,
lrwork, iwork, liwork, info)
lwork=MIN(lwmax,int(work(1)))
lrwork=min(lwmax,int(rwork(1)))
liwork=min(lwmax,iwork(1))
call zheevd('V','U',2,Hamiltonian,2,eigenvals, work, lwork, rwork,
lrwork, iwork, liwork, info)
eigenvects=Hamiltonian
IF (INFO .EQ. 0) then
!write(*,*) "Hamiltonian has been successfully diagonalized!"
ELSE
write(*,*) "Problems in digonalizing the
Hamiltonian................................."
STOP
End IF
!write(*,*) eigenvals
!call sorty(eigenvals,2,order)
write(*,*) eigenvals
call print_matrix("eigenvectors:",2,2,eigenvects)
!write(*, '(A, f8.4,1X,A,2i4)') "t = ", t, "order()=", order
Deallocate(eigenvals)
Deallocate(ceigenvals)
Deallocate(rwork)
Deallocate(iwork)
Deallocate(ifail)
Deallocate(eigenvects)
Deallocate(order)
End Program twobytwo
!!************************************************************************************************
Subroutine sort(eigenvals, n)
!!
!! Using the Gnome algorithm to sort n eigenvalues
!!---------------------------------------------------------------------------
!!
Integer :: n
Real*8 :: eigenvals(n)
!!LOCAL
Real*8 :: eigentmp
Integer :: i,j,k
! Write(*,*) 'Sorting >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>'
i=2
j=3
Do While (i <= n)
IF ( eigenvals(i-1) <= eigenvals (i) ) then
i=j
j=j+1
ELSE
eigentmp=eigenvals(i)
eigenvals(i)=eigenvals(i-1)
eigenvals(i-1)=eigentmp
i=i-1
IF ( i == 1 ) then
i=j
j=j+1
End IF
End IF
End do
End Subroutine sort
!!********************************************************************************
Subroutine sorty(eigenvals, n, order)
!!
!! Using the Gnome algorithm to sort n eigenvalues and berry curvature
!!---------------------------------------------------------------------------
!!
Integer :: n
Real*8 :: eigenvals(n)
Integer :: order(n)
!!LOCAL
Real*8 :: eigentmp
Integer :: i,j,k,o
Do i=1,n
order(i)=i
End Do
!Write(*,*) 'Sorting >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>'
i=2
j=3
Do While (i <= n)
IF ( eigenvals(i-1) <= eigenvals (i) ) then
i=j
j=j+1
ELSE
eigentmp=eigenvals(i)
eigenvals(i)=eigenvals(i-1)
eigenvals(i-1)=eigentmp
o=order(i)
order(i)=order(i-1)
order(i-1)=o
i=i-1
IF ( i == 1 ) then
i=j
j=j+1
End IF
End IF
End do
End Subroutine sorty
!!*****************************************************************************************
Subroutine PRINT_MATRIX(message, M, N, A)
!!
!! print out complex matrix
!!-----------------------------------------------------------------------------
!!
Character*(*) message
Integer M, N, I, J
Complex*16 A(M, N)
Write(*,*)
Write(*,*) message
Do I = 1, M
write(*,9998) ( A(I, J), J=1, N)
End Do
!
9998 FORMAT( 12(:,1X,'(',F16.12,',',F16.12,')'))
Return
End Subroutine PRINT_MATRIX
import numpy
H=[[0.6+0.0j, -1.97537668-0.09386068j],[-1.97537668+0.09386068j, -0.6+0.0j]]
eval,evecs=numpy.linalg.eig(H)
print eval
print evecs
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
