dangerous behavior of lapack routin zheev

Question

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.

Any relation to: https://stackoverflow.com/questions/51322226/inaccurate-zheev-eigen-values-and-vectors. See also the discussion about `kind` with https://stackoverflow.com/questions/51302912/incorrect-inconsistent-results-from-zgeev-lapack — albert, Jul 13 '18 at 12:56
I've related this to an answer with a different subroutine: it still comes down to having the correct arguments. — francescalus, Jul 13 '18 at 13:12

score 3 · Accepted Answer · answered Jul 13 '18 at 13:04

3

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

answered Jul 13 '18 at 13:04

albert

8,285
3
19
32

Shouldn't the compiler complain? – pawel_winzig Jul 14 '18 at 09:36
2

No the compiler has no knowledge of the interface unless you explicitly use and `interface` in your program or it is handled by means of a module (but they are, as far as I know, not provided by LaPack and are also compiler dependent. It is possible to generate interfaces yourself, but for the gfortran compiler I don't know the options, for the Intel compiler see the options `-gen_interfaces -warn interfaces -module`). – albert Jul 14 '18 at 09:45

dangerous behavior of lapack routin zheev

1 Answers1

Linked