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

Question

我正在使用 C 语言进行 LU 分解。我的代码非常简单 算法可以使用两个循环并行化，一个用于更新下三角矩阵，一个用于 更新上三角矩阵，但似乎我错过了一些了解:(

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