C C ________________________________________________________ C | | C | COMPUTE SEVERAL EIGENPAIRS 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 | NE --TOTAL NUMBER OF EIGENPAIRS TO COMPUTE | 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 | W --WORK ARRAY (LENGTH AT LEAST NE(N+2)-2) | 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 | NE --NUMBER OF EIGENPAIRS LEFT TO COMPUTE | C | | C | DIF --INPUT FOR SUBROUTINE WHATIS | C | | C | SIZE --INPUT FOR SUBROUTINE WHATIS | C | | C | BUILTIN FUNCTIONS: ABS,AIMAG,CABS,CMPLX,REAL,SQRT | C | PACKAGE SUBROUTINES: CVC,HSE,PX,RVC,STP,XP | C |________________________________________________________| C SUBROUTINE MPOWER(EX,X,EY,Y,N,NE,DIF,SIZE,NDIGIT,LIMIT,MULT,W) REAL W(1),X(1),Y(1),A,B,C,D,DIF,E,F,G,Q,R,S,SIZE,T COMPLEX EX,EY INTEGER H,I,J,K,L,LIMIT,M,M1,M2,N,NDIGIT,NE,NP,NS INTEGER N2,N3,N4,N5,O,P,P1 EXTERNAL MULT DATA NP/0/ IF ( NE .NE. 0 ) GOTO 10 WRITE(6,*) 'SINCE THE REQUESTED EIGENPAIRS HAVE BEEN COMPUTED,' WRITE(6,*) 'SUBROUTINE MPOWER WILL QUIT' RETURN 10 IF ( NE .EQ. NP ) GOTO 20 H = 0 NS = NE N3 = (NS*(3+N+N-NS))/2 20 O = H + 1 P = O + 1 P1 = P + 1 M = N + 1 M1 = N + H M2 = N + O N2 = N + N N4 = (O*P)/2 + N3 - 2 N5 = (H*(N+N-O))/2 Q = 10.**(-2*NDIGIT) L = 0 J = 1 DO 30 I = 1,N X(I) = X(J) 30 J = J + 2 IF ( H .EQ. 0 ) GOTO 50 CALL XP(X,W(M),H,N) DO 40 I = 1,H 40 X(I) = 0. 50 T = 0. DO 60 I = O,N 60 T = T + X(I)**2 IF ( T .NE. 0. ) GOTO 70 WRITE(6,*) 'STARTING GUESS FOR SUBROUTINE MPOWER CANNOT BE ZERO' WRITE(6,*) 'AND IT CANNOT LIE IN THE SPACE SPANNED BY THE' WRITE(6,*) 'PREVIOUSLY COMPUTED EIGENVECTORS' STOP C ----------------------------------- C |*** NORMALIZE STARTING VECTOR ***| C ----------------------------------- 70 T = 1./SQRT(T) DO 80 I = O,N 80 X(I) = T*X(I) C ------------------------- C |*** START ITERATION ***| C ------------------------- 90 DO 100 I = 1,N 100 Y(I) = X(I) CALL STP(MULT,Y,W,H,M,N,N5) DO 110 I = 1,N 110 Y(I+N) = Y(I) CALL STP(MULT,Y(M),W,H,M,N,N5) 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 120 I = O,N 120 A = A + Y(I+N)**2 IF ( A .EQ. 0. ) GOTO 670 A = 1./SQRT(A) S = 0. DO 130 I = O,N S = S + X(I)*Y(I) J = I + N 130 Y(J) = Y(J)*A A = X(O) B = 1. + ABS(A) IF ( A .LT. 0. ) GOTO 140 A = A + 1. S = S + Y(O) GOTO 150 140 A = A - 1. S = S - Y(O) 150 S = S/B Y(O) = Y(O) - A*S C = 0. D = 0. DO 160 I = P,N T = Y(I) - S*X(I) Y(I) = T C = C + T*X(I) 160 D = D + T*T IF ( D .NE. 0. ) GOTO 170 C = X(P) D = 1. Y(P) = 1. 170 D = 1./SQRT(D) C = C/B Y(O) = -C*A*D DO 180 I = P,N 180 Y(I) = D*(Y(I)-C*X(I)) S = 0. T = 0. DO 190 I = O,N J = I + N S = S + Y(J)*X(I) 190 T = T + Y(J)*Y(I) D = 0. DO 200 I = O,N A = Y(I+N) D = D + (A-S*X(I)-T*Y(I))**2 200 X(I) = A L = L + 2 IF ( L .GE. LIMIT ) GOTO 210 IF ( D .GT. Q ) GOTO 90 C --------------------------------------------------------- C |*** IMPROVE ACCURACY USING ORTHOGONAL STARTING PAIR ***| C --------------------------------------------------------- 210 S = 0. DO 220 I = O,N 220 S = S + X(I)*Y(I) A = X(O) B = 1. + ABS(A) IF ( A .LT. 0. ) GOTO 230 A = A + 1. S = S + Y(O) GOTO 240 230 A = A - 1. S = S - Y(O) 240 S = S/B Y(O) = Y(O) - A*S C = 0. D = 0. DO 250 I = P,N T = Y(I) - S*X(I) Y(I) = T C = C + T*X(I) 250 D = D + T*T IF ( D .NE. 0. ) GOTO 260 C = X(P) D = 1. Y(P) = 1. 260 C = C/B D = 1./SQRT(D) Y(O) = -C*A*D Y(M2) = Y(O) X(M2) = X(O) DO 270 I = P,N J = I + N T = D*(Y(I)-C*X(I)) Y(I) = T Y(J) = T 270 X(J) = X(I) CALL STP(MULT,X,W,H,M,N,N5) CALL STP(MULT,Y,W,H,M,N,N5) IF ( H .EQ. 0 ) GOTO 290 DO 280 I = 1,H X(I) = 0. 280 Y(I) = 0. 290 A = 0. B = 0. C = 0. D = 0. S = 0. T = 0. DO 300 I = O,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) 300 D = D + G*Y(J) IF ( T .EQ. 0. ) GOTO 710 IF ( S .EQ. 0. ) GOTO 760 F = 1./SQRT(S) G = 1./SQRT(T) R = 0. S = 0. T = 0. DO 310 I = O,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) 310 Y(I) = G*Y(I) T = F*F*S + G*G*T S = F*F*R L = L + 2 IF ( L .GE. LIMIT ) GOTO 400 IF ( T .LE. Q ) GOTO 410 IF ( S .GT. Q ) GOTO 210 C -------------------------------- C |*** SINGLE REAL EIGENVALUE ***| C -------------------------------- 320 SIZE = 3*L + 1 330 DIF = SQRT(S) 340 EX = CMPLX(A,0.) EY = (0.,0.) DO 350 I = 1,N 350 X(I+N) = X(I) CALL STP(MULT,X(M),W,H,M,N,N5) CALL RVC(A,X,X(M),W(N3),H) CALL PX(X,W(M),H,N,N5) IF ( O .GE. NS ) GOTO 380 CALL HSE(X(M),W(M),H,O,N) DO 360 I = 1,N 360 Y(I) = 0. Y(O) = 1. N5 = (O*(N+N-P))/2 CALL STP(MULT,Y,W,O,M,N,N5) Y(P) = 0. DO 370 I = 1,P 370 W(N4+I) = Y(I) 380 J = N2 I = N 390 K = J - 1 Y(K) = 0. Y(J) = 0. X(K) = X(I) X(J) = 0. J = J - 2 I = I - 1 IF ( I .GT. 0 ) GOTO 390 NE = NE - 1 NP = NE H = O RETURN C ----------------------------- C |*** ITERATIONS AT LIMIT ***| C ----------------------------- 400 SIZE = 3*L + 2 IF ( S+S .LT. T ) GOTO 330 GOTO 420 410 SIZE = 3*L + 1 420 DIF = SQRT(T) IF ( H .EQ. 0 ) GOTO 440 DO 430 I = M,M1 X(I) = 0. 430 Y(I) = 0. 440 IF ( P .GE. NS ) GOTO 500 CALL HSE(X,W(M),H,O,N) T = 0. DO 450 I = O,N 450 T = T + X(I)*Y(I) DO 460 I = O,N 460 Y(I) = Y(I) - T*X(I) CALL HSE(Y,W(M),O,P,N) DO 470 I = 1,N X(I) = 0. 470 Y(I) = 0. X(O) = 1. N5 = (P*(N+N-P1))/2 CALL STP(MULT,X,W,P,M,N,N5) DO 480 I = 1,P 480 W(I+N4) = X(I) Y(P) = 1. CALL STP(MULT,Y,W,P,M,N,N5) N5 = (H*(N+N-O))/2 Y(P1) = 0. N4 = N4 + P DO 490 I = 1,P1 490 W(I+N4) = Y(I) 500 S = .5*(D-A) T = B*C R = S*S + T IF ( R .LT. 0. ) GOTO 610 C ------------------------------ C |*** TWO REAL EIGENVALUES ***| C ------------------------------ R = ABS(S) + SQRT(R) IF ( S .LT. 0. ) GOTO 510 EX = CMPLX(A-T/R,0.) EY = CMPLX(A+R,0.) GOTO 520 510 EX = CMPLX(A+T/R,0.) EY = CMPLX(A-R,0.) 520 S = CABS(A-EX) + ABS(B) T = CABS(D-EX) + ABS(C) IF ( S .GT. T ) GOTO 540 IF ( T .NE. 0. ) GOTO 530 E = 0. F = 1. A = 1. B = 0. GOTO 580 530 E = C/T F = (EX-D)/T GOTO 550 540 E = (EX-A)/S F = B/S 550 S = CABS(A-EY) + ABS(B) T = CABS(D-EY) + ABS(C) IF ( S .GT. T ) GOTO 570 IF ( T .NE. 0. ) GOTO 560 E = 0. F = 1. A = 1. B = 0. GOTO 580 560 A = C/T B = (EY-D)/T GOTO 580 570 A = (EY-A)/S B = B/S 580 DO 590 I = O,N J = I + N S = F*X(J) + E*Y(J) T = B*X(J) + A*Y(J) X(I) = S Y(I) = T X(J) = S 590 Y(J) = T CALL STP(MULT,X(M),W,H,M,N,N5) T = EX CALL RVC(T,X,X(M),W(N3),H) CALL PX(X,W(M),H,N,N5) CALL STP(MULT,Y(M),W,H,M,N,N5) T = EY CALL RVC(T,Y,Y(M),W(N3),H) CALL PX(Y,W(M),H,N,N5) J = N2 I = N 600 K = J - 1 X(K) = X(I) X(J) = 0. Y(K) = Y(I) Y(J) = 0. I = I - 1 J = J - 2 IF ( I .GT. 0 ) GOTO 600 NE = NE - 2 IF ( NE .LT. 0 ) NE = 0 NP = NE H = P RETURN C ----------------------------- C |*** COMPLEX EIGENVALUES ***| C ----------------------------- 610 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) IF ( S .GT. T ) GOTO 650 B = C/T C = (Q-D)/T A = R/T DO 620 I = O,N J = I + N T = C*X(J) + B*Y(J) S = A*X(J) X(I) = T Y(I) = S X(J) = T 620 Y(J) = S 630 CALL STP(MULT,X(M),W,H,M,N,N5) CALL STP(MULT,Y(M),W,H,M,N,N5) CALL CVC(Q,R,X,Y,X(M),Y(M),W(N3),H) CALL PX(X,W(M),H,N,N5) CALL PX(Y,W(M),H,N,N5) J = N2 I = N 640 K = J - 1 X(J) = Y(I) Y(J) = -Y(I) X(K) = X(I) Y(K) = X(I) J = J - 2 I = I - 1 IF ( I .GT. 0 ) GOTO 640 NE = NE - 2 IF ( NE .LT. 0 ) NE = 0 NP = NE H = P RETURN 650 C = B/S B = (Q-A)/S D = R/S DO 660 I = O,N J = I + N T = C*X(J) + B*Y(J) S = D*Y(J) X(I) = T Y(I) = S X(J) = T 660 Y(J) = S GOTO 630 670 A = 0. DIF = 0. SIZE = 3*L + 1 DO 680 I = O,N 680 IF ( Y(I) .NE. 0. ) GOTO 690 GOTO 340 690 DO 700 I = O,N 700 X(I) = Y(I) GOTO 340 710 IF ( S .NE. 0. ) GOTO 730 DO 720 I = O,N 720 X(I) = X(I+N) GOTO 670 730 S = 1./SQRT(S) A = A*S DO 740 I = O,N 740 X(I) = X(I)*S S = 0. DO 750 I = O,N 750 S = S + (X(I)-A*X(I+N))**2 IF ( S .GT. Q ) GOTO 90 A = 0. GOTO 320 760 T = 1./SQRT(T) DO 770 I = O,N 770 X(I) = T*Y(I) GOTO 90 END C% SUBROUTINE RVC(F,X,Y,V,K) REAL A,B,C,D,F,R,S,X(1),Y(1),V(1) INTEGER I,J,K,L,M,N L = K J = (K*(K+3))/2 - 1 IF ( L .EQ. 0 ) RETURN 10 M = L - 1 D = F - V(J) IF ( M .EQ. 0 ) GOTO 50 B = V(J-L) IF ( B .EQ. 0. ) GOTO 50 A = F - V(J-L-1) C = V(J-1) S = A*D - B*C IF ( S .NE. 0. ) GOTO 30 20 WRITE(6,*) 'MORE THAN TWO EIGENVALUES WITH THE SAME MAGNITUDE' WRITE(6,*) 'HAVE BEEN DETECTED. SUBROUTINE MPOWER MUST STOP' STOP 30 R = (C*Y(L)+D*Y(M))/S S = (B*Y(M)+A*Y(L))/S X(M) = R X(L) = S N = J - L M = N - L J = M - 1 L = L - 2 IF ( L .EQ. 0 ) RETURN DO 40 I = 1,L 40 Y(I) = Y(I) + R*V(M+I) + S*V(N+I) GOTO 10 50 IF ( D .EQ. 0. ) GOTO 20 R = Y(L)/D X(L) = R N = J - L J = N - 1 L = M IF ( L .EQ. 0 ) RETURN DO 60 I = 1,L 60 Y(I) = Y(I) + R*V(N+I) GOTO 10 END C% SUBROUTINE CVC(F,G,W,X,Y,Z,V,K) REAL W(1),X(1),Y(1),Z(1),V(1),A,B,C,D,F,G,P,Q COMPLEX R,S INTEGER I,J,K,L,M,N L = K J = (K*(K+3))/2 - 1 IF ( L .EQ. 0 ) RETURN 10 M = L - 1 D = F - V(J) IF ( M .EQ. 0 ) GOTO 40 B = V(J-L) IF ( B .EQ. 0. ) GOTO 40 A = F - V(J-L-1) C = V(J-1) S = CMPLX(A*D-G*G-B*C,G*(A+D)) IF ( CABS(S) .NE. 0. ) GOTO 20 WRITE(6,*) 'MORE THAN TWO EIGENVALUES WITH THE SAME MAGNITUDE' WRITE(6,*) 'HAVE BEEN DETECTED. SUBROUTINE MPOWER MUST STOP' STOP 20 R = CMPLX(C*Y(L)+D*Y(M)-G*Z(M),C*Z(L)+G*Y(M)+D*Z(M))/S S = CMPLX(B*Y(M)+A*Y(L)-G*Z(L),G*Y(L)+A*Z(L)+B*Z(M))/S A = REAL(R) B = REAL(S) C = AIMAG(R) D = AIMAG(S) W(M) = A W(L) = B X(M) = C X(L) = D N = J - L M = N - L J = M - 1 L = L - 2 IF ( L .EQ. 0 ) RETURN DO 30 I = 1,L P = V(M+I) Q = V(N+I) Y(I) = Y(I) + A*P + B*Q 30 Z(I) = Z(I) + C*P + D*Q GOTO 10 40 R = CMPLX(Y(L),Z(L))/CMPLX(D,G) A = REAL(R) B = AIMAG(R) W(L) = A X(L) = B N = J - L J = N - 1 L = L - 1 IF ( L .EQ. 0 ) RETURN DO 50 I = 1,L Y(I) = Y(I) + A*V(N+I) 50 Z(I) = Z(I) + B*V(N+I) GOTO 10 END C% SUBROUTINE STP(MULT,X,W,K,M,N,N5) REAL X(1),W(1) INTEGER I,K,M,N,N5 CALL PX(X,W(M),K,N,N5) DO 10 I = 1,N 10 W(I) = X(I) CALL MULT(X,W) CALL XP(X,W(M),K,N) RETURN END SUBROUTINE HSE(Y,W,K,L,N) REAL Y(1),W(1),S,T INTEGER I,J,K,L,N J = K*N - (K*(1+K))/2 T = 0. DO 10 I = L,N 10 T = T + Y(I)**2 T = SQRT(T) S = Y(L) IF ( S .LT. 0. ) S = S - T IF ( S .GE. 0. ) S = S + T T = SQRT(ABS(S*T)) IF ( T .NE. 0. ) T = 1./T DO 20 I = L,N Y(I) = T*Y(I) 20 W(I+J) = Y(I) Y(L) = S*T W(L+J) = Y(L) RETURN END SUBROUTINE PX(Y,Z,K,N,N5) REAL Y(1),Z(1),T INTEGER I,J,K,L,M,N,N5 IF ( K .EQ. 0 ) RETURN L = N5 DO 20 M = 1,K J = K - M + 1 L = L - N + J T = 0. DO 10 I = J,N 10 T = T + Z(L+I)*Y(I) DO 20 I = J,N 20 Y(I) = Y(I) - T*Z(L+I) RETURN END SUBROUTINE XP(Y,Z,K,N) REAL Y(1),Z(1),T INTEGER I,J,K,L,N IF ( K .EQ. 0. ) RETURN L = 0 DO 30 J = 1,K T = 0. DO 10 I = J,N 10 T = T + Z(L+I)*Y(I) DO 20 I = J,N 20 Y(I) = Y(I) - T*Z(L+I) 30 L = L + N - J RETURN END