Register to reply 
Problems with zheevby Rositta
Tags: zheev 
Share this thread: 
#1
Feb912, 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 5by5 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...



#2
Feb912, 12:57 PM

Mentor
P: 21,215




#3
Feb1012, 06:34 AM

P: 13

=>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=(x2x1)/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)**21)*(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**29)+cos(x)*(9*xx**3) cb2=sin(x)*(9*xx**3)cos(x)*(4*x**29) db1=sin(x)*(3x**2)3*x*cos(x) db2=cos(x)*(3x**2)+3*x*sin(x) cb=db1*(sin(y)*(4*y**29)+cos(y)*(9*yy**3))cb1*(sin(y)*(3y**2)3*y*cos(y)) db=cb2*(sin(y)*(3y**2)3*y*cos(y))+db2*(sin(y)*(4*y**29)+cos(y)*(9*yy**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*m2) 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 


#4
Feb1012, 02:51 PM

Engineering
Sci Advisor
HW Helper
Thanks
P: 6,953

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. 


#5
Feb1712, 04:53 AM

P: 13




Register to reply 
Related Discussions  
What makes putnamstyle problems 'different' from other problems?  Academic Guidance  3  
Projectile Motion problems (don't want answers to problems)  Calculus & Beyond Homework  1  
Chesmistry 12VSEPER ProblemsQuick Problems  Biology, Chemistry & Other Homework  1  
Can someone double check to see if I did these problems correctly? EFlux problems  Introductory Physics Homework  3 