厄密矩阵对角化中的特征向量不连续性

Question

我需要对一个 2x2 Hermitian 矩阵进行对角化，该矩阵取决于连续变化的参数 x。对于对角化，我使用 EISPACK。当我将特征向量的实部和虚部绘制为 x 的函数时，我注意到它们存在不连续性。特征值计算正常。当我在 Maxima 中绘制特征向量时，解似乎是连续的。我需要连续的特征向量，因为下一步我需要计算它们的导数。

在我用作测试的 f77 代码下方（在 mingw 上使用 gfortran 编译）。

program Eigenvalue

  implicit none

  integer n, m 
  parameter (n=2)
  integer  ierr, matz, i, j 
  double precision x, dx, xf, amp, xin 
  double precision w(n)
  double precision Ar(n,n), Ai(n,n)
  double precision xr(n,n), xi(n,n)
  double precision fm1(2,n)   ! f77
  double precision fv1(n)     ! f77 
  double precision fv2(n)     ! f77       
  double complex psi1a, psi1b, psi2a, psi2b 

  m = 51              

  xf  = 10.d0  
  xin = 0.0d0     

  amp = 2.d0 
  dx = (xf - xin)/(m-1)

  do i = 1, m     
  x = dx*(i-m) + xf   

  Ar(1,1) = dsin(x)**2
  Ar(1,2) = amp*dcos(x)
  Ar(2,1) = amp*dcos(x)
  Ar(2,2) = dcos(x)**2

  Ai(1,1) = 0.0d0
  Ai(1,2) = amp*dsin(x)  
  Ai(2,1) = -amp*dsin(x) 
  Ai(2,2) = 0.0d0

  matz = 1

  call ch ( n, n, ar, ai, w, matz, xr, xi, fv1, fv2, fm1, ierr ) !f77

  write(20,*) x, w(1), w(2)

  write(21,*) x, xr(1,1), xi(1,1)
  write(22,*) x, xr(2,1), xi(2,1)

  write(23,*) x, xr(1,2), xi(1,2)
  write(24,*) x, xr(2,2), xi(2,2)       

!     autovetor 1

  psi1a = cmplx(xr(1,1),xi(1,1)) 
  psi1b = cmplx(xr(1,2),xi(1,2)) 

!     autovetor 2 

  psi2a = cmplx(xr(2,1),xi(2,1)) 
  psi2b = cmplx(xr(2,2),xi(2,2))         

  end do
  end

Answer 1

虽然不是真正的答案，但以下是我在 LAPACK 中使用的代码。我使用了最新版本的 LAPACK 和 BLAS，具有以下编译器选项：

gfortran -Og -std=f2008 -Wall -Wextra {location_of_lapack}/liblapack.a {location_of_blas}/blas_LINUX.a main.f90 -o main

我正在 Mac OS X 上使用来自自制软件的 gfortran 6.3.0 进行编译。

正如 Ian 上面提到的，dcos 之类的东西被 cos 替换了，我使用了 KIND= 公式来确保相同的精度。

Ian在上面也提到了任意阶段。这个问题已经回答here；我已将此解决方案翻译成下面的代码。 "magic" 发生在调用 ZHEEV 之后。

通过此修复，我没有发现任何不连续之处。

program Eigenvalue
  !> This can be used with the f2008 call
  use, intrinsic :: iso_fortran_env
  implicit none

  !> dp contains the kind value for double precision.
  !> Use below if compiling to f2008
  integer, parameter :: dp = REAL64
  !> Use below if compiling with f95 up and comment out iso_fortran_env
  !>integer, parameter :: dp = SELECTED_REAL_KIND(15, 300)
  !> Set wp to the desired precision.
  integer, parameter :: wp = dp

  integer, parameter :: n = 2
  integer :: i, j, k, m
  real(kind=wp) x, dx, xf, amp, xin 
  real(kind=wp), dimension(n) ::  w
  real(kind=wp), dimension(n, n) ::  Ar, Ai
  complex(kind=wp), dimension(n, n) :: A
  complex(kind=wp), dimension(max(1,2*n-1)) :: WORK
  integer, parameter :: lwork = max(1,2*n-1)
  real(kind=wp), dimension(max(1, 3*n-2)) :: RWORK
  integer :: info
  complex(kind=wp) ::  psi1a, psi1b, psi2a, psi2b 
  real(kind=wp) :: mag
  m = 51              

  xf  = 10.0_wp
  xin = 0.0_wp   

  amp = 2.0_wp

  if (m .eq. 1) then
     dx = 0.0_wp
  else
     dx = (xf - xin)/(m-1)
  end if

  do i = 1, m
     x = dx*(i-m) + xf   

     Ar(1,1) = sin(x)**2
     Ar(1,2) = amp*cos(x)
     Ar(2,1) = amp*cos(x)
     Ar(2,2) = cos(x)**2

     Ai(1,1) = 0.0_wp
     Ai(1,2) = amp*sin(x)  
     Ai(2,1) = -amp*sin(x) 
     Ai(2,2) = 0.0_wp

     do j = 1, n
        do k = 1, n
           A(j, k) = cmplx(Ar(j, k), Ai(j, k), kind=wp)
        end do
     end do
     call ZHEEV('V', 'U', N, A, N, W, WORK, LWORK, RWORK, INFO)
     do j = 1, n
        A(:, j) = A(:, j) / A(1, j)
        mag = sqrt(real(A(1, j)*conjg(A(1, j)))+ real(A(2, j)*conjg(A(2, j))))
        A(:, j) = A(:, j)/mag
     end do

     psi1a = A(1, 1)
     psi1b = A(1, 2)
     psi2a = A(2, 1)
     psi2b = A(2, 2)

  end do
end program Eigenvalue

厄密矩阵对角化中的特征向量不连续性

Eigenvector discontinuities in hermitian matrix diagonalization

fortran

gfortran

fortran77

intel-fortran