1      SUBROUTINE IMTQL2(NM,N,D,E,Z,IERR)
2C***BEGIN PROLOGUE  IMTQL2
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 and eigenvectors of symmetric
9C            tridiagonal matrix using implicit QL method.
10C***DESCRIPTION
11C
12C     This subroutine is a translation of the ALGOL procedure IMTQL2,
13C     NUM. MATH. 12, 377-383(1968) by Martin and Wilkinson,
14C     as modified in NUM. MATH. 15, 450(1970) by Dubrulle.
15C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971).
16C
17C     This subroutine finds the eigenvalues and eigenvectors
18C     of a SYMMETRIC TRIDIAGONAL matrix by the implicit QL method.
19C     The eigenvectors of a FULL SYMMETRIC matrix can also
20C     be found if  TRED2  has been used to reduce this
21C     full matrix to tridiagonal form.
22C
23C     On INPUT
24C
25C        NM must be set to the row dimension of two-dimensional
26C          array parameters as declared in the calling program
27C          dimension statement.
28C
29C        N is the order of the matrix.
30C
31C        D contains the diagonal elements of the input matrix.
32C
33C        E contains the subdiagonal elements of the input matrix
34C          in its last N-1 positions.  E(1) is arbitrary.
35C
36C        Z contains the transformation matrix produced in the
37C          reduction by  TRED2, if performed.  If the eigenvectors
38C          of the tridiagonal matrix are desired, Z must contain
39C          the identity matrix.
40C
41C      On OUTPUT
42C
43C        D contains the eigenvalues in ASCENDING order.  If an
44C          error exit is made, the eigenvalues are correct but
45C          UNORDERED for indices 1,2,...,IERR-1.
46C
47C        E has been destroyed.
48C
49C        Z contains orthonormal eigenvectors of the symmetric
50C          tridiagonal (or full) matrix.  If an error exit is made,
51C          Z contains the eigenvectors associated with the stored
52C          eigenvalues.
53C
54C        IERR is set to
55C          ZERO       for normal return,
56C          J          if the J-th eigenvalue has not been
57C                     determined after 30 iterations.
58C
59C     Calls PYTHAG(A,B) for sqrt(A**2 + B**2).
60C
61C     Questions and comments should be directed to B. S. Garbow,
62C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
63C     ------------------------------------------------------------------
64C***REFERENCES  B. T. SMITH, J. M. BOYLE, J. J. DONGARRA, B. S. GARBOW,
65C                 Y. IKEBE, V. C. KLEMA, C. B. MOLER, *MATRIX EIGEN-
66C                 SYSTEM ROUTINES - EISPACK GUIDE*, SPRINGER-VERLAG,
67C                 1976.
68C***ROUTINES CALLED  PYTHAG
69C***END PROLOGUE  IMTQL2
70C
71      INTEGER I,J,K,L,M,N,II,NM,MML,IERR
72      REAL D(N),E(N),Z(NM,N)
73      REAL B,C,F,G,P,R,S,S1,S2
74      REAL PYTHAG
75C
76C***FIRST EXECUTABLE STATEMENT  IMTQL2
77      IERR = 0
78      IF (N .EQ. 1) GO TO 1001
79C
80      DO 100 I = 2, N
81  100 E(I-1) = E(I)
82C
83      E(N) = 0.0E0
84C
85      DO 240 L = 1, N
86         J = 0
87C     .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT ..........
88  105    DO 110 M = L, N
89            IF (M .EQ. N) GO TO 120
90            S1 = ABS(D(M)) + ABS(D(M+1))
91            S2 = S1 + ABS(E(M))
92            IF (S2 .EQ. S1) GO TO 120
93  110    CONTINUE
94C
95  120    P = D(L)
96         IF (M .EQ. L) GO TO 240
97         IF (J .EQ. 30) GO TO 1000
98         J = J + 1
99C     .......... FORM SHIFT ..........
100         G = (D(L+1) - P) / (2.0E0 * E(L))
101         R = PYTHAG(G,1.0E0)
102         G = D(M) - P + E(L) / (G + SIGN(R,G))
103         S = 1.0E0
104         C = 1.0E0
105         P = 0.0E0
106         MML = M - L
107C     .......... FOR I=M-1 STEP -1 UNTIL L DO -- ..........
108         DO 200 II = 1, MML
109            I = M - II
110            F = S * E(I)
111            B = C * E(I)
112            IF (ABS(F) .LT. ABS(G)) GO TO 150
113            C = G / F
114            R = SQRT(C*C+1.0E0)
115            E(I+1) = F * R
116            S = 1.0E0 / R
117            C = C * S
118            GO TO 160
119  150       S = F / G
120            R = SQRT(S*S+1.0E0)
121            E(I+1) = G * R
122            C = 1.0E0 / R
123            S = S * C
124  160       G = D(I+1) - P
125            R = (D(I) - G) * S + 2.0E0 * C * B
126            P = S * R
127            D(I+1) = G + P
128            G = C * R - B
129C     .......... FORM VECTOR ..........
130            DO 180 K = 1, N
131               F = Z(K,I+1)
132               Z(K,I+1) = S * Z(K,I) + C * F
133               Z(K,I) = C * Z(K,I) - S * F
134  180       CONTINUE
135C
136  200    CONTINUE
137C
138         D(L) = D(L) - P
139         E(L) = G
140         E(M) = 0.0E0
141         GO TO 105
142  240 CONTINUE
143C     .......... ORDER EIGENVALUES AND EIGENVECTORS ..........
144      DO 300 II = 2, N
145         I = II - 1
146         K = I
147         P = D(I)
148C
149         DO 260 J = II, N
150            IF (D(J) .GE. P) GO TO 260
151            K = J
152            P = D(J)
153  260    CONTINUE
154C
155         IF (K .EQ. I) GO TO 300
156         D(K) = D(I)
157         D(I) = P
158C
159         DO 280 J = 1, N
160            P = Z(J,I)
161            Z(J,I) = Z(J,K)
162            Z(J,K) = P
163  280    CONTINUE
164C
165  300 CONTINUE
166C
167      GO TO 1001
168C     .......... SET ERROR -- NO CONVERGENCE TO AN
169C                EIGENVALUE AFTER 30 ITERATIONS ..........
170 1000 IERR = L
171 1001 RETURN
172      END
173