将Fortran77代码转换为Matlab代码以找到特征值/向量

Question

我将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 

Answer 1

我看到了一些可能会导致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 

Answer 2

在你的Fortran代码声明一组变量作为REAL。默认情况下，大多数编译器都将它们实现为32位浮点数。 Matlab版本中的相应变量默认为64位浮点数。

没有看到您的输入或输出，很难说这是部分还是全部，这是两个版本输出差异的原因。但是，改变浮点数的精度通常是你报告的问题的原因，尤其是在特征值计算等棘手的数值方法中。

在我写作时，我还观察到你已经翻译了将Fortran相等的浮点值与Matlab进行比较的不良做法。这在Fortran中是不好的做法，在Matlab中是不好的做法。

当你写表达式如

2.*c 

在Fortran中由REAL值2.0乘以标量或数组变量c的每个元件在Matlab小心。在Matlab中，它将标量或矢量c的每个元素乘以整数2，.*是Matlab中的一个运算符（或者，如果您愿意，也可以是运算符的组合），意思是“逐个执行乘法运算”。你已经巧妙地改变了你的程序的语义，并且可能改变了你计算的一些数字的确切值。

这引出了我的最终评论：你称之为Matlab版本只不过是Fortran代码的音译，即从一种语言到另一种语言的逐行复制。你可以在某些时候避开它，但迟早这种天真的方法会让你失望。也许它已经有了。你真的应该把你的Matlab重写成Matlab，就好像你打算编写好的Matlab一样。