C C ________________________________________________________ C | | C |COMPUTE ONE OR TWO EIGENVALUES OF LARGEST MAGNITUDE AND | C | THE CORRESPONDING EIGENVECTORS USING THE POWER METHOD | C | | C | INPUT: | C | | C | X --COMPLEX ARRAY CONTAINING REAL STARTING | C | GUESS FOR THE EIGENVECTOR SPACE | C | | C | N --MATRIX DIMENSION | C | | C | NDIGIT--DESIRED NUMBER OF CORRECT DIGITS | C | | C | LIMIT --MAXIMUM NUMBER OF ITERATIONS | C | | C | MULT --NAME OF SUBROUTINE TO MULTIPLY MATRIX A| C | BY VECTOR (NAME EXTERNAL IN MAIN PROG.)| C | MULT(P,V) STORES IN P THE PRODUCT AV | C | WHERE P AND V ARE REAL ARRAYS | C | | C | OUTPUT: | C | | C | X,Y --COMPLEX ARRAYS CONTAINING EIGENVECTORS | C | | C | EX,EY --COMPLEX VARIABLES CONTAINING THE EIGEN-| C | VALUES (IF EY IS ZERO, THEN IGNORE THE | C | (EY,Y) EIGENPAIR) | C | | C | DIF --INPUT FOR SUBROUTINE WHATIS | C | | C | SIZE --INPUT FOR SUBROUTINE WHATIS | C | | C | BUILTIN FUNCTIONS: ABS,CABS,CMPLX,SQRT | C |________________________________________________________| C SUBROUTINE POWER(EX,X,EY,Y,N,DIF,SIZE,NDIGIT,LIMIT,MULT) REAL X(1),Y(1),A,B,C,D,DIF,E,F,G,Q,R,S,SIZE,T COMPLEX EX,EY INTEGER I,J,K,L,LIMIT,M,N,NDIGIT,N2 M = N + 1 N2 = N + N Q = 10.**(-2*NDIGIT) L = 0 T = 0. DO 10 I = 1,N2,2 10 T = T + X(I)**2 IF ( T .NE. 0. ) GOTO 20 WRITE(6,*) 'STARTING GUESS FOR SUBROUTINE POWER CANNOT BE ZERO' STOP C ----------------------------------- C |*** NORMALIZE STARTING VECTOR ***| C ----------------------------------- 20 T = 1./SQRT(T) J = 1 DO 30 I = 1,N X(I) = T*X(J) 30 J = J + 2 C ------------------------- C |*** START ITERATION ***| C ------------------------- 40 DO 50 I = 1,N 50 Y(I) = X(I) CALL MULT(Y,X) DO 60 I = 1,N 60 Y(I+N) = Y(I) CALL MULT(Y(M),Y) C -------------------------------------------- C |*** COMPUTE DISTANCE FROM X SUB K+1 TO ***| C | SPACE SPANNED BY X SUB K AND X SUB K-1 | C -------------------------------------------- A = 0. DO 70 I = 1,N 70 A = A + Y(I+N)**2 IF ( A .EQ. 0. ) GOTO 440 A = 1./SQRT(A) S = 0. DO 80 I = 1,N S = S + X(I)*Y(I) J = I + N 80 Y(J) = Y(J)*A A = X(1) B = 1. + ABS(A) IF ( A .LT. 0. ) GOTO 90 A = A + 1. S = S + Y(1) GOTO 100 90 A = A - 1. S = S - Y(1) 100 S = S/B Y(1) = Y(1) - A*S C = 0. D = 0. DO 110 I = 2,N T = Y(I) - S*X(I) Y(I) = T C = C + T*X(I) 110 D = D + T*T IF ( D .NE. 0. ) GOTO 120 C = X(2) D = 1. Y(2) = 1. 120 D = 1./SQRT(D) C = C/B Y(1) = -C*A*D DO 130 I = 2,N 130 Y(I) = D*(Y(I)-C*X(I)) S = 0. T = 0. DO 140 I = 1,N J = I + N S = S + Y(J)*X(I) 140 T = T + Y(J)*Y(I) D = 0. DO 150 I = 1,N A = Y(I+N) D = D + (A-S*X(I)-T*Y(I))**2 150 X(I) = A L = L + 2 IF ( L .GE. LIMIT ) GOTO 160 IF ( D .GT. Q ) GOTO 40 C --------------------------------------------------------- C |*** IMPROVE ACCURACY USING ORTHOGONAL STARTING PAIR ***| C --------------------------------------------------------- 160 S = 0. DO 170 I = 1,N 170 S = S + X(I)*Y(I) A = X(1) B = 1. + ABS(A) IF ( A .LT. 0. ) GOTO 180 A = A + 1. S = S + Y(1) GOTO 190 180 A = A - 1. S = S - Y(1) 190 S = S/B Y(1) = Y(1) - A*S C = 0. D = 0. DO 200 I = 2,N T = Y(I) - S*X(I) Y(I) = T C = C + T*X(I) 200 D = D + T*T IF ( D .NE. 0. ) GOTO 210 C = X(2) D = 1. Y(2) = 1. 210 C = C/B D = 1./SQRT(D) Y(1) = -C*A*D Y(M) = Y(1) X(M) = X(1) DO 220 I = 2,N J = I + N T = D*(Y(I)-C*X(I)) Y(I) = T Y(J) = T 220 X(J) = X(I) CALL MULT(X,X(M)) CALL MULT(Y,Y(M)) A = 0. B = 0. C = 0. D = 0. S = 0. T = 0. DO 230 I = 1,N J = I + N F = X(I) G = Y(I) S = S + F*F T = T + G*G A = A + F*X(J) B = B + G*X(J) C = C + F*Y(J) 230 D = D + G*Y(J) IF ( T .EQ. 0. ) GOTO 490 IF ( S .EQ. 0. ) GOTO 540 F = 1./SQRT(S) G = 1./SQRT(T) R = 0. S = 0. T = 0. DO 240 I = 1,N J = I + N R = R + (X(I)-A*X(J))**2 S = S + (X(I)-A*X(J)-C*Y(J))**2 T = T + (Y(I)-B*X(J)-D*Y(J))**2 X(I) = F*X(I) 240 Y(I) = G*Y(I) T = F*F*S + G*G*T S = F*F*R L = L + 2 IF ( L .GE. LIMIT ) GOTO 280 IF ( T .LE. Q ) GOTO 290 IF ( S .GT. Q ) GOTO 160 C -------------------------------- C |*** SINGLE REAL EIGENVALUE ***| C -------------------------------- 250 SIZE = 3*L + 1 260 DIF = SQRT(S) EX = CMPLX(A,0.) EY = (0.,0.) J = N2 I = N 270 K = J - 1 Y(K) = 0. Y(J) = 0. X(K) = X(I) X(J) = 0. J = J - 2 I = I - 1 IF ( J .GT. 0 ) GOTO 270 RETURN C ----------------------------- C |*** ITERATIONS AT LIMIT ***| C ----------------------------- 280 SIZE = 3*L + 2 IF ( S+S .LT. T ) GOTO 260 GOTO 300 290 SIZE = 3*L + 1 300 DIF = SQRT(T) S = .5*(D-A) T = B*C R = S*S + T IF ( R .LT. 0. ) GOTO 400 C ------------------------------ C |*** TWO REAL EIGENVALUES ***| C ------------------------------ R = ABS(S) + SQRT(R) IF ( S .LT. 0. ) GOTO 310 EX = CMPLX(A-T/R,0.) EY = CMPLX(A+R,0.) GOTO 320 310 EX = CMPLX(A+T/R,0.) EY = CMPLX(A-R,0.) 320 S = CABS(A-EX) + ABS(B) T = CABS(D-EX) + ABS(C) IF ( S .GT. T ) GOTO 340 IF ( T .NE. 0. ) GOTO 330 E = 0. F = 1. A = 1. B = 0. GOTO 380 330 E = C/T F = (EX-D)/T GOTO 350 340 E = (EX-A)/S F = B/S 350 S = CABS(A-EY) + ABS(B) T = CABS(D-EY) + ABS(C) IF ( S .GT. T ) GOTO 370 IF ( T .NE. 0. ) GOTO 360 E = 0. F = 1. A = 1. B = 0. GOTO 380 360 A = C/T B = (EY-D)/T GOTO 380 370 A = (EY-A)/S B = B/S 380 J = 0 DO 390 I = M,N2 J = J + 2 K = J - 1 S = F*X(I) + E*Y(I) T = B*X(I) + A*Y(I) Y(J) = 0. X(J) = 0. X(K) = S 390 Y(K) = T RETURN C ----------------------------- C |*** COMPLEX EIGENVALUES ***| C ----------------------------- 400 Q = A + S R = SQRT(-R) EX = CMPLX(Q,R) EY = CMPLX(Q,-R) S = CABS(A-EX) + ABS(B) T = CABS(D-EX) + ABS(C) J = 0 IF ( S .GT. T ) GOTO 420 B = C/T C = (Q-D)/T A = R/T DO 410 I = M,N2 J = J + 2 K = J - 1 T = C*X(I) + B*Y(I) S = A*X(I) X(K) = T Y(K) = T X(J) = S 410 Y(J) = -S RETURN 420 C = B/S B = (Q-A)/S D = R/S DO 430 I = M,N2 J = J + 2 K = J - 1 T = C*X(I) + B*Y(I) S = D*Y(I) X(K) = T Y(K) = T X(J) = S 430 Y(J) = -S RETURN 440 EX = (0.,0.) EY = (0.,0.) DIF = 0. SIZE = 3*L + 1 DO 450 I = 1,N 450 IF ( Y(I) .NE. 0. ) GOTO 470 J = N2 I = N 460 K = J - 1 X(K) = X(I) X(J) = 0. Y(K) = 0. Y(J) = 0. J = J - 2 I = I - 1 IF ( J .GT. 0 ) GOTO 460 RETURN 470 J = N2 I = N 480 K = J - 1 X(K) = Y(I) X(J) = 0. Y(K) = 0. Y(J) = 0. J = J - 2 I = I - 1 IF ( J .GT. 0 ) GOTO 480 RETURN 490 IF ( S .NE. 0. ) GOTO 510 DO 500 I = 1,N 500 X(I) = X(I+N) GOTO 440 510 S = 1./SQRT(S) A = A*S DO 520 I = 1,N 520 X(I) = X(I)*S S = 0. DO 530 I = 1,N 530 S = S + (X(I)-A*X(I+N))**2 IF ( S .GT. Q ) GOTO 40 GOTO 250 540 T = 1./SQRT(T) DO 550 I = 1,N 550 X(I) = T*Y(I) GOTO 40 END