4

我已经从提供的源代码中复制并粘贴了 C 的数字配方用于就地 LU 矩阵分解，问题是它不起作用。

我确定我在做一些愚蠢的事情，但希望有人能够指出我正确的方向；我整天都在努力，看不出我做错了什么。

回复后更新：该项目已完成并正在运行。感谢大家的指导。

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define MAT1 3
#define TINY 1e-20

int h_NR_LU_decomp(float *a, int *indx){
  //Taken from Numerical Recipies for C
  int i,imax,j,k;
  float big,dum,sum,temp;
  int n=MAT1;

  float vv[MAT1];
  int d=1.0;
  //Loop over rows to get implicit scaling info
  for (i=0;i<n;i++) {
    big=0.0;
    for (j=0;j<n;j++)
      if ((temp=fabs(a[i*MAT1+j])) > big) 
        big=temp;
    if (big == 0.0) return -1; //Singular Matrix
    vv[i]=1.0/big;
  }
  //Outer kij loop
  for (j=0;j<n;j++) {
    for (i=0;i<j;i++) {
      sum=a[i*MAT1+j];
      for (k=0;k<i;k++) 
        sum -= a[i*MAT1+k]*a[k*MAT1+j];
      a[i*MAT1+j]=sum;
    }
    big=0.0;
    //search for largest pivot
    for (i=j;i<n;i++) {
      sum=a[i*MAT1+j];
      for (k=0;k<j;k++) sum -= a[i*MAT1+k]*a[k*MAT1+j];
      a[i*MAT1+j]=sum;
      if ((dum=vv[i]*fabs(sum)) >= big) {
        big=dum;
        imax=i;
      }
    }
    //Do we need to swap any rows?
    if (j != imax) {
      for (k=0;k<n;k++) {
        dum=a[imax*MAT1+k];
        a[imax*MAT1+k]=a[j*MAT1+k];
        a[j*MAT1+k]=dum;
      }
      d = -d;
      vv[imax]=vv[j];
    }
    indx[j]=imax;
    if (a[j*MAT1+j] == 0.0) a[j*MAT1+j]=TINY;
    for (k=j+1;k<n;k++) {
      dum=1.0/(a[j*MAT1+j]);
      for (i=j+1;i<n;i++) a[i*MAT1+j] *= dum;
    }
  }
  return 0;
}

void main(){
    //3x3 Matrix
    float exampleA[]={1,3,-2,3,5,6,2,4,3};
    //pivot array (not used currently)
    int* h_pivot = (int *)malloc(sizeof(int)*MAT1);
    int retval = h_NR_LU_decomp(&exampleA[0],h_pivot);
    for (unsigned int i=0; i<3; i++){
      printf("\n%d:",h_pivot[i]);
      for (unsigned int j=0;j<3; j++){
        printf("%.1lf,",exampleA[i*3+j]);
      }
    }
}

WolframAlpha说答案应该是

1,3,-2
2,-2,7
3,2,-2

我越来越：

2,4,3
0.2,2,-2.8
0.8,1,6.5

到目前为止，我已经找到了至少 3 个不同版本的“相同”算法，所以我完全糊涂了。

PS是的，我知道至少有十几个不同的库可以做到这一点，但我更感兴趣的是了解我做错了什么而不是正确答案。

PPS，因为在LU分解中，较低的结果矩阵是统一的，并且使用（我认为）实现的Crouts算法，数组索引访问仍然是安全的，L和U都可以就地相互叠加；因此为此的单个结果矩阵。

4

4 回答 4

10

我认为您的索引本身就有问题。他们有时有不寻常的开始和结束值，而外循环j而不是i让我怀疑。

在您要求任何人检查您的代码之前，这里有一些建议：

  • 仔细检查你的指数
  • 摆脱那些混淆尝试使用sum
  • 使用宏a(i,j)代替a[i*MAT1+j]
  • 编写子函数而不是注释
  • 删除不必要的部分，隔离错误代码

这是遵循这些建议的版本：

#define MAT1 3
#define a(i,j) a[(i)*MAT1+(j)]

int h_NR_LU_decomp(float *a, int *indx)
{
        int i, j, k;
        int n = MAT1;

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

        return 0;
}

它的主要优点是：

  • 它是可读的
  • 有用

但是，它缺乏旋转。根据需要添加子功能。


