拉帕克·鲁廷·哲耶夫的危险行为



我在使用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  

这可以通过将步骤D设置为参数来解决。但我不理解这种行为。它变得更加奇怪:

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

代码运行良好。

来自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是一个标量整数

最新更新