数值配方的LU分解不起作用; 我究竟做错了什么？
我已经从提供的用于就地LU矩阵分解的C的Numerical Recipes的源代码中复制和粘贴，问题是它无法正常工作。
我确信我做的事情很愚蠢，但我很高兴有人能指出我正确的方向; 我一整天都在努力，看不出我做错了什么。
POST-ANSWER UPDATE：项目已完成并正在运行 。 感谢大家的指导。
#include #include #include #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 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= 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算法（我认为）实现，arrays索引访问仍然是安全的，L和U可以相互叠加在一起; 因此，单个结果矩阵。
我认为你的指数本身存在一些错误。 它们有时会有不寻常的开始值和结束值，而外部循环超过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; }
它的主要优点是：
- 它是可读的
- 有用
但它缺乏旋转。 根据需要添加子function。
我的建议：不要在不理解的情况下复制别人的代码。
大多数程序员都是糟糕的程序员。
对于所有神圣的爱，不要将数字接收代码用于任何事情，除非作为教学目的的玩具实现文本中描述的算法 – 而且，实际上，文本并不那么好。 而且，正如您在学习的那样，代码也不是。
当然不要在你自己的代码中加入任何数字Recipies例程 – 许可证是极其严格的限制，特别是考虑到代码质量。 如果你有NR内容，你将无法分发自己的代码。
查看您的系统是否已安装LAPACK库。 它是计算科学和工程中线性代数例程的标准接口，虽然它并不完美，但您可以为任何移动代码的机器找到lapack库，并且可以编译，链接和运行。 如果它尚未安装在您的系统上，您的包管理器（rpm，apt-get，fink，port等）可能知道lapack并可以为您安装它。 如果没有，只要您的系统上有Fortran编译器，就可以从这里下载并编译它，并且可以在同一页面的下面找到标准的C绑定。
为线性代数例程提供标准API非常方便的原因是它们非常常见，但它们的性能与系统有关。 因此，例如， Goto BLAS是线性代数所需的低级操作的x86系统的极其快速的实现; 一旦你有LAPACK工作，你可以安装该库，使尽可能快的一切。
一旦安装了任何类型的LAPACK，对一般矩阵进行LU分解的例程是浮点数的SGETRF，或者双数的DGETRF。 如果您对矩阵的结构有所了解，那么还有其他更快的例程 – 它是对称正定，比如说（SBPTRF），或者它是三对角（STDTRF）。 这是一个很大的图书馆，但是一旦你学会了它，你就会在数值工具箱中拥有一个非常强大的装备。
对我来说最可疑的是标有“寻找最大支点”的部分。 这不仅会搜索，还会更改矩阵A.我发现很难相信这是正确的。
不同版本的LU算法在旋转方面有所不同，因此请确保您理解这一点。 您无法比较不同算法的结果。 更好的检查是查看L乘以U是否等于原始矩阵，或者如果算法进行旋转，则检查其是否排列。 话虽如此，你的结果是错误的，因为行列式是错误的（旋转不会改变行列式，除了符号）。
除了@Philip有很好的建议。 如果您想了解代码，首先要理解LU分解而不需要转动。
要严肃地解释阿尔伯特爱因斯坦：
……一个有手表的男人总是知道确切的时间，但是一个有两个男人的男人永远不会确定……
你的代码肯定没有产生正确的结果，但即使是这样，带有旋转的结果也不会直接对应于没有旋转的结果。 在一个透视解决方案的背景下，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（其中。表示矩阵乘积）。 如果我以另一种方式计算（概念上）相同的分解操作（在这种情况下使用LAPACK例程dgetrf通过numpy）：
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代码的正确运行版本。
因此，在您决定是否有“正确”答案之前，您真的应该花一些时间来理解您复制的代码实现的实际算法。