2013-12-11 28 views
2

我想學習如何對角化這個簡單的矩陣使用LAPACK:Fortran和MATLAB的返回不同的特徵值相同的矩陣

0.8147 0.9058 0.1270 0.9134 
0.6324 0.0975 0.2785 0.5469 
0.9575 0.9649 0.1576 0.9706 
0.9572 0.4854 0.8003 0.1419 

在Matlab中,我只是使用命令EIG(墊),並獲得輸出:

ans = 

    2.4021 
    -0.0346 
    -0.7158 
    -0.4400 

然而,當我嘗試寫一個簡單的Fortran程序對角化同一矩陣,我得到不同的特徵值:

 implicit none 

    real*8, allocatable::dataMat(:,:),dataMatP(:,:),corMat(:,:), 
$ tempMat(:,:),corMat2(:,:) 
    real*8, allocatable::matList(:),rawData(:) 
    real*8, allocatable ::eig(:),diag(:),offdiag(:),tau(:),work(:) 
    real*8 avg1,avg2,SD1,SD2,geneCorSum,genei,genej,temp 
    integer i,j,k,numElements,info,lwork,numGenes,n, 
$ numExperiments,readsize,numAbsent,count,geneTolerance 

    real*8 mean,std 

    n=4 

    allocate(corMat(4,4)) 

corMat(1,1)=0.8147 
corMat(1,2)=0.9058 
corMat(1,3)=0.1270 
corMat(1,4)=0.9134 
corMat(2,1)=0.6234 
corMat(2,2)=0.0975 
corMat(2,3)=0.2785 
corMat(2,4)=0.5469 
corMat(3,1)=0.9575 
corMat(3,2)=0.9649 
corMat(3,3)=0.1576 
corMat(3,4)=0.9706 
corMat(4,1)=0.9572 
corMat(4,2)=0.4854 
corMat(4,3)=0.8003 
corMat(4,4)=0.1419 



    allocate(diag(n)) 
    allocate(offdiag(n-1)) 
    allocate(tau(n-1)) 
    allocate(work(1)) 

    call dsytrd('U',n,corMat,n,diag,offdiag,tau, 
$ work,-1,info) 
    print*,"Returning from Blocksize calculation" 
    if(info.eq.0) then 
    print*,"Work value successfully calculated:",work(1) 
    endif 
    lwork=work(1) 
    deallocate(work) 
    allocate(work(max(1,lwork))) 

    call dsytrd('U',n,corMat,n,diag,offdiag,tau, 
$ work,lwork,info) 
    print*,"Returning from full SSYTRD" 
    if(info.EQ.0) then 
    print*,"Tridiagonal matrix calculated" 
    endif 



    call dsterf(n,diag,offdiag,info) 
    if(info.EQ.0) then 
    print*,"Matrix Diagonalized" 
    endif 


    do i=1,n 
    print*,"lam=",i,diag(i) 
    enddo 

    deallocate(offdiag) 
    deallocate(work) 
    deallocate(tau) 

    end 

這給了我:

lam= 1, -1.0228376083545221 
lam= 2, -0.48858533844019592 
lam= 3, 0.43828991894506536 
lam= 4, 2.2848330351691031 

難道我做錯了什麼,以獲得不同的特徵值?

回答

3

您使用的LAPACK例程假定有一個對稱矩陣,而原始矩陣是而不是

爲了證明這一點,從原來的矩陣創建一個對稱矩陣,利用右上三角部分和運行MATLAB的eig功能:

for i=1:4 
    for j=i:4; 
    xx(i,j) = x(i,j); 
    xx(j,i)=x(i,j); 
    end 
end 

所產生的基質(x是你有原始的矩陣):

xx = 

0.8147 0.9058 0.1270 0.9134 
0.9058 0.0975 0.2785 0.5469 
0.1270 0.2785 0.1576 0.9706 
0.9134 0.5469 0.9706 0.1419 

和原來x的特徵值和對稱xx矩陣:

>> eig(x) 
    ans = 2.4022 -0.0346 -0.7158 -0.4400 

>> eig(xx) 
    ans = -1.0228 -0.4886  0.4383  2.2848 
1

首先,我希望你不只是複製/粘貼默認情況下Matlab打印命令窗口的小數點後四位。其次,corMat(2,1)=0.6234與第一個矩陣中的對應值不同。第三,dsytrd狀態的文檔:

DSYTRD通過正交相似變換減少了實對稱矩陣A到實對稱三對角形式爲T ...

你的矩陣是絕對不對稱isequal(A,A'))。有很多例程處理非對稱矩陣。例如,您可以嘗試使用dgeev

相關問題