C
C      ________________________________________________________
C     |                                                        |
C     |   FIND THE K-TH SMALLEST EIGENVALUE OF A TRIDIAGONAL   |
C     |  MATRIX WHOSE CROSS-DIAGONAL PRODUCTS ARE NONNEGATIVE  |
C     |USING BOTH THE BISECTION METHOD AND A NEWTON-LIKE METHOD|
C     |                                                        |
C     |    INPUT:                                              |
C     |                                                        |
C     |         K     --INDEX OF DESIRED EIGENVALUE (REPLACE K |
C     |                 BY -K TO GET K-TH LARGEST EIGENVALUE)  |
C     |                                                        |
C     |         L     --SUBDIAGONAL (CAN BE IDENTIFIED WITH U) |
C     |                                                        |
C     |         D     --DIAGONAL                               |
C     |                                                        |
C     |         U     --SUPERDIAGONAL                          |
C     |                                                        |
C     |         N     --MATRIX DIMENSION                       |
C     |                                                        |
C     |         W     --WORK ARRAY (LENGTH AT LEAST N)         |
C     |                                                        |
C     |    OUTPUT:                                             |
C     |                                                        |
C     |         E     --EIGENVALUE                             |
C     |                                                        |
C     |    BUILTIN FUNCTIONS: ABS,AMAX1                        |
C     |    PACKAGE SUBROUTINES: CP,EQL,INP,STM                 |
C     |________________________________________________________|
C
      SUBROUTINE TVAL(E,K,L,D,U,N,W)
      REAL D(1),L(1),U(1),W(1),E,S,T,Y,Z
      INTEGER I,K,M,N
      IF ( N .LE. 1 ) GOTO 40
C     -------------------------------------------------------
C     |*** BOUND EIGENVALUES USING GERSCHGORIN'S THEOREM ***|
C     -------------------------------------------------------
      M = N - 1
      Y = D(1)
      Z = D(1)
      S = 0.
      DO 10 I = 1,M
           W(I) = L(I)*U(I)
           IF ( W(I) .LT. 0. ) GOTO 50
           S = S + ABS(U(I))
           T = D(I) + S
           IF ( Z .LT. T ) Z = T
           T = D(I) - S
           IF ( Y .GT. T ) Y = T
10         S = ABS(L(I))
      T = D(N) + S
      IF ( Z .LT. T ) Z = T
      T = D(N) - S
      IF ( Y .GT. T ) Y = T
      M = K
      IF ( K .GT. 0 ) GOTO 20
      M = N + K + 1
20    T = 1.
30    T = .5*T
      S = 1. + T
      IF ( S .GT. 1. ) GOTO 30
      T = 8.*N*T*AMAX1(ABS(Y),ABS(Z))
      IF ( T .EQ. 0. ) GOTO 40
      CALL STM(E,Y,Z,M,T,D,W,N)
      RETURN
40    E = D(1)
      RETURN
50    WRITE(6,*) 'ERROR: SUBROUTINE TVAL CAN ONLY BE APPLIED'
      WRITE(6,*) 'WHEN THE CROSS-DIAGONAL PRODUCTS ARE NONNEGATIVE'
      STOP
      END
