4

インプレース LU 行列分解のために、C の Numerical Recipes の提供されたソース コードから文字どおりコピー アンド ペーストしましたが、問題は機能していません。

私は愚かなことをしていると確信していますが、誰かが私を正しい方向に向けることができれば幸いです。私は一日中取り組んできましたが、何が間違っているのかわかりません。

回答後の更新: プロジェクトは終了し、機能しています。皆様のご指導に感謝いたします。

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define MAT1 3
#define TINY 1e-20

int h_NR_LU_decomp(float *a, int *indx){
  //Taken from Numerical Recipies for C
  int i,imax,j,k;
  float big,dum,sum,temp;
  int n=MAT1;

  float vv[MAT1];
  int d=1.0;
  //Loop over rows to get implicit scaling info
  for (i=0;i<n;i++) {
    big=0.0;
    for (j=0;j<n;j++)
      if ((temp=fabs(a[i*MAT1+j])) > big) 
        big=temp;
    if (big == 0.0) return -1; //Singular Matrix
    vv[i]=1.0/big;
  }
  //Outer kij loop
  for (j=0;j<n;j++) {
    for (i=0;i<j;i++) {
      sum=a[i*MAT1+j];
      for (k=0;k<i;k++) 
        sum -= a[i*MAT1+k]*a[k*MAT1+j];
      a[i*MAT1+j]=sum;
    }
    big=0.0;
    //search for largest pivot
    for (i=j;i<n;i++) {
      sum=a[i*MAT1+j];
      for (k=0;k<j;k++) sum -= a[i*MAT1+k]*a[k*MAT1+j];
      a[i*MAT1+j]=sum;
      if ((dum=vv[i]*fabs(sum)) >= big) {
        big=dum;
        imax=i;
      }
    }
    //Do we need to swap any rows?
    if (j != imax) {
      for (k=0;k<n;k++) {
        dum=a[imax*MAT1+k];
        a[imax*MAT1+k]=a[j*MAT1+k];
        a[j*MAT1+k]=dum;
      }
      d = -d;
      vv[imax]=vv[j];
    }
    indx[j]=imax;
    if (a[j*MAT1+j] == 0.0) a[j*MAT1+j]=TINY;
    for (k=j+1;k<n;k++) {
      dum=1.0/(a[j*MAT1+j]);
      for (i=j+1;i<n;i++) a[i*MAT1+j] *= dum;
    }
  }
  return 0;
}

void main(){
    //3x3 Matrix
    float exampleA[]={1,3,-2,3,5,6,2,4,3};
    //pivot array (not used currently)
    int* h_pivot = (int *)malloc(sizeof(int)*MAT1);
    int retval = h_NR_LU_decomp(&exampleA[0],h_pivot);
    for (unsigned int i=0; i<3; i++){
      printf("\n%d:",h_pivot[i]);
      for (unsigned int j=0;j<3; j++){
        printf("%.1lf,",exampleA[i*3+j]);
      }
    }
}

WolframAlphaは、答えは

1,3,-2
2,-2,7
3,2,-2

私は得ています:

2,4,3
0.2,2,-2.8
0.8,1,6.5

これまでのところ、「同じ」アルゴリズムの少なくとも 3 つの異なるバージョンを見つけたので、完全に混乱しています。

PSはい、これを行うには少なくとも12の異なるライブラリがあることは知っていますが、正しい答えよりも自分が間違っていることを理解することに興味があります。

PPS LU 分解では、下の結果の行列は 1 であり、Crouts アルゴリズムを (私が思うに) 実装されているように使用すると、配列インデックスへのアクセスは依然として安全であり、L と U の両方をインプレースで重ね合わせることができます。したがって、これに対する単一の結果行列です。

4

4 に答える 4

10

聖なるものすべてを愛するために、テキストで説明されているアルゴリズムの目的を教えるためのおもちゃの実装として以外は、数値レシピコードを使用しないでください。実際、テキストはそれほど素晴らしいものではありません。そして、あなたが学んでいるように、どちらもコードではありません。

確かに、独自のコードに数値レシピルーチンを含めないでください。特にコードの品質を考えると、ライセンスはめちゃくちゃ制限されています。そこにNRのものがある場合、独自のコードを配布することはできません。

システムにすでにLAPACKライ​​ブラリがインストールされているかどうかを確認します。これは、計算科学および工学における線形代数ルーチンへの標準インターフェイスであり、完璧ではありませんが、コードを移動したすべてのマシンのlapackライブラリを見つけることができ、コンパイル、リンク、および実行することができます。システムにまだインストールされていない場合は、パッケージマネージャー(rpm、apt-get、fink、portなど)がlapackについて知っている可能性があり、インストールできます。そうでない場合は、システムにFortranコンパイラーがあれば、ここからダウンロードしてコンパイルできます。標準のCバインディングは、同じページのすぐ下にあります。

