廣義對角化與LAPACK

Question

我使用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，我只是嘗試不同的值。

我知道什麼是錯的，因爲我知道什麼是特徵值和這個矩陣的特徵值和我得到錯誤的特徵向量，而我爲了瞭解它的工作方式，然後計算巨大diagonalizations做這樣一個簡單的計算。

有沒有人看到有什麼問題？

謝謝

Answer 1

缺少您發佈的代碼中的某些元素。基本上是WORK，eig，v_mat和t_mat的數組聲明。

無論如何，LWORK參數實際上是WORK矢量的大小。即

DOUBLE PRECISION WORK(100) 
LWORK=100 

Lapack規定最小值爲LWORK=3*N-1。在你的情況下N=2。

對於這個例子，我會建議使用一個大的工作向量（即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) 

此外，你需要經過DSYGV調用查看INF值。如果不是零則發生錯誤。

編輯：固定源代碼

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