2015-04-17 83 views
4

我使用LAPACK庫中的DSYEV和DSYEVD來查找特徵值和特徵向量(編譯語法:gfortran -llapack)。但是,我發現特定矩陣的錯誤特徵值(-0.44,0.35,0.88)。出了什麼問題?LAPACK給我不正確的特徵值

人們可以很容易地看到矩陣有零行列式,所以至少有一個特徵值必須爲零。

這裏是我的代碼(希望這不是太大):

Program Real_Eigenvec 
    implicit none 

    integer, parameter:: n=3 
    integer:: i,j, flag 
    real*8:: A(n,n),X(n,n) 
    real*8:: lambda(n) 
    real*8, parameter:: p=0.5d0/dsqrt(2.d0), q=1.d0-1.d0/dsqrt(2.d0) 


    Print*,'Enter flag: 0 for DSYEV, 1 for DSYEVD' 
    Read*, flag 


    A= transpose(reshape((/ 0.d0, 1.d0, 0.d0, p, q, p, 0.5d0, 0.0d0, 0.5d0 /), shape(A))) 


    print*,'Dimension of the matrix, n=',int(sqrt(float(size(A)))) 

    Print*,'A matrix in full form:' 
    Do i=1,n 
     print 100, (A(i,j),j=1,n) 
    End Do 

    call Eigen(A,lambda,X,n,flag) 

    ! Print the eigenvalues and eigenvectors. 

    PRINT 200 
    DO i = 1, n 
     PRINT 201, lambda(i), (X(i,j), j=1,n) 
    END DO 

100 FORMAT (1X,10(:2X,F10.2)) 
200 FORMAT (/1X, 'Eigenvalue', 16X, 'Eigenvector^T') 
201 FORMAT (1X, F10.2,4X,6(:f10.2)) 

    End Program Real_Eigenvec 



!!!  SUBROUTINES ----------------------------------------- 


    Subroutine Eigen(A,lambda,X,n,flag) 
    implicit none 

    integer:: i,n,flag 
    real*8:: A(n,n),Ap(n,n),X(n,n) 
    real*8:: lambda(n) 
    real*8, allocatable :: work(:) 
     integer, allocatable :: iwork(:) 
    integer:: lwork,liwork,info 

    print*,'n in Eiegen routine=',n 

    lwork=3*n-1  ! DSYEV for flag=0 
    if (flag==1) then ! DSYEVD for flag=1 
     lwork=1+6*n+2*n**2 
    end if 

    liwork=3+5*n 


    allocate(work(lwork)) 
    allocate(iwork(liwork)) 

    Ap=A 

    if (flag==0) then 
     CALL DSYEV ('v', 'l', n, Ap, n, lambda, work, lwork, info) 
    else 
     CALL DSYEVD ('V', 'U', n, Ap, n, lambda, work, & 
       & lwork, iwork, liwork, info) 
     ! For doumentation visit: http://www.netlib.org/lapack/explore-html/d1/da2/dsyevd_8f.html 
    end if 

    X=Ap 

    print*,'info=',info 

    deallocate(work) 
    deallocate(iwork) 

    End Subroutine Eigen 
+0

沒有多少人關注[tag:fortran90],如果你想有人幫助你,最好只在必要時使用[tag:fortran]和版本標籤(這裏不是這種情況)。 –

+0

謝謝@VladimirF。我將爲未來的職位保留這一點。 – hbaromega

回答

7

如LAPACK documentationDSYEV可用於對稱矩陣表示。

DSYEV計算所有特徵值和,任選地,一個真正的對稱矩陣A.

的特徵向量在該例子中矩陣A不是對稱

Dimension of the matrix, n=   3 
A matrix in full form: 
    0.00  1.00  0.00 
    0.35  0.29  0.35 
    0.50  0.00  0.50 

在這種情況下,你應使用DGEEV用於非對稱matrices

DGEEV計算N乘N實非對稱矩陣A,特徵值和可選的左和/或右特徵向量。

一般使用的特徵值是複數,所以你必須提供WRWL。此外,您需要定義是否要離開VL或右側的VR特徵向量。

A * v(j) = lambda(j) * v(j) 
u(j)**H * A = lambda(j) * u(j)**H 

函數的定義是:

DGEEV(JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR, LDVR, WORK, LWORK, INFO) 

我建議使用它像

LWORK = 4*N 
CALL DGEEV('N', 'V', n, A, n, wr, wl, Ap, n, Ap, n, work, lwork, info) 

爲了獲得左,右特徵向量使用

real*8:: A(n,n),VL(n,n),VR(n,n) 
real*8:: wr(n),wl(n) 
lwork = 4*N 
allocate(work(lwork)) 
CALL DGEEV('V', 'V', n, A, n, wr, wl, VL, n, VR, n, work, lwork, info) 

對於你的矩陣imagi所有特徵值的nary部分是零。所以特徵值是(1.00, -0.21, 0.00)

+0

非常感謝。我完全不知道DSYEV只適用於對稱矩陣。然而,使用'CALL DGEEV'('N','V',n,A,n,wr,wl,Ap,n,Ap,n,work,lwork,info)'我仍然收到錯誤消息: **在進入DGEEV時,參數號13有非法值 info = -13'。我已經宣佈了'real * 8 :: wr(n),wi(n)'。 – hbaromega

+1

該消息告訴您檢查參數號13.您是否這樣做?哪一個是它的價值? –

+0

既然你需要特徵向量@hbaromega,你必須設置LWORK> = 4 * N。我會在對 – ztik

相關問題