C语言单矩阵实现带部分主元PA=LU分解结果异常求修复建议

问题根源 #

你代码的核心错误出在消元计算逻辑：选主元完成后，内层循环仅对c[k][k]单个元素做了更新，完全不符合LU分解的高斯消元要求。
正确的逻辑是：计算得到L矩阵第k行第j列的乘数c[k][j]后，需要用该乘数对第k行从j+1到n-1的所有上三角元素做消元更新，而非仅更新对角位置的c[k][k]。

另外你代码还有两个隐藏问题：

• 未引入math.h头文件，调用fabs会触发未定义行为
• 没有初始化置换向量p，如果外部调用时未提前给p赋值为0~n-1的序列，置换结果会完全错误

修复后的代码 #

#include <math.h>

double plupmc(int n, double **c, int *p, double tol) {
  int i, j, k, pivot_ind = 0, temp_ind;
  double pivot, *temp_row;

  // 初始化置换向量，若要求外部初始化可删除本段
  for (i = 0; i < n; i++) {
    p[i] = i;
  }

  for (j = 0; j < n-1; ++j) {
    // 选主元逻辑保留，原有实现正确
    pivot = 0.;
    for (i = j; i < n; ++i) {
      if (fabs(c[i][j]) > fabs(pivot)) {
        pivot = c[i][j];
        pivot_ind = i;
      }
    }

    // 行交换逻辑保留，原有实现正确
    temp_row = c[j];
    c[j] = c[pivot_ind];
    c[pivot_ind] = temp_row;

    temp_ind  = p[j];
    p[j] = p[pivot_ind];
    p[pivot_ind] = temp_ind;

    // 消元逻辑完全重写
    for (k = j+1; k < n; ++k) {
      // 计算L矩阵的严格下三角元素，原有该行逻辑正确
      c[k][j] /= c[j][j];
      // 新增内层循环，更新第k行所有U矩阵元素
      for (int t = j+1; t < n; t++) {
        c[k][t] -= c[k][j] * c[j][t];
      }
    }
  }
  return 0.;
}

验证结果 #

用你给出的测试矩阵[[1,2,3],[4,5,6],[7,8,9]]测试修复后的代码：

• 存储在c矩阵中的计算结果为：

7.000   8.000   9.000
0.143   0.857   1.714
0.571   0.500   0.000

拆分得到的L、U矩阵：

• L = [[1, 0, 0], [0.143, 1, 0], [0.571, 0.5, 1]]
• U = [[7, 8, 9], [0, 0.857, 1.714], [0, 0, 0]]
• 置换向量p = [2,0,1]

三者计算P^-1*L*U的结果与原矩阵完全一致，符合PA=LU的分解要求。

内容的提问来源于stack exchange，提问作者Marc