我将一段用 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
最佳答案
我看到一些可能会给您的 matlab 翻译带来麻烦的事情。一个是您对 fortran 符号的转换。您需要使用 abs(r) 而不仅仅是 r。
我看到的另一个更严重的问题是您尝试重构 goto 导致的流程结构。当我使用 f2matlab(以及我编写的一个未发布的工具 remgoto)转换它时,它产生了以下流程结构。我希望这对您有所帮助。请注意,这一切都未经测试!
function [d,e,n,np,z]=tqli(d,e,n,np,z);
remg([1:2])=true;
for i = 2 : n;
e(i-1) = e(i);
end
e(n) = 0.;
for l = 1 : n;
while (1);
if(remg(2))
if(remg(1))
iter = 0;
end;
remg(1)=true;
for m = l : n - 1;
dd = abs(d(m)) + abs(d(m+1));
if( abs(e(m))+dd==dd )
remg(2)=false;
break;
end;
end;
if(~(remg(2)))
continue;
end;
m = fix(n);
end;
remg(2)=true;
if( m~=l )
if( iter==30 )
disp(['too many iterations in tqli',' -- Hit Return to continue']);
pause ;
end;
iter = fix(iter + 1);
g =(d(l+1)-d(l))./(2..*e(l));
r = pythag(g,1.);
g = d(m) - d(l) + e(l)./(g+(abs(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.;
remg(1)=false;
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;
% Omit lines from here ...
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; k = fix(n+1);
% ... to here when finding only eigenvalues.
end;
if(~(remg(1)))
continue;
end;
d(l) = d(l) - p;
e(l) = g;
e(m) = 0.;
remg(1)=false;
continue;
end;
break;
end;
end;
end %subroutine tqli
关于algorithm - 将 Fortran77 代码转换为 Matlab 代码以查找特征值/向量,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/11710366/