lapack routin zheev 的危险行为

Question

我在使用 lapack 例程时偶然发现了一个奇怪的行为 zheev()。有两个问题我不明白

1) 我的一个全局变量似乎被 zheev() 覆盖了。下面的小程序展示了它：

[用 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

输出为

5.9999999999999995E-004
   0.0000000000000000     
   0.0000000000000000     
   0.0000000000000000     
   0.0000000000000000     
   0.0000000000000000     
   0.0000000000000000     
   0.0000000000000000

这可以通过将 stepD 设置为参数来解决。但我不明白这种行为。越来越奇怪了:

2) 将定义中的 lwork 设置为某个值，如

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

（只需将上面对应的行改一下即可）给出以下结果：

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

怎么会这样？将 lwork 设置为 10 应该没有任何效果，因为稍后它会设置为 -1。一个重要的概念是：如果删除

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

代码运行良好。

Answer 1

来自 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 
    )

这里我们看到 RWORK 是一个双精度类型的数组，但在提供的代码中 rwork 是一个标量整数

lapack routin zheev 的危险行为

dangerous behavior of lapack routin zheev

fortran

matrix

eigenvalue

lapack