我将 fortran 77 中的书面代码转换为 Matlab 代码。此函数使用 QL 算法计算矩阵的特征值和特征向量。由于某些原因,我不能在 matlab 中使用 eig 函数的结果。该方法得到的特征值与 eig 函数得到的特征值不完全相同,有的相同,有的不同。我不知道问题出在哪里。感谢您的任何帮助和建议。如果需要运行和查看结果,我可以提供输入数组。
这是fortran代码:
SUBROUTINE tqli(d,e,n,np,z)
INTEGER n,np
REAL d(np),e(np),z(np,np)
CU USES pythag
INTEGER i,iter,k,l,m
REAL b,c,dd,f,g,p,r,s,pythag
do 11 i=2,n
e(i-1)=e(i)
11 continue
e(n)=0.
do 15 l=1,n
iter=0
1 do 12 m=l,n-1
dd=abs(d(m))+abs(d(m+1))
if (abs(e(m))+dd.eq.dd) goto 2
12 continue
m=n
2 if(m.ne.l)then
if(iter.eq.30)pause 'too many iterations in tqli'
iter=iter+1
g=(d(l+1)-d(l))/(2.*e(l))
r=pythag(g,1.)
g=d(m)-d(l)+e(l)/(g+sign(r,g))
s=1.
c=1.
p=0.
do 14 i=m-1,l,-1
f=s*e(i)
b=c*e(i)
r=pythag(f,g)
e(i+1)=r
if(r.eq.0.)then
d(i+1)=d(i+1)-p
e(m)=0.
goto 1
endif
s=f/r
c=g/r
g=d(i+1)-p
r=(d(i)-g)*s+2.*c*b
p=s*r
d(i+1)=g+p
g=c*r-b
C Omit lines from here ...
do 13 k=1,n
f=z(k,i+1)
z(k,i+1)=s*z(k,i)+c*f
z(k,i)=c*z(k,i)-s*f
13 continue
C ... to here when finding only eigenvalues.
14 continue
d(l)=d(l)-p
e(l)=g
e(m)=0.
goto 1
endif
15 continue
return
END
以下是matlab代码:
function [d,z]=tqli(d,e,n,np,z)
for i=2:n
e(i-1)=e(i);
end
e(n)=0.;
for l=1:n
iter=0;
for m=l:(n-1)
dd=abs(d(m))+abs(d(m+1));
if ((abs(e(m))+dd)==dd)
break
end
end
m=n;
if (m~=l)
if (iter==30)
disp('too many iteration in tqli')
end
iter=iter+1;
g=(d(l+1)-d(l))/(2.*e(l));
r=pythag(g,1.);
g=d(m)-d(l)+e(l)/(g+r*sign(g));
s=1.;
c=1.;
p=0.;
for i=(m-1):-1:l
f=s*e(i);
b=c*e(i);
r=pythag(f,g);
e(i+1)=r;
if(r==0.)
d(i+1)=d(i+1)-p;
e(m)=0.;
break
end
s=f/r;
c=g/r;
g=d(i+1)-p;
r=(d(i)-g)*s+2.*c*b;
p=s*r;
d(i+1)=g+p;
g=c*r-b;
for k=1:n
f=z(k,i+1);
z(k,i+1)=s*z(k,i)+c*f;
z(k,i)=c*z(k,i)-s*f;
end
end
d(l)=d(l)-p;
e(l)=g;
e(m)=0.;
end
end
end