problems with zheev


by Rositta
Tags: zheev
Rositta
Rositta is offline
#1
Feb9-12, 08:02 AM
P: 13
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...
Phys.Org News Partner Science news on Phys.org
Better thermal-imaging lens from waste sulfur
Hackathon team's GoogolPlex gives Siri extra powers
Bright points in Sun's atmosphere mark patterns deep in its interior
Mark44
Mark44 is online now
#2
Feb9-12, 12:57 PM
Mentor
P: 20,982
Quote Quote by Rositta View Post
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
Rositta is offline
#3
Feb10-12, 06:34 AM
P: 13
Quote Quote by Mark44 View Post
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.1 3,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)::oct,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
AlephZero is offline
#4
Feb10-12, 02:51 PM
Engineering
Sci Advisor
HW Helper
Thanks
P: 6,344

problems with zheev


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
Rositta is offline
#5
Feb17-12, 04:53 AM
P: 13
Quote Quote by AlephZero View Post
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 :)


Register to reply

Related Discussions
What makes putnam-style problems 'different' from other problems? Academic Guidance 3
Projectile Motion problems (don't want answers to problems) Calculus & Beyond Homework 1
Chesmistry 12|VSEPER Problems|Quick Problems Biology, Chemistry & Other Homework 1
Can someone double check to see if I did these problems correctly? E-Flux problems Introductory Physics Homework 3