数值配方的LU分解不起作用;我究竟做错了什么?
我从字面上复制并粘贴了C的数字配方的源代码,用于就地LU矩阵分解,问题是它不起作用.
我确信我做的事情很愚蠢,但我很高兴有人能指出我正确的方向;我整天都在努力工作,看不出我做错了什么. 后答案更新:项目是finished and working.感谢大家的指导. #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,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,-2 2,7 3,-2 我越来越: 2,3 0.2,-2.8 0.8,1,6.5 到目前为止,我已经找到了至少3个不同版本的“相同”算法,所以我完全糊涂了. PS是的我知道至少有十几个不同的库可以做到这一点,但我更感兴趣的是理解我做错了什么而不是正确的答案. PPS因为在LU分解中,较低的结果矩阵是单位的,并且使用Crouts算法(我认为)实现,阵列索引访问仍然是安全的,L和U可以在原地相互叠加;因此,单个结果矩阵. 解决方法
我认为你的指数本身存在一些错误.它们有时会有不寻常的开始值和结束值,而外部循环超过j而不是我让我怀疑.
在你要求任何人检查你的代码之前,这里有一些建议: >仔细检查您的指数 这是一个遵循以下建议的版本: #define MAT1 3 #define a(i,j) a[(i)*MAT1+(j)] int h_NR_LU_decomp(float *a,int *indx) { int i,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,i); } return 0; } 它的主要优点是: >它的可读性 但它缺乏旋转.根据需要添加子功能. 我的建议:不要在不理解的情况下复制别人的代码. 大多数程序员都是糟糕的程序员. (编辑:李大同) 【声明】本站内容均来自网络,其相关言论仅代表作者个人观点,不代表本站立场。若无意侵犯到您的权利,请及时与联系站长删除相关内容! |