Solving Problems with zheev Subroutine in Do Loop

[Rositta]

I'm using zheev subroutine in a do loop (up to 2000 times) that contains other subroutines and functions. The matrix to diagonalize is only 5-by-5 however zheev gives the correct eigenvalues for only few cases even after reducing the number of times in the do loop. Any hint would be highly appreciated...

 

[Mark44]

I'm using zheev subroutine in a do loop (up to 2000 times) that contains other subroutines and functions. The matrix to diagonalize is only 5-by-5 however zheev gives the correct eigenvalues for only few cases even after reducing the number of times in the do loop. Any hint would be highly appreciated...

Without seeing your code, there's not we can do to help you diagnose the problem.

 

[Rositta]

Without seeing your code, there's not we can do to help you diagnose the problem.

Sorry.. Here is the main programme:

=>program xggauss
IMPLICIT NONE

! declare variables for part1
integer:: nval,i,k0
integer,parameter::kf=20

real*8,parameter::pi=3.141592653589
real*8:: x1,x2,dx,func,ss,x,r,del
complex*8::sphshift

! declare variables for part2
real::vb,bb
real*8,dimension(1:kf)::eb,k,kk,perturb
complex*8,dimension(1:kf)::delta
real*8::Eval(5)
integer,parameter::LDA=5
complex*16,dimension(LDA,LDA)::Evec
parameter(vb=1.,bb=1.)

external func,sphshift

print*,'lower limit='
read*,x1
print*,'upper limit='
read*,x2
print*,'number of intervals='
read*,nval

dx=(x2-x1)/nval
write(*,'(f6.3,1x,a,t12,a,t23,a/)') x1,'to','qgauss','expected'

do 12 k0=1,kf
k(k0)=k0*0.001
eb(k0)=k(k0)**2
kk(k0)=sqrt(vb+eb(k0))

do 11 i=1,nval
x=x1+i*dx
call ggauss(func,kk(k0),x1,x,ss)
11 continue

perturb(k0)=(2*bb/kk(k0)**2)*(1/(14*sqrt(pi)))*ss

delta(k0)=sphshift(k(k0),kk(k0))

call eigvals(perturb(k0),sphshift(k(k0),kk(k0)),Eval,Evec)
write(*,'(1x,I5,1x,f6.4,1x,f10.8,2x,f10.8,2x,f15.13,2x,f15.13,(2x,f15.13,1X,f20.18))')k0,k(k0),kk(k0),eb(k0),ss,perturb(k0),delta(k0)

CALL PRINT_RMATRIX('Eigenvalues', 1, 5, Eval, 1 )

print*,'********************************************************************************'
12 continue

end program xggauss

!
real*8 function func(x,r)
real*8:: x,r
func=x**4*((x*r)*((3/(x*r)**2-1)*(sin(x*r)/(x*r))-(3*cos(x*r)/(x*r)**2)))**2
end function func

!
complex*8 function sphshift(x,y)
real*8:: x,y,cb1,cb2,cb,db1,db2,db,del2,fb1,fb2

cb1=sin(x)*(4*x**2-9)+cos(x)*(9*x-x**3)
cb2=sin(x)*(9*x-x**3)-cos(x)*(4*x**2-9)

db1=sin(x)*(3-x**2)-3*x*cos(x)
db2=cos(x)*(3-x**2)+3*x*sin(x)

cb=db1*(sin(y)*(4*y**2-9)+cos(y)*(9*y-y**3))-cb1*(sin(y)*(3-y**2)-3*y*cos(y))
db=cb2*(sin(y)*(3-y**2)-3*y*cos(y))+db2*(sin(y)*(4*y**2-9)+cos(y)*(9*y-y**3))

del2=atan(-cb/db)
fb1=cos(del2)*sin(del2)/x
fb2=sin(del2)**2/x
sphshift=cmplx(fb1,fb2)

end function sphshift


and this is the subroutine I'm using to find the eigenvalues:


=>subroutine eigvals(perturb,sphshift,W,Q)
implicit none
real*8,intent(in)::perturb


!declare variable for zheev subroutine
integer::Info,Lwork
integer,parameter::m=5
integer,parameter::LDA=m
integer,parameter::Lwmax=100
complex*16,dimension(LDA,m),intent(out)::Q
complex*16,dimension(1:Lwmax)::Work
real*8,intent(inout)::W(m)
double precision::Rwork(3*m-2)
complex*8::sphshift

complex*16,dimension(m,m):ct,ones
data oct/(1.,0.),(0.,0.),(0.,0.),(0.,0.),(5.,0.),(0.,0.),(-4.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),&
(0.,0.),(6.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(-4.,0.),(0.,0.),(5.,0.),(0.,0.),(0.,0.),(0.,0.),(1.,0.)/,&
ones/(1.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(1.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),&
(0.,0.),(1.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(0.,0.),(1.,0.),(0.,0.),(0,0.),(0.,0.),(0.,0.),(0.,0.),(1.,0.)/

integer::i,j

intrinsic int,min

do i=1,m
do j=1,m
Q(i,j)=perturb*oct(i,j)+sphshift*ones(i,j)
end do
end do

CALL ZHEEV('V', 'L', m, Q, LDA, W, WORK, LWORK, RWORK, INFO )
Lwork=min(LWMAX,int(WORK(1)))
CALL ZHEEV('V', 'L', m, Q, LDA, W, WORK, LWORK, RWORK, INFO )

IF(INFO.GT.0) THEN
WRITE(*,*)'The algorithm failed to compute eigenvalues.'
STOP
END IF

end subroutine eigvals

 

[AlephZero]

I think you are trying to call ZHEEV to find the size of the work array, then call it again to solve the eigenproblem.

You need to set Lwork = -1 every time before the first call. Your code doesn't seem to set it to anything before the first call to ZHEEV. Most likely your code crashes "randomly" because Lwork eventually has a bad value.

 

[Rositta]

I think you are trying to call ZHEEV to find the size of the work array, then call it again to solve the eigenproblem.

You need to set Lwork = -1 every time before the first call. Your code doesn't seem to set it to anything before the first call to ZHEEV. Most likely your code crashes "randomly" because Lwork eventually has a bad value.

Thaks a lot AlephZero. All the eigenvalues are calculated correctly now :)