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