C
C      ________________________________________________________
C     |                                                        |
C     | COMPUTE THE DIAGONALIZATION OF A TRIDIAGONAL MATRIX FOR|
C     |  WHICH THE CROSS-DIAGONAL PRODUCTS ARE NONNEGATIVE AND |
C     |     A SUB I,I+1 = 0 IF AND ONLY IF A SUB I+1,I = 0     |
C     |                                                        |
C     |    INPUT:                                              |
C     |                                                        |
C     |         LV    --LEADING (ROW) DIMENSION OF ARRAY V     |
C     |                                                        |
C     |         L     --SUBDIAGONAL (CAN BE IDENTIFIED WITH U) |
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: SQRT                             |
C     |    PACKAGE SUBROUTINES: TDG                            |
C     |________________________________________________________|
C
      SUBROUTINE TDIAG(E,V,LV,L,D,U,N)
      INTEGER I,J,K,LV,N
      REAL D(1),E(1),L(1),U(1),V(LV,1),S,T
      DO 10 I = 1,N
10         V(I,1) = 0.
      V(1,1) = 1.
      IF ( N .GT. 1 ) GOTO 20
      E(1) = D(1)
      RETURN
C     ------------------------------------------
C     |*** SYMMETRIZE A NONSYMMETRIC MATRIX ***|
C     ------------------------------------------
20    K = 1
      S = L(1)
      DO 50 J = 2,N
           DO 30 I = 1,N
30              V(I,J) = 0.
           IF ( S .NE. 0. ) GOTO 40
           IF ( U(K) .NE. 0. ) GOTO 70
           V(J,J) = 1.
           S = L(J)
           GOTO 50
40         IF ( U(K) .EQ. 0. ) GOTO 70
           T = S/U(K)
           IF ( T .LT. 0. ) GOTO 60
           T = SQRT(T)
           V(J,J) = T
           U(K) = T*U(K)
           IF ( J .EQ. N ) GOTO 50
           U(J) = U(J)/T
           S = T*L(J)
50         K = J
      CALL TDG(E,V,LV,D,U,N)
      RETURN
60    WRITE(6,*) 'ERROR: SUBROUTINE TDIAG CAN ONLY BE USED WHEN'
      WRITE(6,*) 'THE CROSS-DIAGONAL PRODUCTS ARE NONNEGATIVE'
      STOP
70    WRITE(6,*) 'ERROR: SUBROUTINE TDIAG CAN ONLY BE USED IF THE'
      WRITE(6,*) 'CROSS-DIAGONAL PRODUCT A SUB I+1,I TIMES A SUB I,I+1'
      WRITE(6,*) 'VANISHES ONLY WHEN BOTH A SUB I+1,I AND A SUB I,I+1'
      WRITE(6,*) 'ARE ZERO'
      STOP
      END