2

我在使用 Lapack 的子程序 DSYGV 时遇到了一些问题：

DSYGV( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,LWORK, INFO )

这是我要进行的对角化：

v_mat*x = eig*t_mat*x

这是我的代码的关键部分：

program pruebadiago

real, dimension(:,:), allocatable :: v_mat, t_mat

real, dimension(:), allocatable   :: eig,WORK

real , parameter                  :: k=3.0,m=4.0

integer, parameter                :: n=2

integer                           :: i

! EXPECTED EIGENVALUES AND EIGENVECTORS

!eig = 0.286475  ------>  u = (0.262866 , 0.425325)

!eig = 1.96353   ------>  u = (0.425325, -0.262866)

allocate(v_mat(n,n),t_mat(n,n),eig(n))

!--------------------------

v_mat(1,1:2) = (/2.0*k,-k/)

v_mat(2,1:2) = (/-k,k/)

!--------------------------

!--------------------------

t_mat(1,1:2) = (/m,0.0/)

t_mat(2,1:2) = (/0.0,m/)

!---------------------------

!diagonalizacion

call DSYGV( 1, 'v', 'u', n , v_mat, 2 , t_mat, 2, eig, WORK,-1, INF )

LWORK=WORK(1)

allocate(WORK(LWORK))

call  DSYGV( 1, 'v', 'u', n , v_mat, 2 , t_mat, 2, eig, WORK,LWORK, INF )

open(unit=100,file="pruebadiago.dat",status="replace",action="write")

do i = 1,n

write(unit=100,fmt=*) "E",i,"=",eig(i),(v_mat(i,j),j=1,n)

!autofuntzioak zutabeka doaz"(100f12.6)"

enddo

close(unit=100)

deallocate(v_mat,t_mat,eig,WORK)

end program pruebadiago

我想我理解了本文档中给出的所有内容：

http://www.netlib.org/lapack/lapack-3.1.1/html/dsygv.f.html

但是我不理解的论点 LWORK 所以我只是尝试不同的值。

我知道出了点问题，因为我知道这个矩阵的特征值和特征向量是什么，我得到了错误的特征值和特征向量，我正在做一个简单的计算，以了解它的工作方式，然后计算巨大的对角化。

有没有人看到有什么问题？

谢谢

4

1 回答 1

3

您发布的代码中缺少一些元素。基本上是 WORK、eig、v_mat 和 t_mat 的数组声明。

无论如何，LWORK 参数实际上是 WORK 向量的大小。IE

DOUBLE PRECISION WORK(100)
LWORK=100

Lapack 指定为 的最小值LWORK=3*N-1。在你的情况下N=2。

对于这个例子，我建议使用一个大的 WORK 向量（即 100），这样你就不会遇到任何问题。

对于大型矩阵，您应该使用双重调用DSYGV.

  • 第一次打电话LWORK=-1并从中获得建议的尺寸WORK(1)
  • 分配NEW WORK具有建议大小的向量。
  • 最后用 求解特征问题DSYGV。

示例代码是：

CALL DSYGV( 1, 'V', 'U', 2 , v_mat, 2 , t_mat, 2, eig, W, -1, INF )
LWORK=W(1)
ALLOCATE ( WORK(LWORK) )
CALL DSYGV( 1, 'V', 'U', 2 , v_mat, 2 , t_mat, 2, eig, WORK, LWORK, INF )

此外，您需要在调用INF后检查值。DSYGV如果它不为零，则发生错误。

编辑：固定源代码

program pruebadiago

double precision, dimension(:,:), allocatable :: v_mat, t_mat
double precision, dimension(:), allocatable :: eig,WORK
double precision :: W
double precision , parameter :: k=3.0,m=4.0
integer, parameter :: n=2
integer :: i

! EXPECTED EIGENVALUES AND EIGENVECTORS
!eig = 0.286475 ------> u = (0.262866 , 0.425325)
!eig = 1.96353 ------> u = (0.425325, -0.262866)

allocate(v_mat(n,n),t_mat(n,n),eig(n))

!--------------------------
v_mat(1,1:2) = (/2.0*k,-k/)
v_mat(2,1:2) = (/-k,k/)
!--------------------------
t_mat(1,1:2) = (/m,0.0d0/)
t_mat(2,1:2) = (/0.0d0,m/)
!---------------------------

call DSYGV( 1, 'v', 'u', n , v_mat, 2 , t_mat, 2, eig, W,-1, INF )

LWORK=W
allocate(WORK(LWORK))

call DSYGV( 1, 'v', 'u', n , v_mat, 2 , t_mat, 2, eig, WORK,LWORK, INF )
open(unit=100,file="pruebadiago.dat",status="replace",action="write")
do i = 1,n
  write(unit=100,fmt=*) "E",i,"=",eig(i),(v_mat(i,j),j=1,n)
enddo

close(unit=100)
deallocate(v_mat,t_mat,eig,WORK)

end program pruebadiago
于 2013-08-23T07:55:49.607 回答