数值配方的LU分解不起作用; 我究竟做错了什么？

Question

我从字面上复制并粘贴了C的数字配方的源代码,用于就地LU矩阵分解,问题是它不起作用.

我确信我做的事情很愚蠢,但我很高兴有人能指出我正确的方向; 我整天都在努力工作,看不出我做错了什么.

POST-ANSWER UPDATE:项目已完成并正在运行.感谢大家的指导.

#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可以在就地彼此叠加(我想); 因此,单个结果矩阵.

Answer 1

我认为你的指数本身存在一些错误.它们有时会有不寻常的起始值和结束值,外部循环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;
}

它的主要优点是:

• 它是可读的
• 有用

但它缺乏旋转.根据需要添加子功能.

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

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

Answer 2

对于所有神圣的爱,不要将数字接收代码用于任何事情,除非作为教学目的的玩具实现文本中描述的算法 - 而且,实际上,文本并不那么好.而且,正如你在学习的那样,代码也不是.

当然不要在你自己的代码中加入任何数字Recipies例程 - 许可证是极其严格的限制,特别是考虑到代码质量.如果你有NR的东西,你将无法分发自己的代码.

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

为线性代数例程提供标准API非常方便的原因是它们非常常见,但它们的性能与系统有关.因此,例如,Goto BLAS 是线性代数所需的低级操作的x86系统的一种非常快速的实现; 一旦你有LAPACK工作,你可以安装该库,以尽可能快地完成所有事情.

一旦安装了任何类型的LAPACK,对一般矩阵进行LU分解的例程是浮点数的SGETRF,或者双重的DGETRF.如果您对矩阵的结构有所了解,那么还有其他更快的例程 - 它是对称的正定,比如(SBPTRF),或者它是三对角(STDTRF).这是一个很大的图书馆,但是一旦你学会了它,你就会在数值工具箱中拥有一个非常强大的装备.