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