/******************************************************************************* * Find Eigenvalues and Eigenvectors of a symmetric tridiagonal matrix * * using QL method * * ---------------------------------------------------------------------------- * * SAMPLE RUNS * * * * Example #1: (tridiagonal matrix) * * 1 2 0 0 0 * * 2 4 7 0 0 * * A = 0 7 10 8 0 * * 0 0 8 -0.75 -9 * * 0 0 0 -9 10 * * * * Output file ttql2.lst contains: * * * * EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC * * TRIDIAGONAL MATRIX BY QL METHOD. * * * * N = 5 * * * * Input tridiagonal matrix: * * 1.000000 2.000000 0.000000 0.000000 0.000000 * * 2.000000 4.000000 7.000000 0.000000 0.000000 * * 0.000000 7.000000 10.000000 8.000000 0.000000 * * 0.000000 0.000000 8.000000 -0.750000 -9.000000 * * 0.000000 0.000000 0.000000 -9.000000 10.000000 * * * * Error code: 0 * * * * Eigenvalues: * * 1 -9.15659229 * * 2 -0.78071442 * * 3 2.53046878 * * 4 12.53989066 * * 5 19.11694726 * * * * Eigenvectors (in lines]: * * -0.056408 0.286458 -0.522285 1.000000 0.469812 * * 1.000000 -0.890357 0.322363 0.344649 0.287721 * * 1.000000 0.765234 -0.446362 -0.252815 -0.304616 * * 0.097161 0.560616 0.656182 -0.282210 1.000000 * * 0.051876 0.469920 1.000000 0.728439 -0.719095 * * * * ---------------------------------------------------------------------------- * * Reference: From Numath Library By Tuan Dang Trong in Fortran 77 * * [BIBLI 18]. * * * * C++ Release 1.0 By J-P Moreau, Paris * *******************************************************************************/ #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 main() { N = 5; // 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 main diagonal D[1]=ONE; D[2]=4.0; D[3]=10.0; D[4]=-0.75; D[5]=10.0; // define subdiagonal E[1]=ZERO; E[2]=TWO; E[3]=7.0; E[4]=8.0; E[5]=-9.0; // Initialize matrix Z to unity matrix for (I=1; I<=N; I++) for (J=1; J<=N; J++) if (J==I) Z[I][J] = ONE; else Z[I][J] = ZERO; fp = fopen("ttql2.lst","w"); fprintf(fp,"\n EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC\n"); fprintf(fp," TRIDIAGONAL MATRIX BY QL METHOD.\n\n"); fprintf(fp," N = %d\n\n", N); fprintf(fp," Input tridiagonal matrix:\n"); for (I=1; I<=N; I++) { for (J=1; J<=N; J++) { if (J < I-1) fprintf(fp,"%10.6f", ZERO); if (J == I-1) fprintf(fp,"%10.6f", E[I]); if (J == I+1) fprintf(fp,"%10.6f", E[I+1]); if (J == I) fprintf(fp,"%10.6f", D[I]); if (J>I+1 && J<=N) fprintf(fp,"%10.6f", ZERO); } fprintf(fp,"\n"); } 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 ttql2.lst\n\n"); free(vmblock); } // end of file ttql2.cpp