0

Fortran 77 で書かれたコードを Matlab コードに変換しました。この関数は、QL アルゴリズムを使用して行列の固有値と固有ベクトルを計算します。何らかの理由で、eig 関数の結果を matlab で使用できません。このメソッドから取得された固有値は、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
4

2 に答える 2

1

matlabの翻訳で問題を引き起こす可能性のあるものをいくつか見ました。1つは、Fortranサインの変換でした。rだけでなくabs（r）を使用する必要があります。

私が見たもう1つのより深刻な問題は、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
于 2012-07-30T11:44:12.253 に答える
1

Fortran コードでは、一連の変数を として宣言しますREAL。ほとんどのコンパイラは、デフォルトで、これらを 32 ビット浮動小数点数として実装します。Matlab バージョンの対応する変数は、デフォルトで 64 ビット浮動小数点数です。

入力または出力のいずれかを確認しないと、これが 2 つのバージョンの出力の違いの原因の一部または全体であるかどうかを判断するのは困難です。ただし、浮動小数点数の精度を変更すると、特に固有値計算などのトリッキーな数値メソッドで、報告されたような問題の原因になることがよくあります。

私が書いている間、浮動小数点値の等価性を比較するという悪い習慣をFortranからMatlabに翻訳したことにも気付きました。これは Fortran の悪い習慣であり、Matlab の悪い習慣です。

次のような式を記述するときは、Matlab で注意してください。

2.*c

Fortran では、スカラー変数または配列変数の各要素にvalueを乗算cします。Matlab では、スカラーまたはベクトルの各要素を整数で乗算します。これは、Matlab の演算子 (または必要に応じて演算子の組み合わせ) であり、「要素ごとに乗算を実行する」ことを意味します。プログラムのセマンティクスを微妙に変更し、計算する数値の正確な値を変更した可能性があります。REAL2.0c2.*

これは私の最後のコメントにつながります.Matlabバージョンと呼ばれるものは、Fortranコードの音訳、つまりある言語から別の言語への行ごとのコピーにすぎません。場合によってはうまくいくこともありますが、遅かれ早かれ、この素朴なアプローチはあなたをつまずかせるでしょう。おそらくそれはすでに持っています。良い Matlab を書くつもりであるかのように、Matlab を Matlab に書き直す必要があります。

于 2012-07-29T16:24:49.397 に答える