將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一樣。