Fortran77 코드를 Matlab으로 변환 고유 값/벡터를 찾는 코드

Question

Fortran 77에서 작성된 코드를 Matlab 코드로 변환했습니다. 이 함수는 QL 알고리즘을 사용하여 행렬의 고유 값과 고유 벡터를 계산합니다. 몇 가지 이유로 나는 matlab에서 eig 함수의 결과를 사용할 수 없다. 이 방법으로 얻은 고유 값은 eig 함수로 얻은 것과 동일하지 않으며 일부는 동일하지만 일부는 다릅니다. 나는 문제가 어디 있는지 모른다. 어떤 도움이나 제안을 해주셔서 감사합니다. 결과를보고 실행하는 데 필요한 경우 입력 배열을 제공 할 수 있습니다.

여기 포트란 코드 :

 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 번역에 문제를 일으킬 수있는 몇 가지를 보았다. 하나는 포트란 간판의 전환이었습니다. r 대신 abs (r)을 사용해야합니다.

내가 본 다른 심각한 문제는 goto 's를 리팩토링하려는 노력에서 비롯된 흐름 구조입니다. 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

당신이 REAL으로 변수의 무리를 선언합니다. 대부분의 컴파일러는 기본적으로 이것을 32 비트 부동 소수점 숫자로 구현합니다. Matlab 버전의 해당 변수는 기본적으로 64 비트 부동 소수점 숫자입니다.

입력 또는 출력을 보지 않고도 이것이 부분적으로 또는 전체적으로 두 버전의 출력 차이의 원인인지 여부를 구분하기가 어렵습니다. 그러나 부동 소수점 숫자의 정밀도를 변경하는 것은 종종 고유 한 계산과 같은 까다로운 수치 방법에서보고하는 문제의 원인입니다.

필자도 필자는 부동 소수점 값을 Fortran에서 Matlab로 비교하는 나쁜 습관을 번역했다고보고합니다. Fortran에서는 나쁜 습관을 가지고 있으며 Matlab에서는 나쁜 습관을 가지고 있습니다. 이 같은 REAL 값 2.0하여 스칼라 또는 배열 변수 c의 각 요소를 곱 포트란에서

2.*c 

와 같은 표현을 쓸 때

는 매트랩주의하십시오. Matlab에서는 스칼라 또는 벡터 인 c의 각 요소에 정수 2을 곱하고 Matlab에서 .*은 연산자 (또는 원하는 경우 연산자의 조합)로 '요소별로 곱셈을 수행'을 의미합니다. 프로그램의 의미와 아마도 계산 한 숫자의 정확한 값을 미묘하게 변경했습니다.

Matlab 버전이라고하는 것은 한 언어에서 다른 언어로 행 단위로 복사하는 Fortran 코드의 음역입니다. 당신은 시간의 일부를 가지고 도망 갈 수 있습니다. 그러나 조만간이 순진한 접근법이 당신을 떠올리게 할 것입니다. 아마도 그것은 이미 있습니다. 좋은 Matlab을 작성하려는 경우 실제로 Matlab을 Matlab으로 다시 작성해야합니다.