使用TQLI算法的特征值计算失败，出现分段错误
我试图使用我从南加州大学CACS网站获得的TQLI算法计算特征值。 我的测试脚本如下所示：
#include int main() { int i; i = rand(); printf("My random number: %d\n", i); float d[4] = { {1, 2, 3, 4} }; float e[4] = { {0, 0, 0, 0} }; float z[4][4] = { {1.0, 0.0, 0.0, 0.0} , {0.0, 1.0, 0.0, 0.0} , {0.0, 0.0, 1.0, 0.0}, {0.0, 0.0, 0.0, 1.0} }; double *zptr; zptr = &z[0][0]; printf("Element [2][1] of identity matrix: %f\n", z[2][1]); printf("Element [2][2] of identity matrix: %f\n", z[2][2]); tqli(d, e, 4, zptr); printf("First eigenvalue: %f\n", d[0]); return 0; } 当我尝试运行此脚本时，我会收到分段错误错误，您可以在此处看到。 我的代码在什么位置产生这种分段错误。 我相信来自USC的代码没有错误，我很确定错误必须在我对函数的调用中。 但是我无法看到我在arrays设置中犯了什么错误，因为我认为我遵循了说明。
使用TQLI算法的特征值计算失败，出现分段错误
分段错误来自穿过提供的arrays边界。 tqli需要特定的数据准备。
1） 来自CACS的特征码是基于Fortran的，并且从1开始计算索引。
2） tqli期望其矩阵和double量的double指针。
/******************************************************************************/ void tqli(double d[], double e[], int n, double **z) /*******************************************************************************
d ， e应声明为double 。
3）程序需要修改上述function的数据准备。
必须创建基于助手1索引的向量，以便为tqli提供格式正确的数据：
double z[NP][NP] = { {2, 0, 0}, {0, 4, 0}, {0, 0, 2} } ; double **a; double *d,*e,*f; d=dvector(1,NP); // 1-index based vector e=dvector(1,NP); f=dvector(1,NP); a=dmatrix(1,NP,1,NP); // 1-index based matrix for (i=1;i<=NP;i++) // loading data from zero besed `ze` to `a` for (j=1;j<=NP;j++) a[i][j]=z[i-1][j-1];
完整的测试程序如下。 它使用来自CACS的特征码 ：
/******************************************************************************* Eigenvalue solvers, tred2 and tqli, from "Numerical Recipes in C" (Cambridge Univ. Press) by WH Press, SA Teukolsky, WT Vetterling, and BP Flannery *******************************************************************************/ #include #include #include #define NR_END 1 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) double **dmatrix(int nrl, int nrh, int ncl, int nch) /* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */ { int i,nrow=nrh-nrl+1,ncol=nch-ncl+1; double **m; /* allocate pointers to rows */ m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*))); m += NR_END; m -= nrl; /* allocate rows and set pointers to them */ m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double))); m[nrl] += NR_END; m[nrl] -= ncl; for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; /* return pointer to array of pointers to rows */ return m; } double *dvector(int nl, int nh) /* allocate a double vector with subscript range v[nl..nh] */ { double *v; v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double))); return v-nl+NR_END; } /******************************************************************************/ void tred2(double **a, int n, double d[], double e[]) /******************************************************************************* Householder reduction of a real, symmetric matrix a[1..n][1..n]. On output, a is replaced by the orthogonal matrix Q effecting the transformation. d[1..n] returns the diagonal elements of the tridiagonal matrix, and e[1..n] the off-diagonal elements, with e[1]=0. Several statements, as noted in comments, can be omitted if only eigenvalues are to be found, in which case a contains no useful information on output. Otherwise they are to be included. *******************************************************************************/ { int l,k,j,i; double scale,hh,h,g,f; for (i=n;i>=2;i--) { l=i-1; h=scale=0.0; if (l > 1) { for (k=1;k<=l;k++) scale += fabs(a[i][k]); if (scale == 0.0) /* Skip transformation. */ e[i]=a[i][l]; else { for (k=1;k<=l;k++) { a[i][k] /= scale; /* Use scaled a's for transformation. */ h += a[i][k]*a[i][k]; /* Form sigma in h. */ } f=a[i][l]; g=(f >= 0.0 ? -sqrt(h) : sqrt(h)); e[i]=scale*g; h -= f*g; /* Now h is equation (11.2.4). */ a[i][l]=fg; /* Store u in the ith row of a. */ f=0.0; for (j=1;j<=l;j++) { /* Next statement can be omitted if eigenvectors not wanted */ a[j][i]=a[i][j]/h; /* Store u/H in ith column of a. */ g=0.0; /* Form an element of Au in g. */ for (k=1;k<=j;k++) g += a[j][k]*a[i][k]; for (k=j+1;k<=l;k++) g += a[k][j]*a[i][k]; e[j]=g/h; /* Form element of p in temporarily unused element of e. */ f += e[j]*a[i][j]; } hh=f/(h+h); /* Form K, equation (11.2.11). */ for (j=1;j<=l;j++) { /* Form q and store in e overwriting p. */ f=a[i][j]; e[j]=g=e[j]-hh*f; for (k=1;k<=j;k++) /* Reduce a, equation (11.2.13). */ a[j][k] -= (f*e[k]+g*a[i][k]); } } } else e[i]=a[i][l]; d[i]=h; } /* Next statement can be omitted if eigenvectors not wanted */ d[1]=0.0; e[1]=0.0; /* Contents of this loop can be omitted if eigenvectors not wanted except for statement d[i]=a[i][i]; */ for (i=1;i<=n;i++) { /* Begin accumulation of transformation matrices. */ l=i-1; if (d[i]) { /* This block skipped when i=1. */ for (j=1;j<=l;j++) { g=0.0; for (k=1;k<=l;k++) /* Use u and u/H stored in a to form PQ */ g += a[i][k]*a[k][j]; for (k=1;k<=l;k++) a[k][j] -= g*a[k][i]; } } d[i]=a[i][i]; /* This statement remains. */ a[i][i]=1.0; /* Reset row and column of a to identity matrix for next iteration. */ for (j=1;j<=l;j++) a[j][i]=a[i][j]=0.0; } } /******************************************************************************/ 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. */ e[n]=0.0; for (l=1;l<=n;l++) { iter=0; do { for (m=l;m<=n-1;m++) { /* Look for a single small subdiagonal element to split the matrix. */ dd=fabs(d[m])+fabs(d[m+1]); if ((double)(fabs(e[m])+dd) == dd) break; } if (m != l) { if (iter++ == 30) printf("Too many iterations in tqli"); g=(d[l+1]-d[l])/(2.0*e[l]); /* Form shift. */ r=pythag(g,1.0); g=d[m]-d[l]+e[l]/(g+SIGN(r,g)); /* This is dm - ks. */ s=c=1.0; p=0.0; for (i=m-1;i>=l;i--) { /* A plane rotation as in the original QL, followed by Givens */ f=s*e[i]; /* rotations to restore tridiagonal form. */ b=c*e[i]; e[i+1]=(r=pythag(f,g)); if (r == 0.0) { /* Recover from underflow. */ d[i+1] -= p; e[m]=0.0; break; } s=f/r; c=g/r; g=d[i+1]-p; r=(d[i]-g)*s+2.0*c*b; d[i+1]=g+(p=s*r); g=c*rb; /* Next loop can be omitted if eigenvectors not wanted */ for (k=1;k<=n;k++) { /* Form eigenvectors. */ f=z[k][i+1]; z[k][i+1]=s*z[k][i]+c*f; z[k][i]=c*z[k][i]-s*f; } } if (r == 0.0 && i >= l) continue; d[l] -= p; e[l]=g; e[m]=0.0; } } while (m != l); } } /******************************************************************************/ double pythag(double a, double b) /******************************************************************************* Computes (a2 + b2)1/2 without destructive underflow or overflow. *******************************************************************************/ { double absa,absb; absa=fabs(a); absb=fabs(b); if (absa > absb) return absa*sqrt(1.0+(absb/absa)*(absb/absa)); else return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+(absa/absb)*(absa/absb))); } #define NP 3 #define TINY 1.0e-6 double sqrt(double x) { union { int i; double x; } u; ux = x; ui = (1<<29) + (ui >> 1) - (1<<22); return ux; } int main() { int i,j,k; double ze[NP][NP] = { {2, 0, 0}, {0, 4, 0}, {0, 0, 2} } ; double **a; double *d,*e,*f; d=dvector(1,NP); e=dvector(1,NP); f=dvector(1,NP); a=dmatrix(1,NP,1,NP); for (i=1;i<=NP;i++) for (j=1;j<=NP;j++) a[i][j]=ze[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] += (ze[j-1][k-1]*a[k][i]); } printf("%s %3d %s %10.6f\n","\neigenvalue",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]); } } //free_dmatrix(a,1,NP,1,NP); //free_dvector(f,1,NP); //free_dvector(e,1,NP); //free_dvector(d,1,NP); return 0; }
输出：
Eigenvectors for a real symmetric matrix: eigenvalue 1 = 2.000000 vector mtrx*vect. ratio 1.000000 2.000000 2.000000 0.000000 0.000000 div. by 0 0.000000 0.000000 div. by 0 eigenvalue 2 = 4.000000 vector mtrx*vect. ratio 0.000000 0.000000 div. by 0 1.000000 4.000000 4.000000 0.000000 0.000000 div. by 0 eigenvalue 3 = 2.000000 vector mtrx*vect. ratio 0.000000 0.000000 div. by 0 0.000000 0.000000 div. by 0 1.000000 2.000000 2.000000
我希望它最终有助于澄清有关tqli数据准备的tqli 。