我的建议：不要在不理解的情况下复制别人的代码。

大多数程序员都是糟糕的程序员。

于 2011-04-17T01:43:32.797 回答
10

出于对所有神圣事物的热爱，除了作为玩具实现以用于文本中描述的算法的教学目的之外，不要将 Numerical Recipies 代码用于任何东西——而且，实际上，文本并不是那么好。而且，正如您所学习的那样，代码也不是。

当然，不要在您自己的代码中添加任何 Numerical Recipies 例程——许可证非常严格，尤其是考虑到代码质量。如果你有 NR 的东西，你将无法分发你自己的代码。

查看您的系统是否已经安装了LAPACK库。它是计算科学和工程中线性代数例程的标准接口，虽然它并不完美，但您可以为您将代码移动到的任何机器找到 lapack 库，并且您可以编译、链接和运行。如果它尚未安装在您的系统上，您的包管理器（rpm、apt-get、fink、port 等）可能知道 lapack 并可以为您安装它。如果没有，只要您的系统上有 Fortran 编译器，您就可以从这里下载并编译它，标准 C 绑定可以在同一页面的下方找到。

为线性代数例程提供标准 API 如此方便的原因是它们非常常见，但它们的性能非常依赖于系统。因此，例如，Goto BLAS 是 x86 系统的一种非常快速的实现，用于线性代数所需的低级操作；一旦你有 LAPACK 工作，你可以安装那个库来让一切尽可能快。

一旦你安装了任何类型的 LAPACK，对一般矩阵进行 LU 分解的例程是 SGETRF 用于浮点数，或 DGETRF 用于双精度数。如果您对矩阵的结构有所了解，还有其他更快的例程 - 它是对称正定矩阵，例如（SBPTRF），或者它是三对角矩阵（STDTRF）。这是一个很大的库，但是一旦您了解了它的方法，您的数字工具箱中就会拥有一个非常强大的工具。

于 2011-04-17T03:03:36.327 回答
3

对我来说最可疑的是标记为“搜索最大枢轴”的部分。这不仅会搜索，还会改变矩阵 A。我很难相信这是正确的。

LU 算法的不同版本在旋转方面有所不同，因此请确保您理解这一点。您无法比较不同算法的结果。更好的检查是查看 L 乘以 U 是否等于您的原始矩阵，或者如果您的算法确实旋转，则其排列。话虽如此，您的结果是错误的，因为行列式是错误的（旋转不会改变行列式，符号除外）。

除此之外，@Philip 有很好的建议。如果您想了解代码，请先了解 LU 分解而不进行旋转。

于 2011-04-17T10:31:51.353 回答
3

用阿尔伯特爱因斯坦的话说：

...一个有手表的人总是知道确切的时间，但一个有两个的人永远不确定......

您的代码肯定不会产生正确的结果，但即使是这样，带有旋转的结果也不会直接对应于不旋转的结果。在旋转解决方案的上下文中，Alpha 真正给你的可能是这样的：

    1  0  0      1  0  0       1  3 -2
P=  0  1  0   L= 2  1  0  U =  0 -2  7
    0  0  1      3  2  1       0  0 -2

然后它将满足条件 A = PLU（其中 . 表示矩阵乘积）。如果我以另一种方式计算（名义上）相同的分解操作（在这种情况下通过 numpy 使用 LAPACK 例程 dgetrf）：

In [27]: A
Out[27]: 
array([[ 1,  3, -2],
       [ 3,  5,  6],
       [ 2,  4,  3]])

In [28]: import scipy.linalg as la

In [29]: LU,ipivot = la.lu_factor(A)

In [30]: print LU
[[ 3.          5.          6.        ]
 [ 0.33333333  1.33333333 -4.        ]
 [ 0.66666667  0.5         1.        ]]

In [31]: print ipivot
[1 1 2]

在使用 ipivot 进行一点黑魔法之后，我们得到

    0  1  0       1        0    0       3  5       6
P = 0  0  1   L = 0.33333  1    0  U =  0  1.3333 -4
    1  0  0       0.66667  0.5  1       0  0       1

这也满足 A = PLU 。这两种分解都是正确的，但它们是不同的，它们不会对应于 NR 代码的正确运行版本。

所以在你决定你是否有“正确”的答案之前，你真的应该花一点时间来理解你复制的代码实现的实际算法。

于 2011-04-18T09:56:09.640 回答