从数字配方中分解出LU不能正常工作；我究竟做错了什么? [英] LU Decomposition from Numerical Recipes not working; what am I doing wrong?

问题描述

我已经从提供的C的数字食谱的原位复制并粘贴了原位LU矩阵分解的源代码，问题是它不起作用.

I've literally copied and pasted from the supplied source code for Numerical Recipes for C for in-place LU Matrix Decomposition, problem is its not working.

我确定我做的事很愚蠢，但是希望有人能为此指出正确的方向；我整天都在工作，看不到我在做什么错.

I'm sure I'm doing something stupid but would appreciate anyone being able to point me in the right direction on this; I've been working on its all day and can't see what I'm doing wrong.

应答后更新:该项目为

POST-ANSWER UPDATE: The project is finished and working. Thanks to everyone for their guidance.

#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]);
      }
    }
}

我得到了:

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

到目前为止，我已经找到了至少3个相同"算法的不同版本，所以我完全感到困惑.

And so far I have found at least 3 different versions of the 'same' algorithm, so I'm completely confused.

PS是的，我知道至少有十二种不同的库可以执行此操作，但是比正确的答案，我更想了解自己在做什么错.

PS yes I know there are at least a dozen different libraries to do this, but I'm more interested in understanding what I'm doing wrong than the right answer.

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

PPS since in LU Decomposition the lower resultant matrix is unity, and using Crouts algorithm as (i think) implemented, array index access is still safe, both L and U can be superimposed on each other in-place; hence the single resultant matrix for this.

推荐答案

我认为您的索引固有存在问题.它们有时具有异常的开始和结束值，并且j而不是i的外部循环使我感到怀疑.

I think there's something inherently wrong with your indices. They sometimes have unusual start and end values, and the outer loop over j instead of i makes me suspicious.

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

Before you ask anyone to examine your code, here are a few suggestions:

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

• double-check your indices
• get rid of those obfuscation attempts using sum
• use a macro a(i,j) instead of a[i*MAT1+j]
• write sub-functions instead of comments
• remove unnecessary parts, isolating the erroneous code

以下是遵循以下建议的版本:

Here's a version that follows these suggestions:

#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;
}

它的主要优点是:

• 可读
• 有效

不过，它缺乏枢轴作用.根据需要添加子功能.

It lacks pivoting, though. Add sub-functions as needed.

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

My advice: don't copy someone else's code without understanding it.

大多数程序员都是坏程序员.

Most programmers are bad programmers.

这篇关于从数字配方中分解出LU不能正常工作；我究竟做错了什么?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！