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...
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
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.