C C ________________________________________________________ C | | C | FIND ALL EIGENVALUES OF A TRIDIAGONAL MATRIX WITH NON- | C | NEGATIVE CROSS-DIAGONAL PRODUCTS | C | | C | INPUT: | 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 --EIGENVALUES | C | | C | L,D,U --ORIGINAL L, D, AND U ARE UNTOUCHED | C | | C | BUILTIN FUNCTIONS: ABS,SQRT | C | PACKAGE SUBROUTINES: NEWTON,PWK,SORT | C |________________________________________________________| C SUBROUTINE TVALS(E,L,D,U,N,W) REAL E(1),L(1),D(1),U(1),W(1),Q,R,S,T,V INTEGER I,J,K,M,N IF ( N .GT. 1 ) GOTO 10 E(1) = D(1) RETURN 10 M = N - 1 DO 20 I = 1,M 20 W(I) = L(I)*U(I) CALL PWK(E,D,W,N,W) IF ( W(1) .EQ. 1. ) GOTO 70 IF ( W(1) .EQ. 2. ) GOTO 80 T = 1. 30 T = T/2. S = 1. + T IF ( S .GT. 1. ) GOTO 30 T = SQRT(T) DO 40 I = 1,M 40 W(I) = L(I)*U(I) C -------------------------------------------------- C |*** REFINE EIGENVALUES USING NEWTON'S METHOD ***| C -------------------------------------------------- DO 60 I = 1,N Q = 0. J = 0 V = E(I) 50 CALL NEWTON(V,S,K,D,W,N) J = J + 1 IF ( J .GT. 2 ) GOTO 60 S = ABS(S) R = ABS(V) + S IF ( Q .LT. R ) Q = R IF ( S .GT. T*Q ) GOTO 50 60 E(I) = V CALL SORT(E,W,N) RETURN 70 WRITE(6,*) 'THE QR METHOD DID NOT CONVERGE' RETURN 80 WRITE(6,*) 'ERROR: TO USE SUBROUTINE PVALS, THE CROSS-DIAGONAL' WRITE(6,*) 'PRODUCT L(I)*U(I) MUST BE NONNEGATIVE FOR EVERY I' RETURN END C% C C ________________________________________________________ C | | C | FIND ALL EIGENVALUES OF A TRIDIAGONAL MATRIX WITH NON- | C | NEGATIVE CROSS DIAGONAL PRODUCTS USING THE PWK VERSION | C | OF THE QR ALGORITHM | C | | C | INPUT: | C | | C | E --ARRAY FOR EIGENVALUES WHICH CAN BE | C | IDENTIFIED WITH D | C | | C | D --DIAGONAL OF THE 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 | W --WORK ARRAY WHICH CAN BE IDENTIFIED WITH| C | U (LENGTH AT LEAST N) | C | | C | OUTPUT: | C | | C | E --EIGENVALUES | C | | C | D,U --THE ORIGINAL D AND U ARE UNTOUCHED | C | UNLESS W IS IDENTIFIED WITH U OR E IS | C | IDENTIFIED WITH D | C | | C | BUILTIN FUNCTIONS: ABS,SQRT | C |________________________________________________________| C SUBROUTINE PWK(E,D,U,N,W) REAL E(1),D(1),U(1),W(1),A,B,C,G,H,OC,OG,P,R,S,T,V INTEGER I,J,K,KK,L,M,N E(N) = D(N) IF ( N .EQ. 1 ) RETURN K = 1 M = N - 1 L = N + 1 T = 1. DO 20 I = 1,M E(I) = D(I) J = L - I W(J) = U(J-1) IF ( W(J) .GT. 0. ) GOTO 20 IF ( W(J) .EQ. 0. ) GOTO 10 W(1) = 2. RETURN 10 IF ( T .EQ. 1. ) K = J T = 0. 20 CONTINUE W(1) = 0. T = 1. 30 T = T/2. S = 1. + T IF ( S .GT. 1. ) GOTO 30 L = N 40 M = 0 KK = 0 V = 0. 50 P = W(L) S = SQRT(P) H = E(L) R = S + ABS(H) IF ( V .LT. R ) V = R IF ( S .GT. T*V ) GOTO 100 IF ( S .LE. T*R ) GOTO 60 V = R IF ( KK .EQ. 0 ) GOTO 90 60 L = L - 1 IF ( L .EQ. 1 ) RETURN IF ( K .LE. L ) GOTO 40 K = 1 DO 70 I = 2 , L 70 IF ( W(I) .EQ. 0. ) K = I GOTO 40 80 W(L) = S*P E(L) = G + H GOTO 50 90 KK = 1 100 M = M + 1 IF ( M .GT. 40 ) GOTO 130 R = .5*(E(L-1)-H) IF ( R .LT. 0. ) H = H + P/(ABS(R)+SQRT(R*R+P)) IF ( R .GE. 0. ) H = H - P/(ABS(R)+SQRT(R*R+P)) C = 1. S = 0. I = K G = E(I) - H P = G*G C ------------------------------- C |*** THE QR TRANSFORMATION ***| C ------------------------------- 110 J = I I = I + 1 IF ( I .GT. L ) GOTO 80 B = W(I) R = P + B A = S*R IF ( A .EQ. 0. ) K = J W(J) = A OC = C C = P/R S = B/R A = E(I) OG = G G = C*(A-H) - S*OG E(J) = OG + (A-G) IF ( C .EQ. 0. ) GOTO 120 P = G*G/C GOTO 110 120 P = OC*B GOTO 110 130 W(1) = 1. RETURN END