仅使用一个矩阵在 C 中实现具有部分旋转的 LU 分解

Question

为了计算 PA = LU 分解，我设计了以下 C 函数，仅使用一个矩阵来存储和计算数据：

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

  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) {
      c[k][j] /= c[j][j];
      c[k][k] -= c[k][j]*c[j][k];
    }

  }
  return 0.;
}

其中 n 是矩阵的阶，c 是指向矩阵的指针，p 是指向存储在部分旋转系统时完成的排列的向量的指针。变量 tol 目前不相关。该程序将分解的下三角部分和上三角部分都存储在 c 中，其中 U 对应于 c 的上三角部分，L 对应于 c 的严格下三角部分，在对角线上添加 1。对于我能够测试的内容，与部分旋转对应的程序部分工作正常，但是，用于计算矩阵条目的算法没有给出预期的结果，我不明白为什么。例如，如果我尝试计算矩阵的 LU 分解

1. 2. 3.
4. 5. 6.
7. 8. 9.

我明白了

    1.    0.     0.         7. 8. 9.
l : 0.143 1.     0.     u : 0. 2. 1.714*
    0.571 0.214* 1.         0. 0. 5.663*

其乘积不对应于矩阵 c 的任何排列。事实上，错误的条目似乎是标有星号的条目。

如果有任何解决此问题的建议，我将不胜感激。

Answer 1

我发现您的代码存在问题，您在计算实际分解时对行进行规范化的方式存在一些概念性错误：

for (k = j+1; k < n; ++k) {
  c[k][j] /= c[j][j];
  c[k][k] -= c[k][j]*c[j][k];
}

成为：

for (k = j+1; k < n; ++k) {
  temp=c[k][j]/=c[j][j];
  for(int q=j+1;q<n;q++){
        c[k][q] -= temp*c[j][q];
      }
}

返回结果：

7.000000 8.000000 9.000000
0.142857 0.857143 1.714286
0.571429 0.500000 -0.000000

如果您有任何问题，我很乐意为您提供帮助。

在这里完整实现：​​

#include<stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

double plupmc(int n, double **c, int *p, double tol) {
  int i, j, k, pivot_ind = 0, temp_ind;
  int ii, jj;
  double *vv=calloc(n,sizeof(double));
  double pivot, *temp_row;
  double temp;

  for (j = 0; j < n; ++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) {
      temp=c[k][j]/=c[j][j];
      for(int q=j+1;q<n;q++){
            c[k][q] -= temp*c[j][q];
          }
    }
    for(int q=0;q<n;q++){
      for(int l=0;l<n;l++){
        printf("%lf ",c[q][l]);
      }
      printf("\n");
    }

  }
  return 0.;
}

int main() {
  double **x;
  x=calloc(3,sizeof(double));
  for(int i=0;i<3;i++){
    x[i]=calloc(3,sizeof(double));
  }
  memcpy(x[0],(double[]){1,2,3},3*sizeof(double));
  memcpy(x[1],(double[]){4,5,6},3*sizeof(double));
  memcpy(x[2],(double[]){7,8,9},3*sizeof(double));

  int *p=calloc(3,sizeof(int));
  memcpy(p,(int[]){0,1,2},3*sizeof(int));
  plupmc(3,x,p,1);
  for(int i=0;i<3;i++){
    free(x[i]);
  }
  free(p);
  free(x);


}