在用户Groo发现的this question的评论中，南加州大学高级计算与模拟合作研究所（University of Southern California）为C所做的this TQLI implementation中，存在一个非常基本的错误，即所有阵列都被视为基于一个阵列。虽然我觉得很奇怪，一个非常著名的机构在其中一个代码中会有这样一个基本的错误，这让我更加困惑，基本上，你可以在网上找到的TQLI算法和相关tred2算法的每一个其他实现都会犯同样的错误。
实例：
TU Graz
Stanford
真的有可能是所有不同的人都犯了同样的错误，还是我遗漏了什么？是否有基于1的C-Was数组版本？
最佳答案
好问题！来自上述源代码的源代码表明，从index1开始对数组进行计算。
阿尔索
/*******************************************************************************
Eigenvalue solvers, tred2 and tqli, from "Numerical Recipes in C" (Cambridge
Univ. Press) by W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. 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
使用基于索引的数组。
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 Files
在本版本中，
0使用基于一个索引的向量和矩阵。调用
tqli需要特殊的数据准备，这些数据由tqli和vector函数。普通的matrix函数不能直接使用。必须准备数据：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];
零基索引矩阵
float c[10][10]用于填充1基索引矩阵tqli。给出了here的完整示例。
#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;
}
结论如下：
真的有可能那些不同的人
错了还是我遗漏了什么？
算法是正确的。拼图中丢失的关键是正确的数据准备。
是否有基于1的C-Was数组版本？
不。
关于c - 在C中的许多TQLI实现中是否有错误？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/48949875/