C 中的许多 TQLI 实现是否存在错误? [英] Is there an error in many TQLI implementations in C?

查看：23
本文介绍了C 中的许多 TQLI 实现是否存在错误?的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！

问题描述

在

使用基于 0-index 的数组.

关于 TU Graz 和 Stanford 算法他们只需要以特定格式提供输入数据.

这个例子来自:数字食谱第二版.ANSI C 文件

在这个版本中，tqli 使用基于 1-index 的向量和矩阵.调用 tqli 需要特殊的数据准备，由vector 和 matrix 函数.tqli 函数不直接使用普通的float c[10][10].必须准备好数据:

d=vector(1,NP);e=向量(1，NP)；f=向量(1，NP)；a=矩阵(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-base索引矩阵a.

给出了完整的例子 这里.

#define NP 10#define TINY 1.0e-6诠释主要(无效){整数 i,j,k;浮动 *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=向量(1，NP)；e=向量(1，NP)；f=向量(1，NP)；a=矩阵(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("
实对称矩阵的特征向量
");对于 (i=1;i<=NP;i++) {对于 (j=1;j<=NP;j++) {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
","特征值",i," =",d[i]);printf("%11s %14s %9s
","vector","mtrx*vect.","ratio");对于 (j=1;j<=NP;j++) {if (fabs(a[j][i]) < TINY)printf("%12.6f %12.6f %12s
",a[j][i],f[j],"除以 0");别的printf("%12.6f %12.6f %12.6f
",a[j][i],f[j],f[j]/a[j][i]);}printf("按 ENTER 继续...
");(无效) getchar();}free_matrix(a,1,NP,1,NP);自由向量(f，1，NP)；自由向量(e，1，NP)；自由向量(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("
Eigenvectors for a real symmetric matrix
");
    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
","eigenvalue",i," =",d[i]);
        printf("%11s %14s %9s
","vector","mtrx*vect.","ratio");
        for (j=1;j<=NP;j++) {
            if (fabs(a[j][i]) < TINY)
                printf("%12.6f %12.6f %12s
",
                    a[j][i],f[j],"div. by 0");
            else
                printf("%12.6f %12.6f %12.6f
",
                    a[j][i],f[j],f[j]/a[j][i]);
        }
        printf("Press ENTER to continue...
");
        (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屋！

查看全文
相关文章
其他开发最新文章
热门教程
热门工具
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