C
C      ________________________________________________________
C     |                                                        |
C     |            COMPUTE THE DIAGONALIZATION OF A            |
C     |              SYMMETRIC TRIDIAGONAL MATRIX              |
C     |                                                        |
C     |    INPUT:                                              |
C     |                                                        |
C     |         V     --EITHER THE IDENTITY MATRIX OR THE SIM- |
C     |                 ILARITY TRANSFORMATION USED TO TRIDIAG-|
C     |                 ONALIZE THE COEFFICIENT MATRIX         |
C     |                                                        |
C     |         LV    --LEADING (ROW) DIMENSION OF ARRAY V     |
C     |                                                        |
C     |         D     --DIAGONAL                               |
C     |                                                        |
C     |         U     --SUPERDIAGONAL                          |
C     |                                                        |
C     |         N     --MATRIX DIMENSION                       |
C     |                                                        |
C     |    OUTPUT:                                             |
C     |                                                        |
C     |         E     --ARRAY OF EIGENVALUES                   |
C     |                                                        |
C     |         V     --ARRAY OF EIGENVECTORS                  |
C     |                                                        |
C     |    NOTE: THE ALGORITHM DESTROYS THE STARTING D AND U   |
C     |                                                        |
C     |    BUILTIN FUNCTIONS: ABS,AMAX1,SQRT                   |
C     |    PACKAGE SUBROUTINES: EIG3,SORT2                     |
C     |________________________________________________________|
C
      SUBROUTINE TDG(E,V,LV,D,U,N)
      INTEGER G,H,I,J,K,L,LL,LV,L1,M,N,NS
      REAL D(1),E(1),U(1),V(LV,1),A,B,F,P,Q,R,S,T,T0,T1,T2,W,X,Y,Z
      E(1) = D(1)
      IF ( N .LE. 1 ) RETURN
      J = 1
      H = 1
      D(1) = 1.
      T = 0.
      DO 10 I = 2,N
           E(I) = D(I)
           D(I) = 1.
           T = T + ABS(E(H)) + ABS(U(H))
           IF ( U(H) .EQ. 0. ) J = I
10         H = I
      B = 65536.**(-3)
      R = 1.
20    R = .5*R
      S = 1. + R
      IF ( S .GT. 1. ) GOTO 20
      T0 = R + R
      T2 = T0*T0
      NS = 50*N
      K = N
      LL = 0
      L = N - 1
      L1 = 0
      IF ( K .EQ. J ) GOTO 250
      GOTO 60
30    K = L
      L = L - 1
      J = 0
      DO 40 I = 1,L
40         IF ( U(I) .EQ. 0. ) J = I
      J = J + 1
      IF ( K .EQ. J ) GOTO 250
      T = 0.
      R = SQRT(D(J))
C     -----------------------------
C     |*** SET ERROR TOLERANCE ***|
C     -----------------------------
      DO 50 I = J,L
           S = SQRT(D(I+1))
           T = T + D(I)*ABS(E(I)) + R*ABS(U(I))*S
50         R = S
60    T = T + ABS(E(K))
      IF ( T .EQ. 0. ) T = 1.
      T1 = 1./(T*T0)
      GOTO 230
70    LL = LL + 1
      IF ( LL .GT. NS ) GOTO 330
      T = 1.
      DO 80 I = J,K
80         IF ( D(I) .LT. T ) T = D(I)
      IF ( T .GT. B ) GOTO 120
C     ----------------------------
C     |*** RESCALE THE MATRIX ***|
C     ----------------------------
      E(J) = E(J)*D(J)
      S = SQRT(D(J))
      D(J) = 1.
      DO 90 I = 1,N
90         V(I,J) = V(I,J)*S
      IF ( J .GT. 1 ) U(J-1) = U(J-1)*S
      H = J
      G = J + 1
      DO 110 M = G,K
           E(M) = E(M)*D(M)
           T = SQRT(D(M))
           D(M) = 1.
           U(H) = S*U(H)*T
           S = T
           DO 100 I = 1,N
100             V(I,M) = V(I,M)*T
110        H = M
120   T = D(K)
      S = D(L)
C     -----------------------
C     |*** COMPUTE SHIFT ***|
C     -----------------------
      CALL EIG3(F,F,S*E(L),T*E(K),S*U(L),T*U(L))
      H = J - 1
      I = J
      P = D(J)*E(J) - F
      Q = D(J)*U(J)
130   G = H
      H = I
      I = I + 1
      Y = D(H)
      Z = D(I)
