As said in the title, I'm new in fortran and lapack as well, so I hope that my question will not be too trivial.
Anyway, I'm banging my head on that particular problem since two days and google was not helpful this time therefore I'm here.
The code which I have problem with is the following (cleaned a bit for readability):
- Code: Select all
program hubbard_1D
implicit none
!variables
double precision :: n,mz,U,t
!real :: y
integer :: site,Hsize,i,j
double precision, dimension (:,:), allocatable :: Hamilt,Heigvec
double precision, dimension (:), allocatable :: Hdiag,Hoffdiag,Heigval
!program
write (*,*) "********************"
write (*,*) "* 1D Hubbard Model *"
write (*,*) "********************"
write (*,*)
n = 4.
mz = 2.
U = 2.
t = 1.
!Hamiltonian Matrix creation
allocate ( Hamilt ( Hsize , Hsize ) )
allocate ( Hdiag ( Hsize ) )
allocate ( Hoffdiag ( Hsize - 1 ) )
do i = 1,Hsize
do j = 1,Hsize
if (i == j) then
Hamilt(i,j) = U*(n+(-1)**i*mz)/2
else if ((i == j+1) .or. (i == j-1)) then
Hamilt(i,j) = -t
else
Hamilt(i,j) = 0.
end if
end do
end do
do i = 1,Hsize
Hdiag(i) = Hamilt(i,i)
end do
do i = 1,(Hsize - 1)
Hoffdiag(i) = Hamilt(i,i+1)
end do
!Call Lapack routine dstev : eigenvect+eigenval of tridiagonal matrix
call diatri(Hdiag,Hoffdiag,Heigval,Heigvec,Hsize)
write (*,*) "Eigenvalue list:",Heigval
deallocate ( Hamilt )
deallocate ( Hdiag )
deallocate ( Hoffdiag )
end program hubbard_1D
!---------------------------------------------------------!
!Diagonalization of a real symmetric tridiagonal matrix !
!---------------------------------------------------------!
!input: a(n) = diagonal matrix elements !
! b(n) = off-diagonal matrix elements !
! n = size of the matrix !
!output: eig(n) = eigenvalues in ascending order !
! vec(n,n) = eigenvectors !
!---------------------------------------------------------!
!Uses the LAPACK subroutine DSTEV !
!---------------------------------------------------------!
subroutine diatri(a,b,eig,vec,n)
!-------------------------------!
implicit none
integer :: n,inf
real(8) :: a(n),b(n),d(n),e(n),eig(n),vec(n,n),work(max(1,2*n-2))
d=a
e=b
call dstev('V',n,d,e,vec,n,work,inf)
eig=d
end subroutine diatri
!---------------------!
As you can see, I stole the subroutine calling dstev on the net, assuming that it will work better that the one I first wrote (not included in the previous lines).
Apparently it is not the problem.
I compiled this code on at least two different machines (mac OSX using -framework vecLib OR direct -llapack - lblas option ; and a small Unix cluster) using two different compilers (g95 or gfortran) and I don't get errors at compilation time.
BUT
when I execute the program, I get a segfault. I traced it to the call of dstev and...I'm stuck ! No idea what the problem is !
Any help will be warmly welcomed.
Thanks,
Schrall
