我想學習如何對角化這個簡單的矩陣使用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
難道我做錯了什麼,以獲得不同的特徵值?