为何ZHEEV函数对同一矩阵输出不同特征值？Fortran代码问题

先说明：你用的ZHEEV是LAPACK里专门给厄米矩阵算特征值的函数，你的输入矩阵是标准厄米矩阵（ha(2,1)是ha(1,2)的共轭转置），理论上肯定能算出正确结果。循环外调用出问题，十有八九是内存初始化、变量作用域或者参数传递的锅，给你列几个排查的关键点：

1. 查变量的作用域和初始化 #

• 循环外的ha矩阵、存特征值的w数组、工作区work这些变量，是不是声明后没正确初始化就传给ZHEEV了？Fortran里没初始化的局部变量是随机值，但全局变量可能默认是0，要是之前有其他操作污染了内存，就会算错。
• 对比循环内的代码：是不是循环里每次都重新初始化ha？而循环外的ha可能被之前的调用改了——ZHEEV就算只算特征值（JOBZ='N'），也会把输入矩阵覆盖成中间结果，要是你循环外调用一次后没重新赋值ha就再调用，那第二次算的是被改烂的矩阵，肯定出0。

2. 核对ZHEEV的参数传递 #

ZHEEV的Fortran标准调用格式是：

CALL ZHEEV(JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, INFO)

重点盯这几个参数：

• UPLO：选'U'就只读上三角，选'L'读下三角，你的矩阵上下三角是共轭的，选哪个都没问题，但要是传错字符（比如打错成'u'或者其他），函数读错元素就会算错。
• LDA：必须大于等于矩阵阶数N（这里N=2），要是传成1或者其他数，会内存越界读错数据。
• WORK和LWORK：工作区数组大小够不够？要是LWORK设太小，ZHEEV根本没法正常算。可以先调用一次让它返回需要的最小大小，再分配工作区：

  INTEGER :: lwork, info
  DOUBLE COMPLEX :: work(1)
  CALL ZHEEV('N', 'U', 2, ha, 2, w, work, -1, rwork, info)
  lwork = INT(work(1))
  DEALLOCATE(work)
  ALLOCATE(work(lwork))
  CALL ZHEEV('N', 'U', 2, ha, 2, w, work, lwork, rwork, info)
  

  说不定循环内刚好工作区大小够，循环外的LWORK设错了才失败。

3. 检查数据类型对不对 #

所有变量类型必须匹配：

• ha得是DOUBLE COMPLEX（对应LAPACK的Z前缀函数），要是用了单精度的COMPLEX，类型不匹配直接算错。
• w和rwork必须是DOUBLE PRECISION类型，rwork的大小至少是2*N-2（这里N=2，所以至少2个元素）。

给你个正确的示例代码 #

你可以对照自己的代码改：

PROGRAM test_zheev
  IMPLICIT NONE
  INTEGER, PARAMETER :: N = 2
  DOUBLE COMPLEX :: ha(N,N)
  DOUBLE PRECISION :: w(N)
  DOUBLE COMPLEX :: work(1)
  DOUBLE PRECISION :: rwork(2*N-2)
  INTEGER :: lwork, info

  ! 初始化哈密顿矩阵
  ha(1,1) = DCMPLX(0.0d0, 0.0d0)
  ha(1,2) = DCMPLX(0.0d0, -4.8497422611928558d0)
  ha(2,1) = DCMPLX(0.0d0, 4.8497422611928558d0)
  ha(2,2) = DCMPLX(0.0d0, 0.0d0)

  ! 先查需要多大的工作区
  CALL ZHEEV('N', 'U', N, ha, N, w, work, -1, rwork, info)
  lwork = INT(work(1))
  DEALLOCATE(work)
  ALLOCATE(work(lwork))

  ! 正式调用算特征值
  CALL ZHEEV('N', 'U', N, ha, N, w, work, lwork, rwork, info)

  ! 输出结果
  IF (info == 0) THEN
    PRINT *, '特征值：'
    PRINT *, w(1), w(2)
  ELSE
    PRINT *, 'ZHEEV调用失败，INFO=', info
  END IF

END PROGRAM test_zheev

总结 #

优先查这几点：

• 循环外的ha是不是在调用ZHEEV前被改了（比如之前的调用覆盖了它）
• 工作区WORK的大小够不够
• UPLO、LDA这些参数有没有传错

内容的提问来源于stack exchange，提问作者Kieran