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