C
C      ________________________________________________________
C     |                                                        |
C     | APPLY 1 STEP OF NEWTON'S METHOD TO THE CHARACTERISTIC  |
C     |           POLYNOMIAL OF A TRIDIAGONAL MATRIX           |
C     |                                                        |
C     |    INPUT:                                              |
C     |                                                        |
C     |         T     --ESTIMATE FOR EIGENVALUE                |
C     |                                                        |
C     |         D     --DIAGONAL OF COEFFICIENT MATRIX A       |
C     |                                                        |
C     |         U     --THE CROSS DIAGONAL PRODUCTS A SUB I,I+1|
C     |                 TIMES A SUB I+1,I                      |
C     |                                                        |
C     |         N     --MATRIX DIMENSION                       |
C     |                                                        |
C     |    OUTPUT:                                             |
C     |                                                        |
C     |         T     --IMPROVED ESTIMATE FOR EIGENVALUE       |
C     |                                                        |
C     |         S     --NEWTON STEP                            |
C     |                                                        |
C     |         K     --THE NUMBER OF EIGENVALUES LESS THAN OR |
C     |                 EQUAL TO THE ORIGINAL T (ONLY VALID    |
C     |                 WHEN EVERY ELEMENT OF U IS NONNEGATIVE)|
C     |                                                        |
C     |    BUILTIN FUNCTIONS: ABS,SIGN                         |
C     |    PACKAGE FUNCTIONS: ISIG                             |
C     |________________________________________________________|
C
      SUBROUTINE NEWTON(T,S,K,D,U,N)
      REAL A,B,C,E,F,G,H,P,Q,S,T,X,Y,Z
      REAL D(1),U(1)
      INTEGER I,ISIG,J,K,L,N
      K = 0
      IF ( N .GT. 1 ) GOTO 10
      S = D(1) - T
      T = D(1)
      IF ( S .LE. 0. ) K = 1
      RETURN
10    Z = 0.
      Y = 65536.**4
      X = 1./Y
      C = 2.
      B = C*(D(1)-T)
      E = -C
      F = Z
      G = Z
      H = Z
      I = 1
      J = 1
      L = 1
      IF ( ISIG(B) .EQ. L ) GOTO 20
      K = K + 1
      L = -1
      IF ( B .LT. Z ) GOTO 20
      IF ( U(I) .EQ. Z ) GOTO 90
20    P = U(I)
      I = I + 1
      Q = D(I) - T
      A = Q*B - P*C
      IF ( ISIG(A) .EQ. L ) GOTO 30
      K = K + 1
      L = -L
      IF ( A .NE. 0. ) GOTO 30
      IF ( I .GE. N ) GOTO 140
      IF ( U(I) .EQ. 0. ) GOTO 90
30    C = B
      B = A
      A = Q*E - P*F - C
      F = E
      E = A
      A = Q*G - P*H - F
      H = G
      G = A
      IF ( I .GE. N ) GOTO 60
      A = ABS(B)
      IF ( ABS(E) .GT. A ) A = ABS(E)
      IF ( ABS(G) .GT. A ) A = ABS(G)
40    IF ( A .GT. Y ) GOTO 50
      IF ( A .GT. X ) GOTO 20
      IF ( A .EQ. Z ) GOTO 20
      A = A*Y
      B = B*Y
      C = C*Y
      E = E*Y
      F = F*Y
      G = G*Y
      H = H*Y
      GOTO 40
50    A = A*X
      B = B*X
      C = C*X
      E = E*X
      F = F*X
      G = G*X
      H = H*X
      GOTO 40
60    IF ( E .EQ. Z ) GOTO 150
      X = B/E
      IF ( G .EQ. Z ) GOTO 70
      Y = E/(G+G)
      S = X*Y/(X-Y)
      GOTO 80
70    S = -X
80    T = T + S
      RETURN
90    I = I + 1
      B = D(I) - T
      C = 1.
      L = 1
      IF ( ISIG(B) .EQ. L ) GOTO 100
      K = K + 1
      L = -1
      IF ( B .LT. 0. ) GOTO 100
      IF ( I .GE. N ) GOTO 140
      IF ( U(I) .EQ. 0. ) GOTO 90
100   IF ( I .GE. N ) GOTO 140
      J = I
      I = I + 1
      A = (D(I)-T)*B - U(J)*C
      IF ( ISIG(A) .EQ. L ) GOTO 110
      K = K + 1
      L = -L
      IF ( A .NE. 0. ) GOTO 110
      IF ( I .GE. N ) GOTO 140
      IF ( U(I) .EQ. 0. ) GOTO 90
110   C = B
      B = A
120   A = ABS(B)
      IF ( A .GT. Y ) GOTO 130
      IF ( A .GT. X ) GOTO 100
      IF ( A .EQ. Z ) GOTO 100
      C = C*Y
      B = B*Y
      GOTO 120
130   C = C*X
      B = B*X
      GOTO 120
140   S = 0.
      RETURN
150   S = T + SIGN(1.,T)
      T = T + S
      RETURN
      END
      INTEGER FUNCTION ISIG(X)
      REAL X
      IF ( X .LT. 0. ) GOTO 20
      IF ( X .EQ. 0. ) GOTO 10
      ISIG = 1
      RETURN
10    ISIG = 0
      RETURN
20    ISIG = -1
      RETURN
      END