C%
C      ________________________________________________________
C     |                                                        |
C     |   FIND THE K-TH SMALLEST EIGENVALUE OF A TRIDIAGONAL   |
C     |  MATRIX WHOSE CROSS-DIAGONAL PRODUCTS ARE NONNEGATIVE  |
C     |USING BOTH THE BISECTION METHOD AND A NEWTON-LIKE METHOD|
C     |                                                        |
C     |    INPUT:                                              |
C     |                                                        |
C     |         Y,Z   --ENDS OF AN INTERVAL THAT CONTAINS THE  |
C     |                 DESIRED EIGENVALUE                     |
C     |                                                        |
C     |         K     --INDEX OF DESIRED EIGENVALUE            |
C     |                                                        |
C     |         T     --TOLERANCE (ITERATIONS CONTINUE UNTIL   |
C     |                 THE ERROR IN THE EIGENVALUE .LE. T)    |
C     |                                                        |
C     |         D     --DIAGONAL OF THE COEFFICIENT MATRIX A   |
C     |                                                        |
C     |         P     --CROSS-DIAGONAL PRODUCTS A SUB I+1,I    |
C     |                 TIMES A SUB I,I+1                      |
C     |                                                        |
C     |         N     --MATRIX DIMENSION                       |
C     |                                                        |
C     |    OUTPUT:                                             |
C     |                                                        |
C     |         E     --EIGENVALUE                             |
C     |                                                        |
C     |    BUILTIN FUNCTIONS: ABS,AMAX1,SIGN,SQRT              |
C     |    PACKAGE SUBROUTINES: CP,EQL,INP                     |
C     |________________________________________________________|
C
      SUBROUTINE STM(E,Y,Z,K,T,D,P,N)
      REAL A,B,E,FL,FR,L,Q,R,S,T,U,V,W,Y,Z
      INTEGER I,J,K,M,N
      REAL D2,D3,D4,P2,P3,P4,D34,D42,D23
      REAL D(1),F(4),H(4),P(1),X(4)
      DOUBLE PRECISION G(4),C,GL,GR
      L = Y
      R = Z
      IF ( L .LT. R ) GOTO 10
      L = Z
      R = Y
10    IF ( N .EQ. 1 ) GOTO 210
      S = T + T
      U = 1.
20    U = .5*U
      A = 1. + U
      IF ( A .GT. 1. ) GOTO 20
      U = 5.*U
      V = .5*U
      CALL CP(L,FL,GL,I,D,P,N)
      CALL CP(R,FR,GR,J,D,P,N)
      IF ( I .GE. K ) GOTO 190
      IF ( J .LT. K ) GOTO 200
C     --------------------------------
C     |*** ISOLATE THE EIGENVALUE ***|
C     --------------------------------
30    E = R - L
      IF ( E .LE. S ) GOTO 180
      IF ( E .LE. U*(ABS(L)+ABS(R)) ) GOTO 180
      IF ( J .EQ. I+1 ) GOTO 70
40    A = .5*(L+R)
      CALL CP(A,B,C,M,D,P,N)
      IF ( K .LE. M ) GOTO 50
      L = A
      I = M
      FL = B
      GL = C
      GOTO 30
50    R = A
      J = M
      FR = B
      GR = C
      GOTO 30
60    E = R - L
      IF ( E .LE. S ) GOTO 180
      IF ( E .LE. U*(ABS(L)+ABS(R)) ) GOTO 180
70    X(1) = L
      F(1) = FL
      G(1) = GL
      X(2) = R
      F(2) = FR
      G(2) = GR
      CALL EQL(X,F,G,H,2)
      IF ( H(1) .EQ. H(2) ) GOTO 160
C     ---------------------
C     |*** SECANT STEP ***|
C     ---------------------
      A = X(1) - H(1)*(X(1)-X(2))/(H(1)-H(2))
      Q = A
      W = AMAX1(T,V*(ABS(L)+ABS(R)))
      IF ( ABS(A-L) .LT. W ) A = L + W
      IF ( ABS(A-R) .LT. W ) A = R - W
      CALL CP(A,B,C,J,D,P,N)
      IF ( I .GE. J ) GOTO 80
      R = A
      FR = B
      GR = C
      GOTO 90
80    L = A
      FL = B
      GL = C
90    X(3) = A
      F(3) = B
      G(3) = C
      W = R - L
      IF ( W .LE. S ) GOTO 220
      IF ( W .LE. U*(ABS(L)+ABS(R)) ) GOTO 220
      CALL EQL(X,F,G,H,3)
C     --------------------------------------
C     |*** QUADRATIC INTERPOLATION STEP ***|
C     --------------------------------------
      CALL INP(A,X(1),X(2),X(3),H(1),H(2),H(3),L,R)
      B = L
      IF ( ABS(A-L) .GT. ABS(A-R) ) B = R
