下三角矩阵和上三角矩阵给我的错误答案

Question

我正在LU分解中C.My代码非常简单 算法可以使用两个循环并行化一个用于更新下三角矩阵和一个用于更新上三角矩阵，但似乎我想了解的东西:(

 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问题，它提供了一些更有用的提示。