**** for old version, "send otqlrat from seispack"
** From dana!moler Tue, 1 Sep 87 10:15:40 PDT
** New TQLRAT
SUBROUTINE TQLRAT(N,D,E2,IERR)
C
INTEGER I,J,L,M,N,II,L1,MML,IERR
REAL D(N),E2(N)
REAL B,C,F,G,H,P,R,S,T,EPSLON,PYTHAG
C
C This subroutine is a translation of the Algol procedure tqlrat,
C Algorithm 464, Comm. ACM 16, 689(1973) by Reinsch.
C
C This subroutine finds the eigenvalues of a symmetric
C tridiagonal matrix by the rational QL method.
C
C On input
C
C N is the order of the matrix.
C
C D contains the diagonal elements of the input matrix.
C
C E2 contains the squares of the subdiagonal elements of the
C input matrix in its last N-1 positions. E2(1) is arbitrary.
C
C On output
C
C D contains the eigenvalues in ascending order. If an
C error exit is made, the eigenvalues are correct and
C ordered for indices 1,2,...IERR-1, but may not be
C the smallest eigenvalues.
C
C E2 has been destroyed.
C
C IERR is set to
C zero for normal return,
C J if the J-th eigenvalue has not been
C determined after 30 iterations.
C
C Calls PYTHAG for SQRT(A*A + B*B) .
C
C Questions and comments should be directed to Burton S. Garbow,
C Mathematics and Computer Science Div, Argonne National Laboratory
C
C This version dated August 1987.
C Modified by C. Moler to fix underflow/overflow difficulties,
C especially on the VAX and other machines where epslon(1.0e0)**2
C nearly underflows. See the loop involving statement 102 and
C the two statements just before statement 200.
C
C ------------------------------------------------------------------
C
IERR = 0
IF (N .EQ. 1) GO TO 1001
C
DO 100 I = 2, N
100 E2(I-1) = E2(I)
C
F = 0.0E0
T = 0.0E0
E2(N) = 0.0E0
C
DO 290 L = 1, N
J = 0
H = ABS(D(L)) + SQRT(E2(L))
IF (T .GT. H) GO TO 105
T = H
B = EPSLON(T)
C = B * B
if (c .ne. 0.0e0) go to 105
C Spliting tolerance underflowed. Look for larger value.
do 102 i = l, n
h = abs(d(i)) + sqrt(e2(i))
if (h .gt. t) t = h
102 continue
b = epslon(t)
c = b * b
C .......... LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT ..........
105 DO 110 M = L, N
IF (E2(M) .LE. C) GO TO 120
C .......... E2(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
C THROUGH THE BOTTOM OF THE LOOP ..........
110 CONTINUE
C
120 IF (M .EQ. L) GO TO 210
130 IF (J .EQ. 30) GO TO 1000
J = J + 1
C .......... FORM SHIFT ..........
L1 = L + 1
S = SQRT(E2(L))
G = D(L)
P = (D(L1) - G) / (2.0E0 * S)
R = PYTHAG(P,1.0E0)
D(L) = S / (P + SIGN(R,P))
H = G - D(L)
C
DO 140 I = L1, N
140 D(I) = D(I) - H
C
F = F + H
C .......... RATIONAL QL TRANSFORMATION ..........
G = D(M)
IF (G .EQ. 0.0E0) G = B
H = G
S = 0.0E0
MML = M - L
C .......... FOR I=M-1 STEP -1 UNTIL L DO -- ..........
DO 200 II = 1, MML
I = M - II
P = G * H
R = P + E2(I)
E2(I+1) = S * R
S = E2(I) / R
D(I+1) = H + S * (H + D(I))
G = D(I) - E2(I) / G
C Avoid division by zero on next pass
if (g .eq. 0.0e0) g = epslon(d(i))
h = g * (p / r)
200 CONTINUE
C
E2(L) = S * G
D(L) = H
C .......... GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST ..........
IF (H .EQ. 0.0E0) GO TO 210
IF (ABS(E2(L)) .LE. ABS(C/H)) GO TO 210
E2(L) = H * E2(L)
IF (E2(L) .NE. 0.0E0) GO TO 130
210 P = D(L) + F
C .......... ORDER EIGENVALUES ..........
IF (L .EQ. 1) GO TO 250
C .......... FOR I=L STEP -1 UNTIL 2 DO -- ..........
DO 230 II = 2, L
I = L + 2 - II
IF (P .GE. D(I-1)) GO TO 270
D(I) = D(I-1)
230 CONTINUE
C
250 I = 1
270 D(I) = P
290 CONTINUE
C
GO TO 1001
C .......... SET ERROR -- NO CONVERGENCE TO AN
C EIGENVALUE AFTER 30 ITERATIONS ..........
1000 IERR = L
1001 RETURN
END