I'm trying to write a program which computes the hamiltonian of a ising system. The first part of the program works, and i obtain the matrix that i want, but a problem occurs when i try to get the eigenvalues of this complex*16 matrix by using the subroutine zheev. As i read i have to call it 2 times, in order to assign the best possible value to the integer lwork. But for the given matrix, instead of getting the eigenvalues (-2.236, -1, 1, 2.236), as matlab does, i get (132, 0, 1, 0). Could somebody help me?
Thank you very much
module matrices
implicit none
type matrix
complex*16, dimension(:,:), allocatable :: el
integer :: row, col
end type
contains
type(matrix) function assign0(row, col)
integer :: ii, jj, row, col
assign0%row = row
assign0%col = col
allocate(assign0%el(assign0%row, assign0%col))
do ii = 1, assign0%row
do jj = 1, assign0%col
assign0%el(ii, jj) = cmplx(0d0, 0d0)
enddo
enddo
end function assign0
type(matrix) function assignid(row, col)
integer :: ii, jj, row, col
assignid%row = row
assignid%col = col
allocate(assignid%el(assignid%row, assignid%col))
do ii = 1, assignid%row
jj = ii
assignid%el(ii, jj) = cmplx(1d0, 0d0)
enddo
end function assignid
type(matrix) function assignsigx()
allocate(assignsigx%el(assignsigx%row, assignsigx%col))
assignsigx = assign0(2, 2)
assignsigx%el(1, 2) = (1d0, 0d0)
assignsigx%el(2, 1) = (1d0, 0d0)
end function
type(matrix) function assignsigz()
allocate(assignsigz%el(assignsigz%row, assignsigz%col))
assignsigz = assign0(2, 2)
assignsigz%el(1, 1) = (1d0, 0d0)
assignsigz%el(2, 2) = (-1d0, 0d0)
end function
type(matrix) function prodscal(lambda, ma1)
type(matrix) :: ma1
real :: lambda
integer :: ii, jj
prodscal%row = ma1%row
prodscal%col = ma1%col
allocate(prodscal%el(prodscal%row, prodscal%col))
do ii = 1, ma1%row
do jj = 1, ma1%col
prodscal%el(ii, jj) = lambda*ma1%el(ii, jj)
enddo
enddo
end function prodscal
type(matrix) function matsum(ma1, ma2, sign)
type(matrix) :: ma1, ma2
integer :: ii, jj, sign
matsum%row = ma1%row
matsum%col = ma1%col
allocate(matsum%el(matsum%row, matsum%col))
do ii = 1, ma1%row
do jj = 1, ma1%col
if (sign .gt. 0) then
matsum%el(ii, jj) = ma1%el(ii, jj) + ma2%el(ii, jj)
else
matsum%el(ii, jj) = ma1%el(ii, jj) - ma2%el(ii, jj)
endif
enddo
enddo
end function matsum
type(matrix) function tensprod(ma1, ma2)
type(matrix) :: ma1, ma2
integer :: ii, jj, kk, ll
tensprod%row = ma1%row*ma2%row
tensprod%col = ma1%col*ma2%col
allocate(tensprod%el(tensprod%row, tensprod%col))
tensprod = assign0(tensprod%row, tensprod%col)
do ii = 1, ma1%row
do jj = 1, ma1%col
do kk = 1, ma2%row
do ll = 1, ma2%col
tensprod%el((ii - 1)*ma2%row + kk, (jj - 1)*ma2%col + ll) = ma1%el(ii, jj)*ma2%el(kk, ll)
enddo
enddo
enddo
enddo
end function tensprod
end module matrices
program test
use matrices
implicit none
interface
function getdim()
integer :: getdim
end function getdim
end interface
integer :: row, col, ii, jj, N, ind, info, lwork, aus
type(matrix), dimension(2, 2) :: arrmat1, arrmat2
type(matrix), dimension(2) :: element1, element2, element3
type(matrix) :: elementtot, elementtot1, elementtot2
real :: lambda
double precision, dimension(:), allocatable :: w, rwork
complex*16, dimension(:), allocatable :: work
N = 2
allocate(w(15))
allocate(rwork(15))
allocate(work(15))
lwork = -1
aus = 2
lambda = 1d0
do ii = 1, N
ind = ii
do jj = 1, N
if(jj .eq. ind) then
arrmat1(ii, jj) = assignsigz()
else
arrmat1(ii, jj) = assignid(2, 2)
endif
enddo
enddo
do ii = 1, N
ind = ii
do jj = 1, N
if(jj .eq. ind) then
arrmat2(ii, jj) = assignsigx()
else
arrmat2(ii, jj) = assignid(2, 2)
endif
enddo
enddo
do ii = 1, N
element1(ii) = tensprod(arrmat1(ii, 1), arrmat1(ii, 2))
if (N .gt. 2) then
do jj = 3, N
element1(ii) = tensprod(element1(ii), arrmat1(ii, jj))
enddo
endif
enddo
do ii = 1, N
element2(ii) = tensprod(arrmat2(ii, 1), arrmat2(ii, 2))
if (N .gt. 2) then
do jj = 3, N
element2(ii) = tensprod(element2(ii), arrmat2(ii, jj))
enddo
endif
enddo
do ii = 1, N - 1
element3(ii) = assign0(2**N, 2**N)
element3(ii)%el = matmul(element2(ii + 1)%el, element2(ii)%el)
enddo
elementtot = assign0(2**N, 2**N)
elementtot1 = assign0(2**N, 2**N)
elementtot2 = assign0(2**N, 2**N)
do ii = 1, N
elementtot1 = matsum(elementtot1, element1(ii), 1)
enddo
do ii = 1, N - 1
elementtot2 = matsum(elementtot2, element3(ii), 1)
enddo
elementtot2 = prodscal(lambda, elementtot2)
elementtot = matsum(elementtot1, elementtot2, 1)
do ii = 1, 2**N
write(*,'(20G12.4)') elementtot%el(ii, 1:2**N)
enddo
call zheev('N', 'U', 2**N, elementtot%el, 2**N, w, work, lwork, rwork, info)
lwork = int(work(1))
call zheev('N', 'U', 2**N, elementtot%el, 2**N, w, work, lwork, rwork, info)
print*, '***', info
do ii = 1, 2**N
write(*,'(20G12.4)') work(ii)
enddo
stop
end program test
