在C中的许多TQLI实现中是否有错误？ [英] Is there an error in many TQLI implementations in C?
问题描述
在
使用基于索引的数组 0 。
TU Graz和Stanford算法
，它们只需要提供特定格式的输入数据即可。
此示例来自：
数字食谱第二版。 ANSI C文件
在此版本中， tqli 使用基于1索引的向量和矩阵。
调用 tqli 需要特殊的数据准备，这些数据由
vector 和<$ c携带$ c> matrix 函数。 tqli 函数不能直接使用普通的 float c [10] [10] 。必须准备数据：
d = vector（1，NP）;
e = vector（1，NP）;
f = vector（1，NP）;
a = matrix（1，NP，1，NP）; （i = 1; i <= NP; i ++）的
（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;
静态浮点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）; （i = 1; i <= NP; i ++）的
（j = 1; j <= NP; j ++）的
a [i] [j] = c [i-1] [j -1];
tred2（a，NP，d，e）;
tqli（d，e，NP，a）;
printf（实对称矩阵的特征向量 n）； （i = 1; i <= NP; i ++）的
{（j = 1; j <= NP; j ++）的
{{b $ b f [j] = 0.0; （k = 1; k <= NP; k ++）的
f [j] + =（c [j-1] [k-1] * a [k] [i]）;
}
printf（％s％3d％s％10.6f\n，特征值，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.6 f％12s\n，
a [j] [i]，f [j]，除以0）；
else
printf（％12.6f％12.6f％12.6f\n，
a [j] [i]，f [j]，f [j] / a [j ][一世]）;
}
printf（按ENTER键继续... \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）;
返回0；
}
结论如下：
是不是所有这些不同的人都犯了相同的
错误，还是我错过了某些东西？
算法正确。难题中缺少的关键是正确的数据准备。
是否有C版本的数组是基于1的？
否。
In the comments of this question user Groo discovered, that in this TQLI implementation for C done by the Collaboratory For Advanced Computing And Simulations of the University of Southern California there is the very basic mistake that all arrays are treated as if they would be one-based. Although it already seems very strange to me that a very renowned institution would have such a basic mistake in one of there codes it confuses me even more, that basically every other implementation of the TQLI algorithm and the related tred2 algorithm you can find online makes the same mistake.
Examples:
Is it really possible that all those different people made the same mistake or am I missing something? Was there a version of C were arrays were 1-based?
Good question! Source code from above sources indicate that calculations are done on the arrays starting from index 1.
Also
/*******************************************************************************
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
*******************************************************************************/
used by https://www.onlinegdb.com employs 1-based index arrays. See:
/******************************************************************************/
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. */
...
}
You can find newer books were those algorithms are updated:
Numerical Recipes 3rd Edition: The Art of Scientific Computing
By William H. Press
uses 0-index based arrays.
In respect to the TU Graz and Stanford algorithms
they just require supplying input data in the specific format.
This is the example from: Numerical Recipes 2nd ed. ANSI C Files
In this edition the tqli uses vectors and matrices which are 1-index based.
Calling the tqli requires special data preparation which are carries by
vector and matrix functions. The ordinary float c[10][10] is not directly used by tqli function. The data have to be prepared:
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];
Zero based index matrix c[10][10] is used to fill 1-base index matrix a.
Full example is given 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;
}
The conclusion is as follows:
Is it really possible that all those different people made the same mistake or am I missing something?
The algorithms are correct. The missing key in the puzzle was the proper data preparation.
Was there a version of C were arrays were 1-based?
No.
这篇关于在C中的许多TQLI实现中是否有错误？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