C     ------------------------------------
C     |*** APPLY PSEUDO-NEWTON METHOD ***|
C     ------------------------------------
100   Q = A
      W = AMAX1(T,V*(ABS(L)+ABS(R)))
      IF ( ABS(A-L) .LT. W ) GOTO 110
      IF ( ABS(A-R) .GT. W ) GOTO 130
110   IF ( A+A .GT. L+R ) GOTO 120
      A = L + W
      GOTO 130
120   A = R - W
130   IF ( A .LE. L ) GOTO 160
      IF ( A .GE. R ) GOTO 160
      E = .5*E
      IF ( E .LT. ABS(B-A) ) GOTO 160
      CALL CP(A,B,C,J,D,P,N)
      IF ( I .GE. J ) GOTO 140
      R = A
      FR = B
      GR = C
      GOTO 150
140   L = A
      FL = B
      GL = C
150   W = R - L
      IF ( W .LE. S ) GOTO 220
      IF ( W .LE. U*(ABS(L)+ABS(R)) ) GOTO 220
      X(4) = A
      F(4) = B
      G(4) = C
      CALL EQL(X,F,G,H,4)
      IF ( X(1) .LT. L ) GOTO 160
      IF ( X(1) .GT. R ) GOTO 160
      B = X(1)
      D4 = X(4) - B
      D3 = X(3) - B
      D2 = X(2) - B
      D34 = X(3) - X(4)
      D42 = X(4) - X(2)
      D23 = X(2) - X(3)
      P2 = 1./D34
      P3 = 1./D42
      P4 = 1./D23
      P2 = (H(2)-H(1))/(D2*(1.+(D2/D3)*P2*D42+(D2/D4)*P2*D23))
      P3 = (H(3)-H(1))/(D3*(1.+(D3/D2)*P3*D34+(D3/D4)*P3*D23))
      P4 = (H(4)-H(1))/(D4*(1.+(D4/D2)*P4*D34+(D4/D3)*P4*D42))
      P2 = P2 + P3 + P4
      IF ( P2 .EQ. 0. ) GOTO 160
      A = B - H(1)/P2
      GOTO 100
C     --------------------------
C     |*** BISECTION METHOD ***|
C     --------------------------
160   A = .5*(L+R)
      CALL CP(A,B,C,J,D,P,N)
      IF ( I .GE. J ) GOTO 170
      R = A
      FR = B
      GR = C
      GOTO 60
170   L = A
      FL = B
      GL = C
      GOTO 60
180   E = .5*(L+R)
      RETURN
190   E = L
      RETURN
200   E = R
      RETURN
210   E = D(1)
      IF ( L .GT. E ) E = L
      IF ( R .LT. E ) E = R
      RETURN
220   E = Q
      RETURN
      END
C%
      SUBROUTINE INP(A,X,Y,Z,U,V,W,L,R)
      REAL A,B,C,L,P,Q,R,S,T,U,V,W,X,Y,Z
      S = Z - X
      T = (Y-X)/S
      A = (U-V)/T + (W-V)/(1.-T)
      IF ( A .EQ. 0. ) GOTO 40
      B = .5*(W-U)/A - .5
      C = U/A
      T = SQRT(ABS(C))
      IF ( ABS(B) .LT. SIGN(T,C) ) GOTO 60
      T = AMAX1(T,ABS(B))
      IF ( T .EQ. 0. ) GOTO 50
      Q = 1./T
      P = SQRT((Q*B)**2 - Q*C*Q)
      P = T*P
      IF ( ABS(P+B) .GT. ABS(P-B) ) GOTO 10
      Q = P - B
      GOTO 20
10    Q = -(B+P)
20    P = C/Q
      Q = X + S*Q
      P = X + S*P
      IF ( Q .LT. L ) GOTO 30
      IF ( Q .GT. R ) GOTO 30
      A = Q
      RETURN
30    A = P
      RETURN
40    IF ( U .EQ. W ) GOTO 50
      A = X + S*U/(U-W)
      RETURN
50    A = L
      RETURN
60    A = X - S*B
      RETURN
      END
