Fortran 子例程 dsygv returns 无穷大

Question

我对解决广义特征值问题的子程序有疑问：

A * x = lambda * B * x,

其中A和B应该是对称矩阵，B是正定的。我正在尝试解决以下测试用例 2x2

    program prog

      implicit none
      integer, parameter :: n=2
      integer :: work,lwork,info
      real(8) :: a(n,n), b(n,n), e(n)

      a(1,1) = -3
      a(1,2) = 4
      a(2,2) = 3
      b(1,1) = 1
      b(1,2) = 0
      b(2,2) = 1

      call dsygv(1, 'v', 'u', n, a, n, b, n, e, WORK, LWORK, INFO)

      write(6,*) info
      write(6,*) e
      write(6,*) a

   end program

并获得以下输出：

    0
    -5.0000000000000000        5.0000000000000000     
    -Infinity  0.44721359549995793  Infinity  0.89442719099991586

特征值是正确的，其中一个特征向量是正确的，但我无法避免无穷大并获得另一个特征向量。

Answer 1

work 和 lwork 应该是数组（双精度和整数）并且必须根据手册适当调整大小。否则 LAPACK 将使用它不应该使用的部分内存，一切都会崩溃。

可能还有其他的论点也是错误的。使用说明查看源代码：http://www.netlib.org/lapack/lapack-3.1.1/html/dsygv.f.html

给 Fortran 初学者的建议：

不要使用 write(6,*)，而是 write(*,*)，不要使用 real(8)，而是 real(dp)，其中 dp 是一个整数常量，右值（使用selected_real_kind()或其他方式获取值，或者如果你坚持将其设置为8）。

Answer 2

下面我建议对您的代码进行更新，其中：

使用可移植的双精度类型值 (wp);
使用 Print;
使用 LAPACK 推荐的查询机制设置（双精度）工作区work。

这里：

Program prog
  Implicit None
  Integer, Parameter :: n = 2
  Integer, Parameter :: wp = kind(0.0D0)
  Integer :: info, lwork
  Real (Kind=wp) :: a(n, n), b(n, n), dummy_work(1), e(n)
  Real (Kind=wp), Allocatable :: work(:)
  External :: dsygv
  Intrinsic :: int, kind, max

  ! Workspace query:
  lwork = -1
  Call dsygv(1, 'v', 'u', n, a, n, b, n, e, dummy_work, lwork, info)
  lwork = int(dummy_work(1))
  Allocate (work(max(1,lwork)))

  a(1, 1) = -3
  a(1, 2) = 4
  a(2, 2) = 3
  b(1, 1) = 1
  b(1, 2) = 0
  b(2, 2) = 1

  Call dsygv(1, 'v', 'u', n, a, n, b, n, e, work, lwork, info)

  Print *, 'info: ', info
  If (info==0) then
    Print *, 'e: ', e
    Print *, 'a: ', a
  end If

End Program

使用我的编译器得到输出

 info:  0
 e:   -5.0000000000000000   5.0000000000000000
 a:   -0.8944271909999159   0.4472135954999579   0.4472135954999579   0.8944271909999159

请注意，一些供应商为 LAPACK（和 BLAS）例程提供显式接口。如果可用，您应该始终使用它们。（这样做可能会在编译时诊断出您输入错误的工作区数组。）

Fortran 子例程 dsygv returns 无穷大

Fortran subroutine dsygv returns infinities

fortran

lapack