我试图通过对这个简单的矩阵进行对角化来学习如何使用LaPACK:
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(mat),得到输出:
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我是不是做错了什么来得到不同的特征值?
发布于 2013-12-11 18:55:29
您使用的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发布于 2013-12-11 18:56:07
首先,我希望您不只是复制/粘贴Matlab默认从命令窗口打印出来的四位小数点。其次,corMat(2,1)=0.6234与第一个矩阵中的对应值不同。第三,dsytrd的文档声明:
DSYTRD通过正交相似变换将实对称矩阵A降为实对称三对角形式T。
你的矩阵绝对是不对称的(isequal(A,A'))。有处理非对称矩阵的各种例程。例如,您可以尝试使用dgeev。
发布于 2013-12-11 19:36:58
SSYTRD/DSYTRD只适用于对称矩阵。
https://stackoverflow.com/questions/20527031
复制相似问题