在C中的许多TQLI实现中是否存在错误？
用户Groo在这个问题的评论中发现，在南加利福尼亚大学高级计算和模拟合作实验室完成的C的TQLI实现中 ，存在一个非常基本的错误，即所有arrays都被视为一个 -根据。 虽然对我来说已经很奇怪，一个非常着名的机构会在其中一个代码中出现这样一个基本错误，但它让我更加困惑，基本上每个其他的TQLI算法实现和你可以在网上找到的相关tred2算法都能让同样的错误。
例子：
- TU Graz
- 斯坦福
这些不同的人真的有可能犯同样的错误，或者我错过了什么？ 是否有一个版本的Carrays是基于1？
好问题！ 来自上述源的源代码表明从索引1开始对数组进行计算。 也
/******************************************************************************* Eigenvalue solvers, tred2 and tqli, from "Numerical Recipes in C" (Cambridge Univ. Press) by WH Press, SA Teukolsky, WT Vetterling, and BP Flannery *******************************************************************************/
由https://www.onlinegdb.com使用，使用基于1的索引数组。 看到：
/******************************************************************************/ void tqli(double d[], double e[], int n, double **z) /******************************************************************************* QL algorithm with implicit shifts, to determine the eigenvalues and eigenvectors of a real, symmetric, tridiagonal matrix, or of a real, symmetric matrix previously reduced by tred2 sec. 11.2. On input, d[1..n] contains the diagonal elements of the tridiagonal matrix. On output, it returns the eigenvalues. The vector e[1..n] inputs the subdiagonal elements of the tridiagonal matrix, with e[1] arbitrary. On output e is destroyed. When finding only the eigenvalues, several lines may be omitted, as noted in the comments. If the eigenvectors of a tridiagonal matrix are desired, the matrix z[1..n][1..n] is input as the identity matrix. If the eigenvectors of a matrix that has been reduced by tred2 are required, then z is input as the matrix output by tred2. In either case, the kth column of z returns the normalized eigenvector corresponding to d[k]. *******************************************************************************/ { double pythag(double a, double b); int m,l,iter,i,k; double s,r,p,g,f,dd,c,b; for (i=2;i<=n;i++) e[i-1]=e[i]; /* Convenient to renumber the elements of e. */ ... }
您可以找到更新的那些算法的新书：
Numerical Recipes 3rd Edition: The Art of Scientific Computing By William H. Press
使用基于0索引的数组。
In respect to the TU Graz and Stanford algorithms they just require supplying input data in the specific format.
这是一个例子： Numerical Recipes 2nd ed。 ANSI C文件
在这个版本中， tqli使用基于1索引的向量和矩阵。 调用tqli需要特殊的数据准备，这些数据准备由vector和matrix函数携带。 普通float c[10][10]不是由tqli函数直接使用的。 必须准备数据：
d=vector(1,NP); e=vector(1,NP); f=vector(1,NP); a=matrix(1,NP,1,NP); for (i=1;i<=NP;i++) for (j=1;j<=NP;j++) a[i][j]=c[i-1][j-1];
基于零的索引矩阵c[10][10]用于填充1-基索引矩阵a 。
这里给出了完整的例子。
#define NP 10 #define TINY 1.0e-6 int main(void) { int i,j,k; float *d,*e,*f,**a; static float c[NP][NP]={ 5.0, 4.3, 3.0, 2.0, 1.0, 0.0,-1.0,-2.0,-3.0,-4.0, 4.3, 5.1, 4.0, 3.0, 2.0, 1.0, 0.0,-1.0,-2.0,-3.0, 3.0, 4.0, 5.0, 4.0, 3.0, 2.0, 1.0, 0.0,-1.0,-2.0, 2.0, 3.0, 4.0, 5.0, 4.0, 3.0, 2.0, 1.0, 0.0,-1.0, 1.0, 2.0, 3.0, 4.0, 5.0, 4.0, 3.0, 2.0, 1.0, 0.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 4.0, 3.0, 2.0, 1.0, -1.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 4.0, 3.0, 2.0, -2.0,-1.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 4.0, 3.0, -3.0,-2.0,-1.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 4.0, -4.0,-3.0,-2.0,-1.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0}; d=vector(1,NP); e=vector(1,NP); f=vector(1,NP); a=matrix(1,NP,1,NP); for (i=1;i<=NP;i++) for (j=1;j<=NP;j++) a[i][j]=c[i-1][j-1]; tred2(a,NP,d,e); tqli(d,e,NP,a); printf("\nEigenvectors for a real symmetric matrix\n"); for (i=1;i<=NP;i++) { for (j=1;j<=NP;j++) { f[j]=0.0; for (k=1;k<=NP;k++) f[j] += (c[j-1][k-1]*a[k][i]); } printf("%s %3d %s %10.6f\n","eigenvalue",i," =",d[i]); printf("%11s %14s %9s\n","vector","mtrx*vect.","ratio"); for (j=1;j<=NP;j++) { if (fabs(a[j][i]) < TINY) printf("%12.6f %12.6f %12s\n", a[j][i],f[j],"div. by 0"); else printf("%12.6f %12.6f %12.6f\n", a[j][i],f[j],f[j]/a[j][i]); } printf("Press ENTER to continue...\n"); (void) getchar(); } free_matrix(a,1,NP,1,NP); free_vector(f,1,NP); free_vector(e,1,NP); free_vector(d,1,NP); return 0; }
结论如下：
这些不同的人真的有可能犯同样的错误，或者我错过了什么？
算法是正确的。 这个难题中缺少的关键是正确的数据准备。
是否有一个版本的Carrays是基于1？
没有。