線形代数ルーチンへの標準APIがあると非常に便利な理由は、それらが非常に一般的であるが、それらのパフォーマンスはシステムに非常に依存しているためです。したがって、たとえば、Goto BLAS は、線形代数に必要な低レベル操作のx86システムの非常に高速な実装です。LAPACKが機能するようになったら、そのライブラリをインストールして、すべてを可能な限り高速化できます。

LAPACKをインストールすると、一般的な行列のLU分解を実行するルーチンは、floatの場合はSGETRF、doubleの場合はDGETRFになります。行列の構造について何か知っている場合は、他にも高速なルーチンがあります。たとえば、対称正定値(SBPTRF)、または三重対角(STDTRF)です。これは大きなライブラリですが、それを回避する方法を学ぶと、数値ツールボックスに非常に強力なギアが含まれるようになります。

于 2011-04-17T03:03:36.327 に答える
10

あなたのインデックスには本質的に何か問題があると思います。jそれらは時々異常な開始値と終了値を持ち、代わりに外側のループがi疑わしくなります。

誰かにコードを調べてもらう前に、いくつかの提案を以下に示します。

  • インデックスを再確認してください
  • を使用してこれらの難読化の試みを取り除きますsum
  • a(i,j)の代わりにマクロを使用するa[i*MAT1+j]
  • コメントの代わりにサブ関数を書く
  • 不要な部分を削除し、エラーのあるコードを分離します

これらの提案に従うバージョンは次のとおりです。

#define MAT1 3
#define a(i,j) a[(i)*MAT1+(j)]

int h_NR_LU_decomp(float *a, int *indx)
{
        int i, j, k;
        int n = MAT1;

        for (i = 0; i < n; i++) {
                // compute R
                for (j = i; j < n; j++)
                        for (k = 0; k < i-2; k++)
                                a(i,j) -= a(i,k) * a(k,j);

                // compute L
                for (j = i+1; j < n; j++)
                        for (k = 0; k < i-2; k++)
                                a(j,i) -= a(j,k) * a(k,i);
        }

        return 0;
}

その主な利点は次のとおりです。

  • それは読める
  • できます

ただし、ピボットがありません。必要に応じてサブ機能を追加します。


私のアドバイス: 他人のコードを理解せずにコピーしないでください。

ほとんどのプログラマーは下手なプログラマーです。

于 2011-04-17T01:43:32.797 に答える
3

私にとって最も疑わしいのは、「最大ピボットの検索」とマークされた部分です。これは検索するだけでなく、行列 A も変更します。これが正しいとは信じがたいです。

異なるバージョンの LU アルゴリズムではピボットが異なるため、それを理解していることを確認してください。異なるアルゴリズムの結果を比較することはできません。より良いチェックは、L x U が元の行列に等しいかどうか、またはアルゴリズムがピボットを行う場合はその順列に等しいかどうかを確認することです。そうは言っても、行列式が間違っているため、結果は間違っています(ピボットを行っても、符号を除いて行列式は変わりません)。

それとは別に、@ Philipには良いアドバイスがあります。コードを理解したい場合は、ピボットせずに LU 分解を理解することから始めてください。

于 2011-04-17T10:31:51.353 に答える
3

アルバート・アインシュタインをひどく言い換えると:

...時計を持っている人は常に正確な時間を知っていますが、2つ持っている人は決して正確ではありません....

あなたのコードは間違いなく正しい結果を生成していませんが、たとえそうであったとしても、ピボットを使用した結果は、ピボットを使用しない結果に直接対応しません。ピボット ソリューションのコンテキストでは、Alpha が実際に提供したものは、おそらくこれと同等です。

    1  0  0      1  0  0       1  3 -2
P=  0  1  0   L= 2  1  0  U =  0 -2  7
    0  0  1      3  2  1       0  0 -2

これは条件 A = PLU を満たします (ここで . は行列の積を表します)。(概念的に) 同じ分解操作を別の方法で計算すると (この場合は numpy を介して LAPACK ルーチン dgetrf を使用):

In [27]: A
Out[27]: 
array([[ 1,  3, -2],
       [ 3,  5,  6],
       [ 2,  4,  3]])

In [28]: import scipy.linalg as la

In [29]: LU,ipivot = la.lu_factor(A)

In [30]: print LU
[[ 3.          5.          6.        ]
 [ 0.33333333  1.33333333 -4.        ]
 [ 0.66666667  0.5         1.        ]]

In [31]: print ipivot
[1 1 2]

ipivot でちょっとした黒魔術を試してみると、

    0  1  0       1        0    0       3  5       6
P = 0  0  1   L = 0.33333  1    0  U =  0  1.3333 -4
    1  0  0       0.66667  0.5  1       0  0       1

これは A = PLU も満たします。これらの因数分解は両方とも正しいですが、それらは異なっており、正しく機能するバージョンの NR コードには対応していません。

したがって、「正しい」答えを持っているかどうかを判断する前に、コピーしたコードが実装する実際のアルゴリズムを理解するのに少し時間を費やす必要があります。

于 2011-04-18T09:56:09.640 に答える