C%
C     ------------------------------------------------
C     |*** EVALUATE CHARACTERISTIC POLYNOMIAL AND ***|
C     |              SIGN ALTERNATIONS               |
C     ------------------------------------------------
      SUBROUTINE CP(T,B,V,L,D,P,N)
      REAL D(1),P(1)
      REAL A,B,C,E,F,T,U,Z
      DOUBLE PRECISION V
      INTEGER I,J,L,M,N
      L = 0
      Z = 0.
      V = Z
      IF ( N .GT. 1 ) GOTO 10
      B = D(1) - T
      IF ( B .LE. Z ) L = 1
      RETURN
10    F = 65536.**4
      E = 1./F
      U = 16.
      M = 0
      I = 1
20    C = 1.
      B = D(I) - T
30    IF ( B .LE. Z ) GOTO 70
      IF ( I .GE. N ) GOTO 130
40    J = I
      I = I + 1
      A = (D(I)-T)*B - P(J)*C
      C = B
      B = A
50    A = ABS(B)
      IF ( A .GT. F ) GOTO 60
      IF ( A .GT. E ) GOTO 30
      IF ( A .EQ. Z ) GOTO 70
      C = C*F
      B = B*F
      V = V - U
      GOTO 50
60    C = C*E
      B = B*E
      V = V + U
      GOTO 50
70    L = L + 1
      IF ( I .GE. N ) GOTO 130
      IF ( B .LT. Z ) GOTO 90
      IF ( P(I) .GT. 0. ) GOTO 90
      I = I + 1
      M = 1
      V = Z
      GOTO 20
80    IF ( B .GE. Z ) GOTO 120
      IF ( I .GE. N ) GOTO 130
90    J = I
      I = I + 1
      A = (D(I)-T)*B - P(J)*C
      C = B
      B = A
100   A = ABS(B)
      IF ( A .GT. F ) GOTO 110
      IF ( A .GT. E ) GOTO 80
      IF ( A .EQ. Z ) GOTO 120
      C = C*F
      B = B*F
      V = V - U
      GOTO 100
110   C = C*E
      B = B*E
      V = V + U
      GOTO 100
120   L = L + 1
      IF ( I .GE. N ) GOTO 130
      IF ( B .GT. Z ) GOTO 40
      IF ( P(I) .GT. 0. ) GOTO 40
      I = I + 1
      M = 1
      V = Z
      GOTO 20
130   IF ( M .EQ. 1 ) GOTO 160
      IF ( B .EQ. 0. ) GOTO 160
      A = 1./U
      IF ( ABS(B) .LT. A ) GOTO 150
140   IF ( ABS(B) .LT. 1. ) RETURN
      B = B*A
      V = V + 1.D0
      GOTO 140
150   B = B*U
      V = V - 1.D0
      IF ( ABS(B) .LT. A ) GOTO 150
      RETURN
160   V = Z
      B = Z
      RETURN
      END
C%
      SUBROUTINE EQL(X,F,G,H,N)
      REAL X(1),F(1),H(1),R,S,Z
      DOUBLE PRECISION G(1),U
      INTEGER I,J,K,L,M,N,NM1
      Z = 0.
      I = 1
      R = X(N)
      S = F(N)
      U = G(N)
      IF ( S .EQ. Z ) GOTO 30
10    IF ( F(I) .EQ. Z ) GOTO 20
      IF ( U .LT. G(I) ) GOTO 30
      IF ( U .GT. G(I) ) GOTO 20
      IF ( ABS(F(I)) .GE. ABS(S) ) GOTO 30
20    I = I + 1
      IF ( I .LT. N ) GOTO 10
      GOTO 50
30    M = N + I
      NM1 = N - 1
      DO 40 J = I,NM1
           K = M - J
           L = K - 1
           X(K) = X(L)
           F(K) = F(L)
40         G(K) = G(L)
      X(I) = R
      F(I) = S
      G(I) = U
50    U = G(N)
      DO 80 I = 1,N
           IF ( F(I) .EQ. Z ) GOTO 60
           IF ( G(I)-U .GT. -99.D0 ) GOTO 70
60         H(I) = Z
           GOTO 80
70         H(I) = F(I)*16.D0**(G(I)-U)
80    CONTINUE
      RETURN
      END