C     -------------------------
C     |*** GIVENS ROTATION ***|
C     -------------------------
      IF ( ABS(P) .GT. ABS(Q) ) GOTO 140
      R = Y/Z
      S = P/Q
      T = R*S*S
      IF ( T .LT. 1. ) GOTO 190
      T = T/(1.+T)
      R = Z/Y
      S = Q/P
      GOTO 150
140   R = Z/Y
      S = Q/P
      T = R*S*S
      IF ( T .GT. 1. ) GOTO 180
      T = 1./(1.+T)
150   D(H) = Y*T
      D(I) = Z*T
      Y = S
      X = R*S
C     ------------------------------
C     |*** PROCESS EIGENVECTORS ***|
C     ------------------------------
      DO 160 M = 1,N
           T = V(M,H)
           S = V(M,I)
           V(M,H) = T + X*S
160        V(M,I) = S - Y*T
      IF ( H .EQ. J ) GOTO 170
      T = U(G) + X*A
      U(G) = T
      IF ( T .EQ. 0. ) J = H
170   T = E(H)
      P = U(H)
      Q = T + X*P
      R = P - Y*T
      W = E(I) - Y*P
      S = P + X*E(I)
      U(H) = S - Y*Q
      E(H) = Q + X*S
      E(I) = W - R*Y
      IF ( I .EQ. K ) GOTO 230
      Q = U(I)
      A = X*Q
      P = Z*W - F
      Q = Q*Z
      GOTO 130
180   T = T/(1.+T)
      R = Y/Z
      S = P/Q
      GOTO 200
190   T = 1./(1.+T)
200   D(H) = Z*T
      D(I) = Y*T
      Y = S
      X = R*S
C     ------------------------------
C     |*** PROCESS EIGENVECTORS ***|
C     ------------------------------
      DO 210 M = 1,N
           T = V(M,H)
           S = V(M,I)
           V(M,H) = X*T + S
210        V(M,I) = Y*S - T
      IF ( H .EQ. J ) GOTO 220
      T = X*U(G) + A
      U(G) = T
      IF ( T .EQ. 0. ) J = H
220   T = E(H)
      P = U(H)
      Q = X*T + P
      R = Y*P - T
      W = Y*E(I) - P
      S = X*P + E(I)
      U(H) = Y*S - Q
      E(H) = X*Q + S
      E(I) = Y*W - R
      IF ( I .EQ. K ) GOTO 230
      A = U(I)
      Q = A
      U(I) = Y*A
      P = Z*W - Y*F
      Q = Q*Z
      GOTO 130
230   T = D(K)
      S = D(L)
      Q = E(K)
      P = U(L)
C     ------------------------------
C     |*** TEST FOR CONVERGENCE ***|
C     ------------------------------
      IF ( ((T*P)*T1)*((S*P)*T1) .GT. 1. ) GOTO 70
      L1 = L1 + 1
      IF ( L1 .GT. 30 ) GOTO 240
      R = AMAX1(ABS(P),ABS(Q))
      IF ( R .EQ. 0. ) GOTO 240
      IF ( ABS((S*P)*(P/R)) .GT. T2*ABS((T*Q)*(Q/R)) ) GOTO 70
240   L1 = 0
      E(K) = Q*T
      K = L
      L = L - 1
      IF ( K .GT. J ) GOTO 230
250   E(J) = E(J)*D(J)
      IF ( J .GT. 2 ) GOTO 30
      IF ( J .EQ. 1 ) GOTO 270
      E(1) = E(1)*D(1)
C     ----------------------------------
C     |*** RESCALE THE EIGENVECTORS ***|
C     ----------------------------------
270   DO 280 J = 1,N
           T = SQRT(D(J))
           DO 280 I = 1,N
280             V(I,J) = V(I,J)*T
C     --------------------------------
C     |*** REORDER THE EIGENPAIRS ***|
C     --------------------------------
      CALL SORT2(E,D,U,N)
      DO 290 I = 1,N
           J = D(I)
290        U(J) = I
      DO 320 J = 1,N
           K = D(J)
           I = U(J)
           D(I) = K
           U(K) = I
           T = E(J)
           E(J) = E(K)
           E(K) = T
           S = 0.
           DO 300 I = 1,N
                T = V(I,K)
                V(I,K) = V(I,J)
                V(I,J) = T
300             S = S + T*T
           S = 1./SQRT(S)
           DO 310 I = 1,N
310             V(I,J) = S*V(I,J)
320   CONTINUE
      D(1) = N
      RETURN
330   K = N - K + 1
      WRITE(6,*) 'SINCE THE STOPPING CRITERION NOT SATISFIED'
      WRITE(6,*) 'AFTER',NS,'ITERATIONS, WE STOP WHILE COMPUTING'
      WRITE(6,*) 'EIGENVALUE NUMBER',K
      D(1) = K
      RETURN
      END