'************************************************************************ '* 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 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: (symmetric matrix) * '* * '* 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]. * '* * '* Basic Release 1.1 By J-P Moreau, Paris * '* (www.jpmoreau.fr) * '* -------------------------------------------------------------------- * '* Release 1.1: added results of example #2. * '************************************************************************ 'PROGRAM TEST_TRED2 DEFDBL A-H, O-Z DEFINT I-N N = 4 DIM Z(N, N), D(N), E(N) 'define input symmetric matrix Z(1,1)= 4#: Z(1,2)=-2#: Z(1,3)=-1#: Z(1,4)= 0# Z(2,1)=-2#: Z(2,2)= 4#: Z(2,3)= 0#: Z(2,4)=-1# Z(3,1)=-1#: Z(3,2)= 0#: Z(3,3)= 4#: Z(3,4)=-2# Z(4,1)= 0#: Z(4,2)=-1#: Z(4,3)=-2#: Z(4,4)= 4# OPEN "ttred2.lst" FOR OUTPUT AS #2 F$ = "###.######" F1$ = "###.########" CLS PRINT #2, PRINT #2, " EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC MATRIX" PRINT #2, " BY HOUSEHOLDER REDUCTION AND QL METHOD." PRINT #2, PRINT #2, " N = "; N PRINT #2, PRINT #2, " Input symmetric matrix:" FOR I = 1 TO N FOR J = 1 TO N PRINT #2, USING F$; Z(I,J); NEXT J PRINT #2, NEXT I PRINT #2, GOSUB 2000 'call TRED2(Z,N,D,E) 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 ttred2.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. ' ' INPUTS: ' 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 (N,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 ASCENDING 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 2000 'tred2(Z, N, D, E) '----------------------------------------------------------------- ' HOUSEHOLDER REDUCTION OF MATRIX Z 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: ' Z (R*8) TABLE (N,N) STORING THE REAL ELEMENTS OF INPUT ' SYMMETRIC MATRIX ' N (I4) SIZE OF Z ' OUTPUTS: ' D (R*8) MAIN DIAGONAL(N) OF THE TRIDIAGONAL MATRIX ' E (R*8) SUB-DIAGONAL (N) OF THE TRIDIAGONAL MATRIX ' Z (R*8) TABLE (N,N) STORING THE ELEMENTS OF THE ' TRANSFORMATION MATRIX. '----------------------------------------------------------------- For i = n To 2 Step -1 L = i - 1 H = 0# Xscale = 0# If L > 1 Then For k = 1 To L Xscale = Xscale + Abs(Z(i, k)) Next k If Abs(Xscale) < 2E-16 Then E(i) = Z(i, L) Else For k = 1 To L Z(i, k) = Z(i, k) / Xscale H = H + Z(i, k) * Z(i, k) Next k F = Z(i, L) If F > 0# Then G = -Sqr(H) Else G = Sqr(H) End If E(i) = Xscale * G H = H - F * G Z(i, L) = F - G F = 0# For j = 1 To L Z(j, i) = Z(i, j) / H G = 0# For k = 1 To j G = G + Z(j, k) * Z(i, k) Next k For k = j + 1 To L G = G + Z(k, j) * Z(i, k) Next k E(j) = G / H F = F + E(j) * Z(i, j) Next j hh = F / (H + H) For j = 1 To L F = Z(i, j) E(j) = E(j) - hh * F G = E(j) For k = 1 To j Z(j, k) = Z(j, k) - F * E(k) - G * Z(i, k) Next k Next j End If Else E(i) = Z(i, L) End If D(i) = H Next i D(1) = 0# E(1) = 0# For i = 1 To n L = i - 1 If D(i) <> 0# Then For j = 1 To L G = 0# For k = 1 To L G = G + Z(i, k) * Z(k, j) Next k For k = 1 To L Z(k, j) = Z(k, j) - G * Z(k, i) Next k Next j End If D(i) = Z(i, i) Z(i, i) = 1# For j = 1 To L Z(j, i) = 0# Z(i, j) = 0# Next j Next i return 'end of file ttred2.bas