Problems with zheev

  1. 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...
     
  2. jcsd
  3. Mark44

    Staff: Mentor

    Without seeing your code, there's not we can do to help you diagnose the problem.
     
  4. 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)::eek: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
     
  5. AlephZero

    AlephZero 7,298
    Science Advisor
    Homework Helper

    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.
     
  6. Thaks a lot AlephZero. All the eigenvalues are calculated correctly now :)
     
Know someone interested in this topic? Share this thead via email, Google+, Twitter, or Facebook

Have something to add?