1      SUBROUTINE TQLRAT(N,D,E2,IERR)
2C***BEGIN PROLOGUE  TQLRAT
3C***DATE WRITTEN   760101   (YYMMDD)
4C***REVISION DATE  830518   (YYMMDD)
5C***CATEGORY NO.  D4A5,D4C2A
6C***KEYWORDS  EIGENVALUES,EIGENVECTORS,EISPACK
7C***AUTHOR  SMITH, B. T., ET AL.
8C***PURPOSE  Computes eigenvalues of symmetric tridiagonal matrix
9C            a rational variant of the QL method.
10C***DESCRIPTION
11C
12C     This subroutine is a translation of the ALGOL procedure TQLRAT,
13C     ALGORITHM 464, COMM. ACM 16, 689(1973) by Reinsch.
14C
15C     This subroutine finds the eigenvalues of a SYMMETRIC
16C     TRIDIAGONAL matrix by the rational QL method.
17C
18C     On Input
19C
20C        N is the order of the matrix.
21C
22C        D contains the diagonal elements of the input matrix.
23C
24C        E2 contains the squares of the subdiagonal elements of the
25C          input matrix in its last N-1 positions.  E2(1) is arbitrary.
26C
27C      On Output
28C
29C        D contains the eigenvalues in ascending order.  If an
30C          error exit is made, the eigenvalues are correct and
31C          ordered for indices 1,2,...IERR-1, but may not be
32C          the smallest eigenvalues.
33C
34C        E2 has been destroyed.
35C
36C        IERR is set to
37C          Zero       for normal return,
38C          J          if the J-th eigenvalue has not been
39C                     determined after 30 iterations.
40C
41C     Calls PYTHAG(A,B) for sqrt(A**2 + B**2).
42C
43C     Questions and comments should be directed to B. S. Garbow,
44C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
45C     ------------------------------------------------------------------
46C***REFERENCES  B. T. SMITH, J. M. BOYLE, J. J. DONGARRA, B. S. GARBOW,
47C                 Y. IKEBE, V. C. KLEMA, C. B. MOLER, *MATRIX EIGEN-
48C                 SYSTEM ROUTINES - EISPACK GUIDE*, SPRINGER-VERLAG,
49C                 1976.
50C***ROUTINES CALLED  PYTHAG
51C***END PROLOGUE  TQLRAT
52C
53      INTEGER I,J,L,M,N,II,L1,MML,IERR
54      REAL D(N),E2(N)
55      REAL B,C,F,G,H,P,R,S,MACHEP
56      REAL PYTHAG
57C
58      DATA MACHEP/1.0E0/
59C***FIRST EXECUTABLE STATEMENT  TQLRAT
60      IF (MACHEP .NE. 1.0E0) GO TO 10
61C
62C   --- This code fails to compute MACHEP correctly on IBM machines. ---
63C   --- Replaced by call to R1MACH on 15 Jun 94 by Ron Boisvert.     ---
64C
65C   05 MACHEP = 0.5E0*MACHEP
66C      IF (1.0E0 + MACHEP .GT. 1.0E0) GO TO 05
67C      MACHEP = 2.0E0*MACHEP
68C
69      MACHEP = R1MACH(4)
70C
71   10 IERR = 0
72      IF (N .EQ. 1) GO TO 1001
73C
74      DO 100 I = 2, N
75  100 E2(I-1) = E2(I)
76C
77      F = 0.0E0
78      B = 0.0E0
79      E2(N) = 0.0E0
80C
81      DO 290 L = 1, N
82         J = 0
83         H = MACHEP * (ABS(D(L)) + SQRT(E2(L)))
84         IF (B .GT. H) GO TO 105
85         B = H
86         C = B * B
87C     .......... LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT ..........
88  105    DO 110 M = L, N
89            IF (E2(M) .LE. C) GO TO 120
90C     .......... E2(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
91C                THROUGH THE BOTTOM OF THE LOOP ..........
92  110    CONTINUE
93C
94  120    IF (M .EQ. L) GO TO 210
95  130    IF (J .EQ. 30) GO TO 1000
96         J = J + 1
97C     .......... FORM SHIFT ..........
98         L1 = L + 1
99         S = SQRT(E2(L))
100         G = D(L)
101         P = (D(L1) - G) / (2.0E0 * S)
102         R = PYTHAG(P,1.0E0)
103         D(L) = S / (P + SIGN(R,P))
104         H = G - D(L)
105C
106         DO 140 I = L1, N
107  140    D(I) = D(I) - H
108C
109         F = F + H
110C     .......... RATIONAL QL TRANSFORMATION ..........
111         G = D(M)
112         IF (G .EQ. 0.0E0) G = B
113         H = G
114         S = 0.0E0
115         MML = M - L
116C     .......... FOR I=M-1 STEP -1 UNTIL L DO -- ..........
117         DO 200 II = 1, MML
118            I = M - II
119            P = G * H
120            R = P + E2(I)
121            E2(I+1) = S * R
122            S = E2(I) / R
123            D(I+1) = H + S * (H + D(I))
124            G = D(I) - E2(I) / G
125            IF (G .EQ. 0.0E0) G = B
126            H = G * P / R
127  200    CONTINUE
128C
129         E2(L) = S * G
130         D(L) = H
131C     .......... GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST ..........
132         IF (H .EQ. 0.0E0) GO TO 210
133         IF (ABS(E2(L)) .LE. ABS(C/H)) GO TO 210
134         E2(L) = H * E2(L)
135         IF (E2(L) .NE. 0.0E0) GO TO 130
136  210    P = D(L) + F
137C     .......... ORDER EIGENVALUES ..........
138         IF (L .EQ. 1) GO TO 250
139C     .......... FOR I=L STEP -1 UNTIL 2 DO -- ..........
140         DO 230 II = 2, L
141            I = L + 2 - II
142            IF (P .GE. D(I-1)) GO TO 270
143            D(I) = D(I-1)
144  230    CONTINUE
145C
146  250    I = 1
147  270    D(I) = P
148  290 CONTINUE
149C
150      GO TO 1001
151C     .......... SET ERROR -- NO CONVERGENCE TO AN
152C                EIGENVALUE AFTER 30 ITERATIONS ..........
153 1000 IERR = L
154 1001 RETURN
155      END
156