dangerous behavior of lapack routin zheev [duplica

Question

This question already has an answer here:

• Incorrect/inconsistent results from zgeev() LAPACK 2 answers

I stumbled upon an odd behavior when using the lapack routine zheev(). There are two issues which I do not understand

1) One of my global variables seems to be overwritten by zheev(). The following small program shows it:

[compiled with gfortran -o test test.f90 -llapack -lblas]

program test    
implicit none

integer, parameter :: dp = 8
integer, parameter :: dim = 3

integer :: l

real(dp), parameter :: kmin = -0.03_dp
real(dp), parameter :: kmax = 0.03_dp
integer, parameter  :: steps = 100
real(dp) :: stepD = (kmax - kmin)/real(steps)
complex(dp) :: matrix(dim,dim)=0.

integer :: info
integer :: rwork=3*dim-2, lwork, lwmax=100
real(dp) :: evals(dim) 
complex(dp) :: work(3*dim-2)

lwork=-1
call zheev('N','U',size(matrix,1), matrix, dim, evals, work, lwork, &
            rwork, info) 
lwork = min( lwmax, int( work(1) ))

do l = 0, 3
    write(*,*) stepD
    call zheev('N', 'U', size(matrix,1), matrix, dim, evals, work, &
    lwork, rwork, info) 
    write(*,*) stepD
end do

end program test

The output is

5.9999999999999995E-004
   0.0000000000000000     
   0.0000000000000000     
   0.0000000000000000     
   0.0000000000000000     
   0.0000000000000000     
   0.0000000000000000     
   0.0000000000000000  

This can be cured by setting stepD to a parameter. But I do not understand this behavior. It is getting even more weird:

2) Setting lwork in the definition to some value like

integer :: rwork=3*dim-2, lwork=10, lwmax=100

(simply change the corresponding line above) gives the following result:

   5.9999999999999995E-004
   5.9999991208314896E-004
   5.9999991208314896E-004
   5.9999991208314896E-004
   5.9999991208314896E-004
   5.9999991208314896E-004
   5.9999991208314896E-004
   5.9999991208314896E-004

How can this happen? Setting lwork to 10 should have no effect since later on it is set to -1. An important notion is: If one removes

lwork=-1
call zheev('N','U',size(matrix,1), matrix, dim, evals, work, lwork, &
            rwork, info) 
lwork = min( lwmax, int( work(1) ))

the code works fine.

Answer 1

From the documentation of zheev (http://www.netlib.org/lapack/explore-html/df/d9a/group__complex16_h_eeigen_gaf23fb5b3ae38072ef4890ba43d5cfea2.html#gaf23fb5b3ae38072ef4890ba43d5cfea2):

subroutine zheev    (   character   JOBZ,
        character   UPLO,
        integer     N,
        complex*16, dimension( lda, * )     A,
        integer     LDA,
        double precision, dimension( * )    W,
        complex*16, dimension( * )      WORK,
        integer     LWORK,
        double precision, dimension( * )    RWORK,
        integer     INFO 
    )   

here we see that RWORK is an array of type double precision, but in the supplied code rwork is a scalar integer