'******************************************************************************** '* 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]. * '* * '* Basic Release 1.0 By J-P Moreau, Paris * '* (www.jpmoreau.fr) * '******************************************************************************** 'PROGRAM TEST_TQL2 DEFDBL A-H, O-Z DEFINT I-N N = 5 DIM Z(N, N), D(N), E(N) 'define main diagonal D(1) = 1#: D(2) = 4#: D(3) = 10#: D(4) = -.75: D(5) = 10# 'define subdiagonal E(1) = 0#: E(2) = 2#: E(3) = 7#: E(4) = 8#: E(5) = -9# 'Initialize matrix Z to unity matrix FOR I = 1 TO N FOR J = 1 TO N IF J = I THEN Z(I, J) = 1# ELSE Z(I, J) = 0# END IF NEXT J NEXT I OPEN "ttql2.lst" FOR OUTPUT AS #2 F$ = "###.######" F1$ = "###.########" CLS PRINT #2, PRINT #2, " EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC" PRINT #2, " TRIDIAGONAL MATRIX BY QL METHOD." PRINT #2, PRINT #2, " N = "; N PRINT #2, PRINT #2, " Input tridiagonal matrix:" FOR I = 1 TO N FOR J = 1 TO N IF J < I - 1 THEN PRINT #2, USING F$; 0#; IF J = I - 1 THEN PRINT #2, USING F$; E(I); IF J = I + 1 THEN PRINT #2, USING F$; E(I + 1); IF J = I THEN PRINT #2, USING F$; D(I); IF J > I + 1 AND J <= N THEN PRINT #2, USING F$; 0#; NEXT J PRINT #2, NEXT I PRINT #2, GOSUB 1000 'call TQL2(N,D,E,Z,IER) PRINT #2, " Error code: "; IER PRINT #2, 'print eigenvalues PRINT #2, " Eigenvalues:" FOR I = 1 TO N PRINT #2, " "; I; PRINT #2, USING F1$; D(I) NEXT I PRINT #2, 'print normalized eigenvectors (with respect to unity) PRINT #2, " Eigenvectors (in lines):" FOR J = 1 TO N zmax = Z(1, J) FOR I = 2 TO N IF ABS(Z(I, J)) > ABS(zmax) THEN zmax = Z(I, J) NEXT I FOR I = 1 TO N PRINT #2, USING F$; Z(I, J) / zmax; NEXT I PRINT #2, NEXT J PRINT #2, CLOSE #2 PRINT PRINT " Results in file ttql2.lst" PRINT END 'of main program 1000 'SUBROUTINE TQL2(N,D,E,Z,IER) '------------------------------------------------------------------------- ' QL METHOD TO DETERMINE THE EIGENVALUES AND EIGENVECTORS OF: ' ' 1) A SYMMETRIC TRIDIAGONAL MATRIX. ' 2) A FULL SYMMETRIC MATRIX AFTER A PREVIOUS CALL TO TRED2. ' ' CALLING MODE: ' CALL TQL2(NM,N,D,E,Z,IER) ' INPUTSS: ' N (I4) SIZE OF Z ' D (R*8) MAIN DIAGONAL (N) OF THE TRIDIAGONAL MATRIX ' E (R*8) SUB-DIAGONAL (N) OF THE TRIDIAGONAL MATRIX ' Z (R*8) TABLE (NM,N) STORING THE UNITY MATRIX IF THE TRIDIAGONAL ' MATRIX IS DEFINED BY D AND E, CASE #1. ' FOR CASE #2, IT CONTAINS THE ELEMENTS OF THE TRANSFORMATION ' MATRIX AFTER A CALL TO TRED2. ' OUTPUTS: ' D (R*8) EIGENVALUES ' Z (R*8) EIGENVECTORS ' IER (I4) ERROR CODE = 0, CONVERGENCE OK. ' = L, NO CONVERGENCE FOR THE Lth EIGENVALUE ' ' REFERENCE: ' J.H.WILKINSON,-C.REINSCH,R.S.MARTIN ' HANDBOOK FOR AUTOMATIC COMPUTATION, VOL.2, LINEAR ALGEBRA ' SPRINGER-VERLAG 1971. '------------------------------------------------------------------------- EPS = 0#: JM = 30 IER = 0 IF N = .1 THEN RETURN ' MACHINE EPSILON IF EPS <> 0# THEN GOTO 12 EPS = 1# 10 EPS = EPS / 2# EPS1 = 1# + EPS IF EPS1 > 1# THEN GOTO 10 12 FOR I = 2 TO N E(I - 1) = E(I) NEXT I E(N) = 0# F = 0# B = 0# FOR L = 1 TO N J = 0 H = EPS * (ABS(D(L)) + ABS(E(L))) IF B < H THEN B = H ' SEEK SMALLEST ELEMENT OF SUBDIAGONAL FOR M = L TO N IF ABS(E(M)) <= B THEN GOTO 18 NEXT M 18 IF M = L THEN GOTO 26 ' START ITERATION 20 IF J = JM THEN GOTO 36 J = J + 1 ' SHIFT G = D(L) P = (D(L + 1) - G) / (2# * E(L)) R = SQR(P * P + 1#) IF P < 0# THEN Sign = -ABS(R) ELSE Sign = ABS(R) END IF D(L) = E(L) / (P + Sign) H = G - D(L) FOR I = L + 1 TO N D(I) = D(I) - H NEXT I F = F + H ' QL TRANSFORMATION P = D(M) C = 1# S = 0# FOR I = M - 1 TO L STEP -1 G = C * E(I) H = C * P IF ABS(P) >= ABS(E(I)) THEN C = E(I) / P R = SQR(C * C + 1#) E(I + 1) = S * P * R S = C / R C = 1# / R ELSE C = P / E(I) R = SQR(C * C + 1#) E(I + 1) = S * E(I) * R S = 1# / R C = C * S END IF P = C * D(I) - S * G D(I + 1) = H + S * (C * G + S * D(I)) ' ELEMENTS OF EIGENVECTORS FOR K = 1 TO N H = Z(K, I + 1) Z(K, I + 1) = S * Z(K, I) + C * H Z(K, I) = Z(K, I) * C - S * H NEXT K NEXT I E(L) = S * P D(L) = C * P IF ABS(E(L)) > B THEN GOTO 20 ' CONVERGENCE 26 D(L) = D(L) + F NEXT L ' SORT EIGENVALUES AND EIGENVECTORS ' IN ASVENDING ORDER FOR L = 2 TO N I = L - 1 K = I P = D(I) FOR J = L TO N IF D(J) >= P THEN GOTO 30 K = J P = D(J) 30 NEXT J IF K = I THEN GOTO 34 D(K) = D(I) D(I) = P FOR J = 1 TO N P = Z(J, I) Z(J, I) = Z(J, K) Z(J, K) = P NEXT J 34 NEXT L RETURN ' NO CONVERGENCE 36 IER = L RETURN 'end of file ttql2.bas