下三角行列と上三角行列は私に間違った答えを与える

Question

私はLUの分解をCで行っています。私のコードは非常に簡単です アルゴリズムは2つのループを使用して1つの下位三角行列を更新します。行列が、私が何かを理解し欠場思わ:(

 for (i=0 ; i<N ; i++){ 
// A[i][i]=1; 
    for (j=i+1 ;j<N ;j++){ 
     A[j][i] = A[j][i]/A[i][i]; //*Update L*// 
    } 
    for (j=i+1;j<N;j++){ 
     for(k=i+1 ;k<N;k++){ 

      A[j][k] = A[j][k] - A[i][k] * A[j][i];//*Update U*// 
     } 
    } 
} 

    printf("\n Matrix after U transformation: \n"); 
    print_matrix(); 

for (i=0 ; i<N ; i++){ 
    A[i][i]=1; 
    for (j=i+1 ;j<N ;j++){ 
     A[j][i] = A[j][i]/A[i][i]; //*Update L*// 
    } 
    for (j=i+1;j<N;j++){ 
     for(k=i+1 ;k<N;k++){ 

      A[j][k] = A[j][k] - A[i][k] * A[j][i];//*Update U*// 
     } 
     } 
    } 

    printf("\n Matrix after L transformation: \n"); 
    print_matrix(); 

This is what I should to get ?! what I'm doing wrong 

L = 

1.0000   0   0   0   0 
0.2000 1.0000   0   0   0 
0.2000 0.1667 1.0000   0   0 
0.2000 0.1667 0.1429 1.0000   0 
0.2000 0.1667 0.1429 0.1250 1.0000 


U = 

50.0000 10.0000 10.0000 10.0000 10.0000 
    0 48.0000 8.0000 8.0000 8.0000 
    0   0 46.6667 6.6667 6.6667 
    0   0   0 45.7143 5.7143 
    0   0   0   0 45.0000 

が、私が得たもの,,,, Lは次のようにあるべきではない？！

Source Matrix : 
50  10  10  10  10 
10  50  10  10  10 
10  10  50  10  10 
10  10  10  50  10 
10  10  10  10  50 

Matrix after U transformation: 
50  10  10  10  10 
    0  48  8  8  8 
    0  0  47  7  7 
    0  0  0  46  6 
    0  0  0  0  45 

Matrix after L transformation: 
    1  10  10  10  10 
    0  1  6  6  6 
    0  -2  1  16  16 
    0  -2  9  1 -129 
    0  -2  9 -134  1 

おかげ

Answer 1

あなたのUマトリックスは正しいですが、これは整数で浮動小数点ではありません。あなたのL行列の対角も正しいです（値を設定しています）が、残りはそうではありません。 （いくつかの括弧を追加し、それを少し変更）である「LU Decomposition from Numerical Recipes not working; what am I doing wrong?」の答えに対するコード確認後：

for (i = 0; i < N; i++) { 
    // compute U 
    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]; 
     } 
    } 
} 

を私は問題であるべき、あなたのコード内のループを抜けていることに気づきました。また、いくつかのより有用なヒントを提供する上記のSOの質問を見てください。