/*********************************************************************** * Find Eigenvalues and Eigenvectors of a real symmetric matrix by * * using Householder reduction and QL method. * * -------------------------------------------------------------------- * * SAMPLE RUNS: * * * * Example #1: (symmetric matrix) * * * * 4 -2 -1 0 * * A = -2 4 0 -1 * * -1 0 4 -2 * * 0 -1 -2 4 * * * * Output file ttred2.lst contains: * * * * EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC REAL MATRIX * * BY HOUSEHOLDER REDUCTION AND QL METHOD. * * * * N = 4 * * * * Input symmetric matrix: * * 4.000000 -2.000000 -1.000000 0.000000 * * -2.000000 4.000000 0.000000 -1.000000 * * -1.000000 0.000000 4.000000 -2.000000 * * 0.000000 -1.000000 -2.000000 4.000000 * * * * Error code: 0 * * * * Eigenvalues: * * 1 1.00000000 * * 2 3.00000000 * * 3 5.00000000 * * 4 7.00000000 * * * * Eigenvectors (in lines): * * 1.000000 1.000000 1.000000 1.000000 * * 1.000000 1.000000 -1.000000 -1.000000 * * 1.000000 -1.000000 1.000000 -1.000000 * * -1.000000 1.000000 1.000000 -1.000000 * * * * * * Example #2: * * * * 1 2 3 -7 12 * * 2 4 7 3 -1 * * A = 3 7 10 8 4 * * -7 3 8 -0.75 -9 * * 12 -1 4 -9 10 * * * * * * EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC MATRIX * * BY HOUSEHOLDER REDUCTION AND QL METHOD. * * * * N = 5 * * * * Input symmetric matrix: * * 1.000000 2.000000 3.000000 -7.000000 12.000000 * * 2.000000 4.000000 7.000000 3.000000 -1.000000 * * 3.000000 7.000000 10.000000 8.000000 4.000000 * * -7.000000 3.000000 8.000000 -0.750000 -9.000000 * * 12.000000 -1.000000 4.000000 -9.000000 10.000000 * * * * Error code: 0 * * * * Eigenvalues: * * 1 -10.48654517 * * 2 -7.77457973 * * 3 0.46334953 * * 4 18.29182060 * * 5 23.75595477 * * * * Eigenvectors (in lines): * * 0.467172 -0.007947 -0.507827 1.000000 0.264432 * * 1.000000 -0.326100 0.209620 -0.147689 -0.815422 * * 0.208812 1.000000 -0.478027 -0.275000 -0.216915 * * 0.093053 0.594911 1.000000 0.455666 0.050740 * * 0.705820 -0.012677 0.131486 -0.527499 1.000000 * * * * -------------------------------------------------------------------- * * Reference: From Numath Library By Tuan Dang Trong in Fortran 77 * * [BIBLI 18]. * * * * C++ Release 1.1 By J-P Moreau, Paris * * -------------------------------------------------------------------- * * Release 1.1: added results of example #2. * ***********************************************************************/ #include #include REAL *D, *E; REAL **Z; REAL zmax; int I,IERR,J,N; FILE *fp; void *vmblock = NULL; REAL Sign(REAL a, REAL b) { if (b ONE) goto e10; e12: for (I=2; I<=N; I++) E[I-1] = E[I]; E[N] = ZERO; F = ZERO; B = ZERO; for (L=1; L<=N; L++) { J = 0; H = EPS*(ABS(D[L])+ABS(E[L])); if (B < H) B = H; // SEEK SMALLEST ELEMENT OF SUBDIAGONAL for (M=L; M<=N; M++) if (ABS(E[M]) <= B) goto e18; e18: if (M == L) goto e26; // START ITERATION e20: if (J == JM) goto e36; J = J+1; // SHIFT G = D[L]; P = (D[L+1]-G)/(TWO*E[L]); R = SQRT(P*P+ONE); D[L] = E[L]/(P + Sign(R,P)); H = G-D[L]; for (I=L+1; I<=N; I++) D[I] = D[I] - H; F = F + H; // QL TRANSFORMATION P = D[M]; C = ONE; S = ZERO; for (I=M-1; I>=L; I--) { G = C*E[I]; H = C*P; if (ABS(P) >= ABS(E[I])) { C = E[I]/P; R = SQRT(C*C+ONE); E[I+1] = S*P*R; S = C/R; C = ONE/R; } else { C = P/E[I]; R = SQRT(C*C+ONE); E[I+1] = S*E[I]*R; S = ONE/R; C = C*S; } P = C*D[I]-S*G; D[I+1] = H+S*(C*G+S*D[I]); // ELEMENTS OF EIGENVECTORS for (K=1; K<=N; K++) { H = Z[K][I+1]; Z[K][I+1] = S*Z[K][I]+C*H; Z[K][I] = Z[K][I]*C-S*H; } } E[L] = S*P; D[L] = C*P; if (ABS(E[L]) > B) goto e20; // CONVERGENCE e26: D[L] = D[L] + F; } // SORT EIGENVALUES AND EIGENVECTORS // IN ASVENDING ORDER for (L=2; L<=N; L++) { I = L-1; K = I; P = D[I]; for (J=L; J<=N; J++) { if (D[J] >= P) goto e30; K = J; P = D[J]; e30: ;} if (K == I) goto e34; D[K] = D[I]; D[I] = P; for (J=1; J<=N; J++) { P = Z[J][I]; Z[J][I] = Z[J][K]; Z[J][K] = P; } e34: ;} return; // NO CONVERGENCE e36: *IER = L; } //TQL2 void TRED2(double **a, int n, double *d, double *e) { /*----------------------------------------------------------------- HOUSEHOLDER REDUCTION OF MATRIX a TO A TRIDIAGONAL FORM. Algorithm: Martin et al., Num. Math. 11, 181-195, 1968. Ref: Smith et al., Matrix Eigensystem Routines -- EISPACK Guide Springer-Verlag, 1976, pp. 489-494. W H Press et al., Numerical Recipes in C, Cambridge U P, 1988, pp. 373-374. ------------------------------------------------------------------- INPUTS: a (R*8) TABLE (N,N) STORING THE REAL ELEMENTS OF INPUT SYMMETRIC MATRIX n (I4) SIZE OF a OUTPUTS: d (R*8) MAIN DIAGONAL(N) OF THE TRIDIAGONAL MATRIX e (R*8) SUB-DIAGONAL (N) OF THE TRIDIAGONAL MATRIX a (R*8) TABLE (N,N) STORING THE ELEMENTS OF THE TRANSFORMATION MATRIX. '----------------------------------------------------------------*/ 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) e[i] = a[i][l]; else { for (k = 1; k <= l; k++) { a[i][k] /= scale; h += a[i][k] * a[i][k]; } f = a[i][l]; g = (f>0 ? -sqrt(h) : sqrt(h)); e[i] = scale * g; h -= f * g; a[i][l] = f - g; f = 0.0; for (j = 1; j <= l; j++) { a[j][i] = a[i][j]/h; g = 0.0; 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; f += e[j] * a[i][j]; } hh = f / (h + h); for (j = 1; j <= l; j++) { f = a[i][j]; e[j] = g = e[j] - hh * f; for (k = 1; k <= j; k++) a[j][k] -= (f * e[k] + g * a[i][k]); } } } else e[i] = a[i][l]; d[i] = h; } d[1] = 0.0; e[1] = 0.0; for (i = 1; i <= n; i++) { l = i - 1; if (d[i]) { for (j = 1; j <= l; j++) { g = 0.0; for (k = 1; k <= l; k++) 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]; a[i][i] = 1.0; for (j = 1; j <= l; j++) a[j][i] = a[i][j] = 0.0; } } //tred2() void main() { N = 4; // allocate memory for matrix A and vectors B, X, D (index 0 not used) vmblock = vminit(); Z = (REAL **) vmalloc(vmblock, MATRIX, N+1, N+1); D = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0); E = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0); if (! vmcomplete(vmblock)) LogError ("No Memory", 0, __FILE__, __LINE__); // define input symmetric matrix Z[1][1]= 4.0; Z[1][2]=-2.0; Z[1][3]=-1.0; Z[1][4]= 0.0; Z[2][1]=-2.0; Z[2][2]= 4.0; Z[2][3]= 0.0; Z[2][4]=-1.0; Z[3][1]=-1.0; Z[3][2]= 0.0; Z[3][3]= 4.0; Z[3][4]=-2.0; Z[4][1]= 0.0; Z[4][2]=-1.0; Z[4][3]=-2.0; Z[4][4]= 4.0; fp = fopen("ttred2.lst","w"); fprintf(fp,"\n EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC REAL MATRIX\n"); fprintf(fp," BY HOUSEHOLDER REDUCTION AND QL METHOD.\n\n"); fprintf(fp," N = %d\n\n", N); fprintf(fp," Input symmetric matrix:\n"); for (I=1; I<=N; I++) { for (J=1; J<=N; J++) fprintf(fp,"%10.6f", Z[I][J]); fprintf(fp,"\n"); } TRED2(Z,N,D,E); TQL2(N,D,E,Z,&IERR); fprintf(fp,"\n Error code: %d\n\n", IERR); // print eigenvalues fprintf(fp," Eigenvalues:\n"); for (I=1; I<=N; I++) fprintf(fp,"%3d%15.8f\n", I, D[I]); // print normalized eigenvectors (with respect to unity) fprintf(fp,"\n Eigenvectors (in lines):\n"); for (J=1; J<=N; J++) { zmax=Z[1][J]; for (I=2; I<=N; I++) if (ABS(Z[I][J]) > ABS(zmax)) zmax=Z[I][J]; for (I=1; I<=N; I++) fprintf(fp,"%10.6f", Z[I][J]/zmax); fprintf(fp,"\n"); } fclose(fp); printf("\n Results in file ttred2.lst\n\n"); free(vmblock); } // end of file ttred2.cpp