      SUBROUTINE IMTQL2 (NM, N, D, E, Z, IERR, JOB)
      IMPLICIT NONE
C
C A TRANSLATION OF THE ALGOL PROCEDURE IMTQL2, NUM. MATH. 12, 377-383 (1968)
C   BY MARTIN AND WILKINSON, AS MODIFIED IN NUM. MATH. 15, 450 (1970) BY
C   DUBRULLE.  HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248 (1971).
C
C FINDS THE EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC TRIDIAGONAL MATRIX
C   BY THE IMPLICIT QL METHOD.  THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX
C   CAN ALSO BE FOUND IF TRED2 HAS BEEN USED TO REDUCE THIS FULL MATRIX TO
C   TRIDIAGONAL FORM.
C
C ON INPUT:
C   NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS
C      AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT.
C
C   N IS THE ORDER OF THE MATRIX.
C
C   D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
C
C   E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX IN ITS LAST
C      N-1 POSITIONS.  E(1) IS ARBITRARY.
C
C   Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE REDUCTION BY TRED2,
C      IF PERFORMED.  IF THE EIGENVECTORS OF THE TRIDIAGONAL MATRIX ARE
C      DESIRED, Z MUST CONTAIN THE IDENTITY MATRIX.
C
C ON OUTPUT:
C   D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN ERROR EXIT IS MADE,
C      THE EIGENVALUES ARE CORRECT BUT UNORDERED FOR INDICES 1,2,...,IERR-1.
C
C   E HAS BEEN DESTROYED.
C
C   Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC TRIDIAGONAL (OR
C      FULL) MATRIX.  IF AN ERROR EXIT IS MADE, Z CONTAINS THE EIGENVECTORS
C      ASSOCIATED WITH THE STORED EIGENVALUES.
C
C   IERR = ZERO FOR NORMAL RETURN,
C          J    IF J-TH EIGENVALUE HAS NOT BEEN DETERMINED AFTER 30 ITERATIONS.
C
C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C   APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
      INTEGER NM, N, IERR, JOB
      DOUBLE PRECISION D(N), E(N), Z(NM,N)
C
      INTEGER I, J, K, L, M, II, MML
      DOUBLE PRECISION B, C, F, G, P, R, S
C
      DOUBLE PRECISION FLOP
C
C
      IERR = 0
      IF (N.EQ.1) GO TO 1001
      DO 100 I = 2, N
        E(I-1) = E(I)
100   CONTINUE
C
      E(N) = 0.0D0
      DO 240 L = 1, N
        J = 0
C
C ***      LOOK FOR SMALL SUB-DIAGONAL ELEMENT
105     CONTINUE
        DO 110 M = L, N
          IF (M.EQ.N) GO TO 120
          P = FLOP (DABS (D(M))+DABS (D(M+1)))
          S = FLOP (P+DABS (E(M)))
          IF (P.EQ.S) GO TO 120
110     CONTINUE
120     CONTINUE
        P = D(L)
        IF (M.EQ.L) GO TO 240
        IF (J.EQ.30) GO TO 1000
        J = J+1
C
C ***      FORM SHIFT
        G = FLOP ((D(L+1)-P)/(2.0D0*E(L)))
        R = FLOP (DSQRT (G*G+1.0D0))
        G = FLOP (D(M)-P+E(L)/(G+DSIGN (R, G)))
        S = 1.0D0
        C = 1.0D0
        P = 0.0D0
        MML = M-L
C
C ***      FOR I = M-1 STEP -1 UNTIL L DO
        DO 200 II = 1, MML
          I = M-II
          F = FLOP (S*E(I))
          B = FLOP (C*E(I))
          IF (DABS (F).LT.DABS (G)) GO TO 150
          C = FLOP (G/F)
          R = FLOP (DSQRT (C*C+1.0D0))
          E(I+1) = FLOP (F*R)
          S = FLOP (1.0D0/R)
          C = FLOP (C*S)
          GO TO 160
C
150       CONTINUE
          S = FLOP (F/G)
          R = FLOP (DSQRT (S*S+1.0D0))
          E(I+1) = FLOP (G*R)
          C = FLOP (1.0D0/R)
          S = FLOP (S*C)
160       CONTINUE
          G = FLOP (D(I+1)-P)
          R = FLOP ((D(I)-G)*S+2.0D0*C*B)
          P = FLOP (S*R)
          D(I+1) = G+P
          G = FLOP (C*R-B)
          IF (JOB.EQ.0) GO TO 185
C
C ***      FORM VECTOR
          DO 180 K = 1, N
            F = Z(K,I+1)
            Z(K,I+1) = FLOP (S*Z(K,I)+C*F)
            Z(K,I) = FLOP (C*Z(K,I)-S*F)
180       CONTINUE
185       CONTINUE
200     CONTINUE
        D(L) = FLOP (D(L)-P)
        E(L) = G
        E(M) = 0.0D0
        GO TO 105
240   CONTINUE
C
C ***      ORDER EIGENVALUES AND EIGENVECTORS
      DO 300 II = 2, N
        I = II-1
        K = I
        P = D(I)
        DO 260 J = II, N
          IF (D(J).GE.P) GO TO 260
          K = J
          P = D(J)
260     CONTINUE
        IF (K.EQ.I) GO TO 300
        D(K) = D(I)
        D(I) = P
        IF (JOB.EQ.0) GO TO 285
        DO 280 J = 1, N
          P = Z(J,I)
          Z(J,I) = Z(J,K)
          Z(J,K) = P
280     CONTINUE
285     CONTINUE
300   CONTINUE
      GO TO 1001
C
C ***      SET ERROR -- NO CONVERGENCE TO AN EIGENVALUE AFTER 30 ITERATIONS
1000  CONTINUE
      IERR = L
C
1001  CONTINUE
      RETURN
      END
