      subroutine eigen(isize,a,evout,vout)

!  compute eigenvalues and eigenvectors for a square real matrix

!   inputs
!     isize - dimension of input square matrix
!     a     - input real matrix
!
!   outputs
!     evout - vector of eigenvalues, length isize
!     vout  - matrix of eigenvectors, columnwise, isize x isize

      real*4 a(isize,isize),ev(isize),v(isize,isize)
      real *4 wi(isize),zi(isize,isize)
      integer ivec(isize)
      real*4 evout(isize),vout(isize,isize)
!      real*4 q(3),asave(3,3)

c$$$      print *,'a'
c$$$      do i=1,isize
c$$$         print *,'q ',(a(i,j),j=1,3)
c$$$      end do
c$$$      asave=a

      call eigenp(isize,isize,a,15.,ev,wi,v,zi,ivec)

c$$$      print *
c$$$      print *,'q ',ev
c$$$      print *
c$$$      do i=1,isize
c$$$         print *,'q ',(v(i,j),j=1,3)
c$$$      end do
c$$$
c$$$      print *
c$$$      do j=1,3
c$$$      q(1)=v(1,j)
c$$$      q(2)=v(2,j)
c$$$      q(3)=v(3,j)
c$$$
c$$$      print *,'consistency ',matmul(asave,q),ev(j)*q
c$$$      end do
c$$$      print *

!  sort and make eigenvalues positive
      do k=1,isize

      emax=-1.
      do i=1,isize
         if(abs(ev(i)).ge.emax)then
            emax=abs(ev(i))
            imax=i
         end if
      end do
      evout(k)=ev(imax)
      do j=1,isize
         vout(j,k)=v(j,imax)
      end do
      ev(imax)=0.

      end do

      return
      end


      SUBROUTINE EIGENP ( N, NM, A, T, EVR, EVI, VECR, VECI, INDIC )

c*********************************************************************72

      IMPLICIT NONE

      INTEGER NM

      REAL A(NM,1)
      DOUBLE PRECISION D1
      DOUBLE PRECISION D2
      DOUBLE PRECISION D3
      REAL ENORM
      REAL EPS
      REAL EVI(NM)
      REAL EVR(NM)
      REAL EX
      INTEGER I
      INTEGER INDIC(NM)
      INTEGER IVEC
      INTEGER IWORK(100)
      INTEGER J
      INTEGER K
      INTEGER K1
      INTEGER KON
      INTEGER L
      INTEGER L1
      INTEGER LOCAL(100)
      INTEGER M
      INTEGER N
      DOUBLE PRECISION PRFACT(100)
      REAL R
      REAL R1
      REAL SUBDIA(100)
      REAL T
      REAL VECI(NM,1)
      REAL VECR(NM,1)
      REAL WORK(100)
      REAL WORK1(100)
      REAL WORK2(100)
C
C  THIS SUBROUTINE FINDS ALL THE EIGENVALULES AND THE
C  EIGENVECTORS OF A REAL GENERAL MATRIX OF ORDER N.
C
C  FIRST IN THE SUBROUTINE SCALE THE MATRIX IS SCALED SO THAT
C  THE CORRESPONDING ROWS AND COLUMNS ARE APPROXIMATELY 
C  BALANCED, AND THEN THE MATRIX IS NORMALIZED SO THAT THE
C  VALUE OF THE EUCLIDEAN NORM OF THE MATRIX IS EQUAL TO ONE.
C
C  THE EIGENVALUES ARE COMPUTED BY THE QR DOUBLE-STEP METHOD
C  IN THE SUBROUTINE HESQR.
C
C  THE EIGENVECTORS ARE COMPUTED BY INVERSE ITERATION IN
C  THE SUBROUTINE REALVE, FOR THE REAL EIGENVALUES, OR IN THE
C  SUBROUTINE COMPVE, FOR THE COMPLEX EIGENVALUES.
C
C  THE ELEMENTS OF THE MATRIX ARE TO BE STORED IN THE FIRST N
C  ROWS AND COLUMNS OF THE TWO DIMENSIONAL ARRAY A.  THE
C  ORIGINAL MATRIX IS DESTROYED BY THE SUBROUTINE.
C
C  N IS THE ORDER OF THE MATRIX.
C
C  NM DEFINES THE FIRST DIMENSION OF THE TWO DIMENSIONAL
C  ARRAYS A, VECR, VECI, AND THE DIMENSION OF THE ONE
C  DIMENSIONAL ARRAYS EVR, EVI AND INDIC.  THEREFORE, THE
C  CALLING PROGRAM SHOULD CONTAIN THE FOLLOWING DECLARATIONS:
C    REAL A(NM,NN)
C    REAL EVI(NM)
C    REAL EVR(NM)
C    INTEGER INDIC(NM)
C    REAL VECI(NM,NN)
C    REAL VECR(NM,NN)
C  WHERE NM AND NN ARE ANY NUMBERS EQUAL TO OR GREATER THAN N.
C  THE UPPER LIMIT FOR NM IS EQUAL TO 100, BUT MAY BE 
C  INCREASED TO THE VALUE MAX BY REPLACING THE DIMENSION
C  STATEMENTS:
C    INTEGER IWORK(100)
C    INTEGER LOCAL(100)
C    DOUBLE PRECISION PRFACT(100)
C    REAL SUBDIA(100)
C    REAL WORK(100)
C    REAL WORK1(100)
C    REAL WORK2(100)
C  IN THE SUBROUTINE EIGENP BY
C    INTEGER IWORK(MAX)
C    INTEGER LOCAL(MAX)
C    DOUBLE PRECISION PRFACT(MAX)
C    REAL SUBDIA(MAX)
C    REAL WORK(MAX)
C    REAL WORK1(MAX)
C    REAL WORK2(MAX)
C  NM AND NN ARE OF COURSE BOUNDED BY THE SIZE OF THE STORE.
C
C  THE REAL PARAMETER T MUST BE SET EQUAL TO THE NUMBER OF
C  BINARY DIGITS IN THE MANTISSA OF A SINGLE PRECISION
C  FLOATING-POINT NUMBER.
C
C  THE REAL PARTS OF THE N COMPUTED EIGENVALUES WILL BE FOUND
C  IN THE FIRST N PLACES OF THE ARRAY EVI.
C  THE REAL COMPONENTS OF THE NORMALIZED EIGENVECTOR I
C  (I=1,2,...,N) CORRESPONDING TO THE EIGENVALUE STORED IN
C  EVR(I) AND EVI(I) WILL BE FOUND IN THE FIRST N PLACES OF
C  THE COLUMN I OF THE TWO DIMENSIONAL ARRAY VECR AND THE
C  IMAGINARY COMPONENTS IN THE FIRST N PLACES OF THE COLUMN I
C  OF THE TWO DIMENSIONAL ARRAY VECI.
C
C  THE REAL EIGENVECTOR IS NORMALIZED SO THAT THE SUM OF THE
C  SQUARES OF THE COMPONENTS IS EQUAL TO ONE.
C
C  THE COMPLEX EIGENVECTOR IS NORMALIZED SO THAT THE
C  COMPONENT WITH THE LARGEST VALUE IN MODULUS HAS ITS REAL
C  PART EQUAL TO ONE AND THE IMAGINARY PART EQUAL TO ZERO.
C
C  THE ARRAY INDIC INDICATES THE SUCCESS OF THE SUBROUTINE
C  EIGENP AS FOLLOWS:
C
C    VALUE OF INDIC(I)    EIGENVALUE I    EIGENVECTOR I
C
C           0             NOT FOUND       NOT FOUND
C           1                 FOUND       NOT FOUND
C           2                 FOUND           FOUND
C
      IF ( N .NE. 1 ) GO TO 1
      EVR(1) = A(1,1)
      EVI(1) = 0.0E+00
      VECR(1,1) = 1.0E+00
      VECI(1,1) = 0.0E+00
      INDIC(1) = 2
      GO TO 25

1     CALL SCALE ( N, NM, A, VECI, PRFACT, ENORM )
C
C  THE COMPUTATION OF THE EIGENVALUES OF THE NORMALIZED
C  MATRIX.
C
      EX = EXP ( - T * ALOG ( 2.0E+00 ) )
      CALL HESQR ( N, NM, A, VECI, EVR, EVI, SUBDIA, INDIC, EPS, EX )
C
C  THE POSSIBLE DECOMPOSITION OF THE UPPER-HESSENBERG MATRIX
C  INTO THE SUBMATRICESS OF LOWER ORDER IS INDICATED IN THE
C  ARRAY LOCAL.  THE DECOMPOSITION OCCURS WHEN SOME
C  SUBDIAGONAL ELEMENTS ARE IN MODULUS LESS THAN A SMALL
C  POSITIVE NUMBER EPS DEFINED IN THE SUBROUTINE HESQR.  THE
C  AMOUNT OF WORK IN THE EIGENVECTOR PROBLEM MAY BE
C  DIMINISHED IN THIS WAY.
C  
      J = N
      I = 1
      LOCAL(1) = 1
      IF ( J .EQ. 1 ) GO TO 4
2     IF ( ABS ( SUBDIA(J-1) ) .GT. EPS ) GO TO 3
      I = I + 1
      LOCAL(I) = 0
3     J = J - 1
      LOCAL(I) = LOCAL(I) + 1
      IF ( J .NE. 1 ) GO TO 2
C
C  THE EIGENVECTOR PROBLEM.
C
4     K = 1
      KON = 0
      L = LOCAL(1)
      M = N
      DO 10 I = 1, N
        IVEC = N - I + 1
        IF ( I .LE. L ) GO TO 5
        K = K + 1
        M = N - L
        L = L + LOCAL(K)
5       IF ( INDIC(IVEC) .EQ. 0 ) GO TO 10
        IF ( EVI(IVEC) .NE. 0.0E+00 ) GO TO 8
C
C  TRANSFER OF AN UPPER-HESSENBERG MATRIX OF ORDER M FROM
C  THE ARRAYS VECI AND SUBDIA INTO THE ARRAY A.
C
        DO 7 K1 = 1, M
          DO 6 L1 = K1, M
6           A(K1,L1) = VECI(K1,L1)
          IF ( K1 .EQ. 1 ) GO TO 7
          A(K1,K1-1) = SUBDIA(K1-1)
7       CONTINUE
C
C  THE COMPUTATION OF THE REAL EIGENVECTOR IVEC OF THE UPPER-
C  HESSENBERG MATRIX CORRESPONDING TO THE REAL EIGENVALUE
C  EVR(IVEC).
C
        CALL REALVE ( N, NM, M, IVEC, A, VECR, EVR, EVI, IWORK,
     &    WORK, INDIC, EPS, EX )
        GO TO 10
C
C  THE COMPUTATION OF THE COMPLEX EIGENVECTOR IVEC OF THE
C  UPPER-HESSENBERG MATRIX CORRESPONDING TO THE COMPLEX
C  EIGENVALUE EVR(IVEC)+I*EVI(IVEC).  IF THE VALUE OF KON IS
C  NOT EQUAL TO ZERO, THEN THIS COMPLEX EIGENVECTOR HAS
C  ALREADY BEEN FOUND FROM ITS CONJUGATE.
C
8       IF ( KON .NE. 0 ) GO TO 9
        KON = 1
        CALL COMPVE ( N, NM, M, IVEC, A, VECR, VECI, EVR, EVI, INDIC,
     &    IWORK, SUBDIA, WORK1, WORK2, WORK, EPS, EX )
        GO TO 10
9       KON = 0
10    CONTINUE
C
C  THE RECONSTRUCTION OF THE MATRIX USED IN THE REDUCTION OF
C  MATRIX A TO AN UPPER-HESSENBERG FORM BY HOUSEHOLDER METHOD.
C
      DO 12 I = 1, N
        DO 11 J = I, N
          A(I,J) = 0.0E+00
11        A(J,I) = 0.0E+00
12      A(I,I) = 1.0E+00

      IF ( N .LE. 2 ) GO TO 15
      M = N - 2
      DO 14 K = 1, M
        L = K + 1
        DO 14 J = 2, N
          D1 = 0.0E+00
          DO 13 I = L, N
            D2 = VECI(I,K)
13          D1 = D1 + D2 * A(J,I)
          DO 14 I = L, N
14          A(J,I) = A(J,I) - VECI(I,K) * D1
C
C  THE COMPUTATION OF THE EIGENVECTORS OF THE ORIGINAL NON-
C  SCALED MATRIX.
C
15    KON = 1
      DO 24 I = 1, N
        L = 0
        IF ( EVI(I) .EQ. 0.0E+00 ) GO TO 16
        L = 1
        IF ( KON .EQ. 0 ) GO TO 16
        KON = 0
        GO TO 24
16      DO 18 J = 1, N
          D1 = 0.0E+00
          D2 = 0.0E+00
          DO 17 K = 1, N
            D3 = A(J,K)
            D1 = D1 + D3 * VECR(K,I)
            IF ( L .EQ. 0 ) GO TO 17
            D2 = D2 + D3 * VECR(K,I-1)
17        CONTINUE
          WORK(J) = D1 / PRFACT(J)
          IF ( L .EQ. 0 ) GO TO 18
          SUBDIA(J) = D2 / PRFACT(J)
18      CONTINUE
C
C  THE NORMALIZATION OF THE EIGENVECTORS AND THE COMPUTATION
C  OF THE EIGENVALUES OF THE ORIGINAL NON-NORMALIZED MATRIX.
C
        IF ( L .EQ. 1 ) GO TO 21
        D1 = 0.0E+00
        DO 19 M = 1, N
19        D1 = D1 + WORK(M)**2
        D1 = SQRT ( D1 )

        DO 20 M = 1, N
          VECI(M,I) = 0.0E+00
20        VECR(M,I) = WORK(M) / D1

        EVR(I) = EVR(I) * ENORM
        GO TO 24

21      KON = 1
        EVR(I) = EVR(I) * ENORM
        EVR(I-1) = EVR(I)
        EVI(I) = EVI(I) * ENORM
        EVI(I-1) = -EVI(I)

        R = 0.0E+00
        DO 22 J = 1, N
          R1 = WORK(J)**2 + SUBDIA(J)**2
          IF ( R .GE. R1 ) GO TO 22
          R = R1
          L = J
22      CONTINUE

        D3 = WORK(L)
        R1 = SUBDIA(L)
        DO 23 J = 1, N
          D1 = WORK(J)
          D2 = SUBDIA(J)
          VECR(J,I) = ( D1 * D3 + D2 * R1 ) / R
          VECI(J,I) = ( D2 * D3 - D1 * R1 ) / R
          VECR(J,I-1) = VECR(J,I)
23        VECI(J,I-1) = -VECI(J,I)

24    CONTINUE

25    RETURN
      END
      SUBROUTINE SCALE ( N, NM, A, H, PRFACT, ENORM )

c*********************************************************************72

      IMPLICIT NONE

      INTEGER NM

      REAL A(NM,1)
      REAL BOUND1
      REAL BOUND2
      DOUBLE PRECISION COLUMN
      REAL ENORM
      DOUBLE PRECISION FACTOR
      DOUBLE PRECISION FNORM
      REAL H(NM,1)
      INTEGER I
      INTEGER ITER
      INTEGER J
      INTEGER N
      INTEGER NCOUNT
      DOUBLE PRECISION PRFACT(NM)
      DOUBLE PRECISION Q
      DOUBLE PRECISION ROW
C
C  THIS SUBROUTINE STORES THE MATRIX OF THE ORDER N FROM THE
C  ARRAY A INTO THE ARRAY H.  AFTERWARD, THE MATRIX IN THE
C  ARRAY A IS SCALED SO THAT THE QUOTIENT OF THE ABSOLUTE SUM
C  OF THE OFF-DIAGONAL ELEMENTS OF COLUMN I AND THE ABSOLUTE
C  SUM OF THE OFF-DIAGONAL ELEMENTS OF ROW I LIES WITHIN THE
C  VALUES OF BOUND1 AND BOUND2.
C
C  THE COMPONENT I OF THE EIGENVECTOR OBTAINED BY USING THE
C  SCALED MATRIX MUST BE DIVIDED BY THE VALUE FOUND IN THE
C  PRFACT(I) OF THE ARRAY PRFACT.  IN THIS WAY, THE EIGENVECTOR
C  OF THE NON-SCALED MATRIX IS OBTAINED.
C
C  AFTER THE MATRIX IS SCALED, IT IS NORMALIZED SO THAT THE
C  VALUE OF THE EUCLIDEAN NORM IS EQUAL TO ONE.
C  IF THE PROCESS OF SCALING WAS NOT SUCCESSFUL, THE ORIGINAL
C  MATRIX FROM THE ARRAY H WOULD BE STORED BAK INTO A AND
C  THE EIGENPROBLEM WOULD BE SOLVED BY USING THIS MATRIX.
C
C  NM DEFINES THE FIRST DIMENSION OF THE ARRAYS A AND H.  NM
C  MUST BE GREATER OR EQUAL TO N.
C
C  THE EIGENVALUES OF THE NORMALIZED MATRIX MUST BE 
C  MULTIPLIED BY THE SCALAR ENORM IN ORDER THAT THEY BECOME
C  THE EIGENVALUES OF THE NON-NORMALIZED MATRIX.
C
      DO 2 I = 1, N
        DO 1 J = 1, N
1         H(I,J) = A(I,J)
2       PRFACT(I) = 1.0E+00

      BOUND1 = 0.75E+00
      BOUND2 = 1.33E+00
      ITER = 0
3     NCOUNT = 0

      DO 8 I = 1, N

        COLUMN = 0.0E+00
        ROW = 0.0E+00
        DO 4 J = 1, N
          IF ( I .EQ. J ) GO TO 4
          COLUMN = COLUMN + ABS ( A(J,I) )
          ROW = ROW + ABS ( A(I,J) )
4       CONTINUE

        IF ( COLUMN .EQ. 0.0E+00 ) GO TO 5
        IF ( ROW .EQ. 0.0E+00 ) GO TO 5
        Q = COLUMN / ROW
        IF ( Q .LT. BOUND1 ) GO TO 6
        IF ( Q .GT. BOUND2 ) GO TO 6
5       NCOUNT = NCOUNT + 1
        GO TO 8
6       FACTOR = DSQRT ( Q )

        DO 7 J = 1, N
          IF ( I .EQ. J ) GO TO 7
          A(I,J) = A(I,J) * FACTOR
          A(J,I) = A(J,I) / FACTOR
7       CONTINUE

        PRFACT(I) = PRFACT(I) * FACTOR

8     CONTINUE

      ITER = ITER + 1
      IF ( ITER .GT. 30 ) GO TO 11
      IF ( NCOUNT .LT. N ) GO TO 3

      FNORM = 0.0E+00
      DO 9 I = 1, N
        DO 9 J = 1, N
          Q = A(I,J)
9         FNORM = FNORM + Q * Q
      FNORM = DSQRT ( FNORM )

      DO 10 I = 1, N
        DO 10 J = 1, N
10        A(I,J) = A(I,J) / FNORM

      ENORM = FNORM
      GO TO 13

11    DO 12 I = 1, N
        DO 12 J = 1, N
12        A(I,J) = H(I,J)

      ENORM = 1.0E+00

13    RETURN
      END
      SUBROUTINE HESQR ( N, NM, A, H, EVR, EVI, SUBDIA, INDIC, EPS, EX )

c*********************************************************************72

      IMPLICIT NONE

      INTEGER NM

      REAL A(NM,1)
      REAL EPS
      REAL EX
      REAL H(NM,1)
      REAL EVI(NM)
      REAL EVR(NM)
      INTEGER I
      INTEGER INDIC(NM)
      INTEGER J
      INTEGER K
      INTEGER L
      INTEGER M
      INTEGER MAXST
      INTEGER M1
      INTEGER N
      INTEGER NS
      REAL R
      DOUBLE PRECISION S
      REAL SHIFT
      DOUBLE PRECISION SR
      DOUBLE PRECISION SR2
      REAL SUBDIA(NM)
      REAL T
      DOUBLE PRECISION X
      DOUBLE PRECISION Y
      DOUBLE PRECISION Z
C
C  THIS SUBROUTINE FINDS ALL THE EIGENVALUES OF A REAL
C  GENERAL MATRIX.  THE ORIGINAL MATRIX A OF ORDER N IS
C  REDUCED TO THE UPPER-HESSENBERG FORM H BY MEANS OF
C  SIMILARITY TRANSFORMATIONS (HOUSEHOLDER METHOD).  THE MATRIX
C  H IS PRESERVED IN THE UPPER HALF OF THE ARRAY H AND IN THE
C  ARRAY SUBDIA.  THE SPECIAL VECTORS USED IN THE DEFINITION
C  OF THE HOUSEHOLDER TRANSFORMATION MATRICES ARE STORED IN
C  THE LOWER PART OF THE ARRAY H.
C
C  NM IS THE FIRST DIMENSION OF THE ARRAYS A AND H.  NM MUST
C  BE EQUAL TO OR GREATER THAN N.
C
C  THE REAL PARTS OF THE N EIGENVALUES WILL BE FOUND IN THE
C  FIRST N PLACES OF THE ARRAY EVR, AND
C  THE IMAGINARY PARTS IN THE FIRST N PLACES OF THE ARRAY EVI.
C
C  THE ARRAY INDIC INDICATES THE SUCCESS OF THE ROUTINE AS
C  FOLLOWS:
C
C    VALUE OF INDIC(I)  EIGENVALUE I
C           0             NOT FOUND
C           1                 FOUND
C
C  EPS IS A SMALL POSITIVE NUMBER THAT NUMERICALLY REPRESENTS
C  ZERO IN THE PROGRAM.  EPS = (EUCLIDEAN NORM OF H)*EX, WHERE
C  EX = 2**(-T).  T IS THE NUMBER OF BINARY DIGITS IN THE
C  MANTISSA OF A FLOATING POINT NUMBER.
C
C  REDUCTION OF THE MATRIX A TO AN UPPER-HESSENBERG FORM H.
C  THERE ARE N-2 STEPS.
C
!      IF ( N - 2 ) 14, 1, 2
      if(n-2.lt.0)go to 14
      if(n-2.eq.0)go to 1
      if(n-2.gt.0)go to 2

1     SUBDIA(1) = A(2,1)
      GO TO 14

2     M = N - 2

      DO 12 K = 1, M

        L = K + 1

        S = 0.0E+00
        DO 3 I = L, N
          H(I,K) = A(I,K)
3         S = S + ABS ( A(I,K) )

        IF ( S .NE. ABS ( A(K+1,K) ) ) GO TO 4
        SUBDIA(K) = A(K+1,K)
        H(K+1,K) = 0.0E+00
        GO TO 12
4       SR2 = 0.0E+00

        DO 5 I = L, N
          SR = A(I,K)
          SR = SR / S
          A(I,K) = SR
5         SR2 = SR2 + SR * SR

        SR = DSQRT ( SR2 )
        IF ( A(L,K) .LT. 0.0 ) GO TO 6
        SR = -SR
6       SR2 = SR2 - SR * A(L,K)
        A(L,K) = A(L,K) - SR
        H(L,K) = H(L,K) - SR * S
        SUBDIA(K) = SR * S
        X = S * DSQRT ( SR2 )

        DO 7 I = L, N
          H(I,K) = H(I,K) / X
7         SUBDIA(I) = A(I,K) / SR2
C
C  PREMULTIPLICATION BY THE MATRIX PR.
C
        DO 9 J = L, N
          SR = 0.0E+00
          DO 8 I = L, N
8           SR = SR + A(I,K) * A(I,J)
          DO 9 I = L, N
9           A(I,J) = A(I,J) - SUBDIA(I) * SR
C
C  POSTMULTIPLICATION BY THE MATRIX PR.
C
        DO 11 J = 1, N
          SR = 0.0E+00
          DO 10 I = L, N
10          SR = SR + A(J,I) * A(I,K)
          DO 11 I = L, N
11          A(J,I) = A(J,I) - SUBDIA(I) * SR

12    CONTINUE

      DO 13 K = 1, M
13      A(K+1,K) = SUBDIA(K)
C
C  TRANSFER OF THE UPPER HALF OF THE MATRIX A INTO THE
C  ARRAY H, AND THE CALCULATION OF THE SMALL POSITIVE NUMBER EPS.
C
      SUBDIA(N-1) = A(N,N-1)
14    EPS = 0.0E+00
      DO 15 K = 1, N
        INDIC(K) = 0
        IF ( K .NE. N ) EPS = EPS + SUBDIA(K)**2
        DO 15 I = K, N
          H(K,I) = A(K,I)
15        EPS = EPS + A(K,I)**2

      EPS = EX * SQRT ( EPS )
C
C  THE QR ITERATIVE PROCESS.  THE UPPER-HESSENBERG MATRIX H IS
C  REDUCED TO THE UPPER-MODIFIED TRIANGULAR FORM.
C
C  DETERMINATION OF THE SHIFT OF ORIGIN FOR THE FIRST STEP OF
C  THE QR ITERATIVE PROCESS.
C
      SHIFT = A(N,N-1)
      IF ( N .LE. 2 ) SHIFT = 0.0E+00
      IF ( A(N,N) .NE. 0.0E+00 ) SHIFT = 0.0E+00
      IF ( A(N-1,N) .NE. 0.0E+00 ) SHIFT = 0.0E+00
      IF ( A(N-1,N-1) .NE. 0.0E+00 ) SHIFT = 0.0E+00
      M = N
      NS = 0
      MAXST = N * 10
C
C  TESTING IF THE UPPER HALF OF THE MATRIX IS EQUAL TO ZERO.
C  IF IT IS EQUAL TO ZERO, THE QR PROCESS IS NOT NECESSARY.
C
      DO 16 I = 2, N
        DO 16 K = I, N
          IF ( A(I-1,K) .NE. 0.0E+00 ) GO TO 18
16        CONTINUE

      DO 17 I = 1, N
        INDIC(I) = 1
        EVR(I) = A(I,I)
17      EVI(I) = 0.0E+00

      GO TO 37
C
C  START THE MAIN LOOP OF THE QR PROCESS.
C
18    K = M - 1
      M1 = K
      I = K
C
C  FIND ANY DECOMPOSITIONS OF THE MATRIX.
C  JUMP TO 34 IF THE LAST SUBMATRIX OF THE DECOMPOSITION IS
C  OF THE ORDER ONE.
C  JUMP TO 35 IF THE LAST SUBMATRIX OF THE DECOMPOSITION IS
C  OF THE ORDER TWO.
C
!      IF ( K ) 37, 34, 19
      if(k.lt.0)go to 37
      if(k.eq.0)go to 34
      if(k.gt.0)go to 19

19    IF ( ABS ( A(M,K) ) .LE. EPS ) GO TO 34
      IF ( M - 2 .EQ. 0 ) GO TO 35
20    I = I - 1
      IF ( ABS ( A(K,I) ) .LE. EPS ) GO TO 21
      K = I
      IF ( K .GT. 1 ) GO TO 20
21    IF ( K .EQ. M1 ) GO TO 35
C
C  TRANSFORMATION OF THE MATRIX OF THE ORDER GREATER THAN TWO.
C
      S = A(M,M) + A(M1,M1) + SHIFT
      SR = A(M,M) * A(M1,M1) - A(M,M1)*A(M1,M) + 0.25E+00 * SHIFT**2
      A(K+2,K) = 0.0E+00
C
C  CALCULATE X1, Y1, Z1 FOR THE SUBMATRIX OBTAINED BY THE
C  DECOMPOSITION.
C
      X = A(K,K) * ( A(K,K) - S ) + A(K,K+1) * A(K+1,K) + SR
      Y = A(K+1,K) * ( A(K,K) + A(K+1,K+1) - S )
      R = DABS ( X ) + DABS ( Y )
      IF ( R .EQ. 0.0E+00 ) SHIFT = A(M,M-1)
      IF ( R .EQ. 0.0E+00 ) GO TO 21
      Z = A(K+2,K+1) * A(K+1,K)
      SHIFT = 0.0E+00
      NS = NS + 1
C
C  THE LOOP FOR ONE STEP OF THE QR PROCESS.
C
      DO 33 I = K, M1
        IF ( I .EQ. K ) GO TO 22
C
C  CALCULATE XR, YR, ZR.
C
        X = A(I,I-1)
        Y = A(I+1,I-1)
        Z = 0.0E+00
        IF ( I + 2 .GT. M ) GO TO 22
        Z = A(I+2,I-1)
22      SR2 = DABS ( X ) + DABS ( Y ) + DABS ( Z )
        IF ( SR2 .EQ. 0.0E+00 ) GO TO 23
        X = X / SR2
        Y = Y / SR2
        Z = Z / SR2
23      S = DSQRT ( X * X + Y * Y + Z * Z )
        IF ( X .LT. 0.0E+0 ) GO TO 24
        S = -S
24      IF ( I .EQ. K ) GO TO 25
        A(I,I-1) = S * SR2
25      IF ( SR2 .NE. 0.0E+00 ) GO TO 26
        IF ( I + 3 .GT. M ) GO TO 33
        GO TO 32
26      SR = 1.0E+00 - X / S
        S = X - S
        X = Y / S
        Y = Z / S
C
C  PREMULTIPLICATION BY THE MATRIX PR.
C
        DO 28 J = I, M
          S = A(I,J) + A(I+1,J) * X
          IF ( I + 2 .GT. M ) GO TO 27
          S = S + A(I+2,J) * Y
27        S = S * SR
          A(I,J) = A(I,J) - S
          A(I+1,J) = A(I+1,J) - S * X
          IF ( I + 2 .GT. M ) GO TO 28
          A(I+2,J) = A(I+2,J) - S * Y
28      CONTINUE
C
C  POSTMULTIPLICATION BY THE MATRIX PR.
C
        L = I + 2
        IF ( I .LT. M1 ) GO TO 29
        L = M

29      DO 31 J = K, L
          S = A(J,I) + A(J,I+1) * X
          IF ( I + 2 .GT. M ) GO TO 30
          S = S + A(J,I+2) * Y
30        S = S * SR
          A(J,I) = A(J,I) - S
          A(J,I+1) = A(J,I+1) - S * X
          IF ( I + 2 .GT. M ) GO TO 31
          A(J,I+2) = A(J,I+2) - S * Y
31      CONTINUE

        IF ( I + 3 .GT. M ) GO TO 33
        S = -A(I+3,I+2) * Y * SR
32      A(I+3,I) = S
        A(I+3,I+1) = S * X
        A(I+3,I+2) = S * Y + A(I+3,I+2)
33    CONTINUE

      IF ( NS .GT. MAXST ) GO TO 37
      GO TO 18
C
C  COMPUTE THE LAST EIGENVALUE.
C
34    EVR(M) = A(M,M)
      EVI(M) = 0.0E+00
      INDIC(M) = 1
      M = K
      GO TO 18
C
C  COMPUTE THE EIGENVALUES OF THE LAST 2X2 MATRIX OBTAINED BY
C  THE DECOMPOSITION.
C
35    R = 0.5E+00 * ( A(K,K) + A(M,M) )
      S = 0.5E+00 * ( A(M,M) - A(K,K) )
      S = S * S + A(K,M) * A(M,K)
      INDIC(K) = 1
      INDIC(M) = 1
      IF ( S .LT. 0.0E+00 ) GO TO 36
      T = DSQRT ( S )
      EVR(K) = R - T   
      EVR(M) = R + T
      EVI(K) = 0.0E+00
      EVI(M) = 0.0E+00
      M = M - 2
      GO TO 18

36    T = DSQRT ( -S )
      EVR(K) = R
      EVI(K) = T
      EVR(M) = R
      EVI(M) = -T
      M = M - 2
      GO TO 18

37    RETURN
      END
      SUBROUTINE REALVE ( N, NM, M, IVEC, A, VECR, EVR, EVI,
     &  IWORK, WORK, INDIC, EPS, EX )

c*********************************************************************72

      IMPLICIT NONE

      INTEGER NM

      REAL A(NM,1)
      REAL BOUND
      REAL EPS
      REAL EVALUE
      REAL EVI(NM)
      REAL EVR(NM)
      REAL EX
      INTEGER I
      INTEGER INDIC(NM)
      INTEGER ITER
      INTEGER IVEC
      INTEGER IWORK(NM)
      INTEGER J
      INTEGER K
      INTEGER L
      INTEGER M
      INTEGER N
      INTEGER NS
      REAL PREVIS
      REAL R
      REAL R1
      DOUBLE PRECISION S
      DOUBLE PRECISION SR
      REAL T
      REAL VECR(NM,1)
      REAL WORK(NM)
C
C  THIS SUBROUTINE FINDS THE REAL EIGENVECTOR OF THE REAL
C  UPPER-HESSENBERG MATRIX IN THE ARRAY A, CORRESPONDING TO
C  THE REAL EIGENVALUE STORED IN EVR(IVEC).  THE INVERSE
C  ITERATION METHOD IS USED.
C
C  NOTE THE MATRIX IN A IS DESTROYED BY THE SUBROUTINE.
C
C  N IS THE ORDER OF THE UPPER HESSENBERG MATRIX.
C
C  NM DEFINES THE FIRST DIMENSION OF THE TWO DIMENSIONAL
C  ARRAYS A AND VECR.  NM MUST BE EQUAL TO OR GREATER THAN N.
C
C  M IS THE ORDER OF THE SUBMATRIX OBTAINED BY A SUITABLE
C  DECOMPOSITION OF THE UPPER-HESSENBERG MATRIX IF SOME
C  SUBDIAGONAL ELEMENTS ARE EQUAL TO ZERO.  THE VALUE OF M IS
C  CHOSEN SO THAT THE LAST N-M COMPONENTS OF THE EIGENVECTOR
C  ARE ZERO.
C
C  IVEC GIVES THE POSITION OF THE EIGENVALUE IN THE ARRAY EVR
C  FOR WHICH THE CORRESPONDING EIGENVECTOR IS COMPUTED.
C  
C  THE ARRAY EVI WOULD CONTAIN THE IMAGINARY PARTS OF THE N
C  EIGENVALUES IF THEY EXISTED.
C
C  THE M COMPONENTS OF THE COMPUTED REAL EIGENVECTOR WILL BE
C  FOUND IN THE FIRST M PLACES OF THE COLUMN IVEC OF THE TWO
C  DIMENSIONAL ARRAY VECR.
C
C  IWORK AND WORK ARE THE WORKING STORES USED DURING THE
C  GAUSSIAN ELIMINATION AND BACKSUBSTITUTION PROCESS.
C
C  THE ARRAY INDIC INDICATES THE SUCCESS OF THE ROUTINE AS
C  FOLLOWS:
C    VALUE OF INDIC(I)  EIGENVECTOR I
C           1             NOT FOUND
C           2                 FOUND
C
C  EPS IS A SMALL POSITIVE NUMBER THAT NUMERICALLY REPRESENTS
C  ZERO IN THE PROGRAM.  EPS = (EUCLIDEAN NORM OF A)*EX, WHERE
C  EX = 2**(-T), T IS THE NUMBER OF BINARY DIGITS IN THE
C  MANTISSA OF A FLOATING POINT NUMBER.
C
      VECR(1,IVEC) = 1.0E+00
      IF ( M .EQ. 1 ) GO TO 24
C
C  SMALL PERTURBATION OF EQUAL EIGENVALUES TO OBTAIN A FULL
C  SET OF EIGENVECTORS.
C
      EVALUE = EVR(IVEC)
      IF ( IVEC .EQ. M ) GO TO 2
      K = IVEC + 1
      R = 0.0E+00
      DO 1 I = K, M
        IF ( EVALUE .NE. EVR(I) ) GO TO 1
        IF ( EVI(I) .NE. 0.0E+00 ) GO TO 1
        R = R + 3.0E+00
1     CONTINUE
      EVALUE = EVALUE + R * EX

2     DO 3 K = 1, M
3       A(K,K) = A(K,K) - EVALUE
C
C  GAUSSIAN ELIMINATION OF THE UPPER-HESSENBERG MATRIX A.  ALL
C  ROW INTERCHANGES ARE INDICATED IN THE ARRAY IWORK.  ALL THE
C  MULTIPLIERS ARE STORED AS THE SUBDIAGONAL ELEMENTS OF A.
C
      K = M - 1
      DO 8 I = 1, K
        L = I + 1
        IWORK(I) = 0
        IF ( A(I+1,I) .NE. 0.0E+00 ) GO TO 4
        IF ( A(I,I) .NE. 0.0E+00 ) GO TO 8
        A(I,I) = EPS
        GO TO 8
4       IF ( ABS ( A(I,I) ) .GE. ABS ( A(I+1,I) ) ) GO TO 6
        IWORK(I) = J
        DO 5 J = I, M
          R = A(I,J)
          A(I,J) = A(I+1,J)
5         A(I+1,J) = R
6       R = -A(I+1,I) / A(I,I)
        A(I+1,I) = R
        DO 7 J = L, M
7         A(I+1,J) = A(I+1,J) + R * A(I,J)
8     CONTINUE

      IF ( A(M,M) .NE. 0.0E+00 ) GO TO 9
      A(M,M) = EPS
C
C  THE VECTOR (1,1,...,1) IS STORED IN THE PLACE OF THE RIGHT
C  HAND SIDE COLUMN VECTOR.
C
9     DO 11 I = 1, N
        IF ( I .GE. M ) GO TO 10
        WORK(I) = 1.0E+00
        GO TO 11
10      WORK(I) = 0.0E+00
11    CONTINUE
C
C  THE INVERSE ITERATION IS PERFORMED ON THE MATRIX UNTIL THE
C  INFINITY NORM OF THE RIGHT-HAND SIDE VECTOR IS GREATER
C  THAN THE BOUND DEFINED AS 0.01 / ( N * EX ).
C
      BOUND = 0.01E+00 / ( EX * FLOAT ( N ) )
      NS = 0
      ITER = 1
C
C  THE BACKSUBSTITUTION.
C
12    R = 0.0E+00
      DO 15 I = 1, M
        J = M - I + 1
        S = WORK(J)
        IF ( J .EQ. M ) GO TO 14
        L = J + 1
        DO 13 K = L, M
          SR = WORK(K)
13        S = S - SR * A(J,K)
14      WORK(J) = S / A(J,J)
        T = ABS ( WORK(J) )
        IF ( R .GE. T ) GO TO 15
        R = T
15    CONTINUE
C
C  THE COMPUTATION OF THE RIGHT-HAND SIDE VECTOR FOR THE NEW
C  ITERATION STEP.
C
      DO 16 I = 1, M
16      WORK(I) = WORK(I) / R
C
C  THE COMPUTATION OF THE RESIDUALS AND COMPARISON OF THE
C  RESIDUALS OF THE TWO SUCCESSIVE STEPS OF THE INVERSE
C  ITERATION.  IF THE INFINITY NORM OF THE RESIDUAL VECTOR
C  IS GREATER THAN THE INFINITY NORM OF THE PREVIOUS RESIDUAL
C  VECTOR, THE COMPUTED EIGENVECTOR OF THE PREVIOUS STEP IS
C  TAKEN AS THE FINAL EIGENVECTOR.
C
      R1 = 0.0E+00
      DO 18 I = 1, M
        T = 0.0E+00
        DO 17 J = I, M
17        T = T + A(I,J) * WORK(J)
        T = ABS ( T )
        IF ( R1 .GE. T ) GO TO 18
        R1 = T
18    CONTINUE

      IF ( ITER .EQ. 1 ) GO TO 19
      IF ( PREVIS .LE. R1 ) GO TO 24
19    DO 20 I = 1, M
20      VECR(I,IVEC) = WORK(I)
      PREVIS = R1
      IF ( NS .EQ. 1 ) GO TO 24
      IF ( ITER .GT. 6 ) GO TO 25
      ITER = ITER + 1
      IF ( R .LT. BOUND ) GO TO 21
      NS = 1
C
C  GAUSSIAN ELIMINATION OF THE RIGHT HAND SIDE VECTOR.
C
21    K = M - 1
      DO 23 I = 1, K
        R = WORK(I+1)
        IF ( IWORK(I) .EQ. 0 ) GO TO 22
        WORK(I+1) = WORK(I) + WORK(I+1) * A(I+1,I)
        WORK(I) = R
        GO TO 23
22      WORK(I+1) = WORK(I+1) + WORK(I) * A(I+1,I)
23    CONTINUE

      GO TO 12

24    INDIC(IVEC) = 2
25    IF ( M .EQ. N ) GO TO 27
      J = M + 1
      DO 26 I = J, N
26      VECR(I,IVEC) = 0.0E+00
27    RETURN
      END
      SUBROUTINE COMPVE ( N, NM, M, IVEC, A, VECR, H, EVR, EVI, INDIC,
     &  IWORK, SUBDIA, WORK1, WORK2, WORK, EPS, EX )

c*********************************************************************72

      IMPLICIT NONE

      INTEGER NM

      REAL A(NM,1)
      REAL B
      REAL BOUND
      DOUBLE PRECISION D
      DOUBLE PRECISION D1
      REAL EPS
      REAL ETA
      REAL EVI(NM)
      REAL EVR(NM)
      REAL EX
      REAL FKSI
      REAL H(NM,1)
      INTEGER I
      INTEGER I1
      INTEGER I2
      INTEGER INDIC(NM)
      INTEGER ITER
      INTEGER IVEC
      INTEGER IWORK(NM)
      INTEGER J
      INTEGER K
      INTEGER L
      INTEGER M
      INTEGER N
      INTEGER NS
      REAL PREVIS
      REAL R
      REAL S
      REAL SUBDIA(NM)
      REAL U
      REAL V
      REAL VECR(NM,1)
      REAL WORK(NM)
      REAL WORK1(NM)
      REAL WORK2(NM)
C
C  THIS SUBROUTINE FINDS THE COMPLEX EIGENVECTOR OF THE REAL
C  UPPER-HESSENBERG MATRIX OF ORDER N CORRESPONDING TO THE
C  COMPLEX EIGENVALUE WITH THE REAL PART IN EVR(IVEC) AND THE
C  CORRESPONDING IMAGINARY PART IN EVI(IVEC).  THE INVERSE
C  ITERATION METHOD IS USED, MODIFIED TO AVOID THE USE OF
C  COMPLEX ARITHMETIC.
C
C  THE MATRIX ON WHICH THE INVERSE ITERATION IS PERFORMED IS
C  BUILT UP INTO THE ARRAY A BY USING THE UPPER-HESSENBERG
C  MATRIX PRESERVED IN THE UPPER HALF OF THE ARRAY H AND IN
C  THE ARRAY SUBDIA.
C
C  NM DEFINES THE FIRST DIMENSION OF THE TWO DIMENSIONAL
C  ARRAYS A, VECR, AND H.  NM MUST BE EQUAL TO OR GREATER
C  THAN N.
C
C  M IS THE ORDER OF THE SUBMATRIX OBTAINED BY A SUITABLE
C  DECOMPOSITION OF THE UPPER-HESSENBERG MATRIX IF SOME
C  SUBDIAGONAL ELEMENTS ARE EQUAL TO ZERO.  THE VALUE OF M IS
C  CHOSEN SO THAT THE LAST N-M COMPONENTS OF THE COMPLEX
C  EIGENVECTOR ARE ZERO.
C
C  THE REAL PARTS OF THE FIRST M COMPONENTS OF THE COMPUTED
C  COMPLEX EIGENVECTOR WILL BE FOUND IN THE FIRST M PLACES OF
C  THE COLUMN WHOSE TOP ELEMENT IS VECR(1,IVEC), AND THE
C  CORRESPONDING IMAGINARY PARTS OF THE FIRST M COMPONENTS OF
C  THE COMPLEX EIGENVECTOR WILL BE FOUND IN THE FIRST M
C  PLACES OF THE COLUMN WHOSE TOP ELEMENT IS VECR(1,IVEC-1).
C
C  THE ARRAY INDIC INDICATES THE SUCCESS OF THE ROUTINE AS
C  FOLLOWS:
C    VALUE OF INDIC(I)  EIGENVECTOR I
C           1             NOT FOUND
C           2                 FOUND
C
C  THE ARRAYS IWORK, WORK1, WORK2 AND WORK ARE THE WORKING
C  STORES USED DURING THE INVERSE ITERATION PROCESS.
C
C  EPS IS A SMALL POSITIVE NUMBER THAT NUMERICALLY REPRESENTS
C  ZERO IN THE PROGRAM.  EPS = (EUCLIDEAN NORM OF H)*EX, WHERE
C  EX = 2**(-T).  T IS THE NUMBER OF BINARY DIGITS IN THE
C  MANTISSA OF A FLOATING POINT NUMBER.
C
      FKSI = EVR(IVEC)
      ETA = EVI(IVEC)
C
C  THE MODIFICATION OF THE EIGENVALUE ( FKSI + I * ETA ) IF MORE
C  EIGENVALUES ARE EQUAL.
C
      IF ( IVEC .EQ. M ) GO TO 2
      K = IVEC + 1
      R = 0.0E+00
      DO 1 I = K, M
        IF ( FKSI .NE. EVR(I) ) GO TO 1
        IF ( ABS ( ETA ) .NE. ABS ( EVI(I) ) ) GO TO 1
        R = R + 3.0E+00
1     CONTINUE
      R = R * EX
      FKSI = FKSI + R
      ETA = ETA + R
C
C  THE MATRIX ((H-FKSI*I)*(H-FKSI*I)+(ETA*ETA)*I) IS
C  STORED INTO THE ARRAY A.
C
2     R = FKSI * FKSI + ETA * ETA
      S = 2.0E+00 * FKSI
      L = M - 1

      DO 5 I = 1, M
        DO 4 J = I, M
          D = 0.0E+00
          A(J,I) = 0.0E+00
          DO 3 K = I, J
3           D = D + H(I,K) * H(K,J)
4         A(I,J) = D - S * H(I,J)
5       A(I,I) = A(I,I) + R

      DO 9 I = 1, L
        R = SUBDIA(I)
        A(I+1,I) = -S * R
        I1 = I + 1
        DO 6 J = 1, I1
6         A(J,I) = A(J,I) + R * H(J,I+1)
        IF ( I .EQ. 1 ) GO TO 7
        A(I+1,I-1) = R * SUBDIA(I-1)
7       DO 8 J = I, M
8         A(I+1,J) = A(I+1,J) + R * H(I,J)
9     CONTINUE

C
C  THE GAUSSIAN ELIMINATION OF THE MATRIX
C  ((H-FKSI*I)*(H-FKSI*I)+(ETA*ETA)*I) IN THE ARRAY A.  THE
C  ROW INTERCHANGES THAT OCCUR ARE INDICATED IN THE ARRAY
C  IWORK.  ALL THE MULTIPLIERS ARE STORED IN THE FIRST AND IN
C  THE SECOND SUBDIAGONAL OF THE ARRAY A.
C
      K = M - 1
      DO 18 I = 1, K
        I1 = I + 1
        I2 = I + 2
        IWORK(I) = 0
        IF ( I .EQ. K ) GO TO 10
        IF ( A(I+2,I) .NE. 0.0E+00 ) GO TO 11
10      IF ( A(I+1,I) .NE. 0.0E+00 ) GO TO 11
        IF ( A(I,I) .NE. 0.0E+00 ) GO TO 18
        A(I,I) = EPS
        GO TO 18

11      IF ( I .EQ. K ) GO TO 12
        IF ( ABS ( A(I+1,I) ) .GE. ABS ( A(I+2,I) ) ) GO TO 12
        IF ( ABS ( A(I,I) ) .GE. ABS ( A(I+2,I) ) ) GO TO 16
        L = I + 2
        IWORK(I) = 2
        GO TO 13

12      IF ( ABS ( A(I,I) ) .GE. ABS ( A(I+1,I) ) ) GO TO 15
        L = I + 1
        IWORK(I) = 1

13      DO 14 J = I, M
          R = A(I,J)
          A(I,J) = A(L,J)
14        A(L,J) = R

15      IF ( I .NE. K ) GO TO 16
        I2 = I1
16      DO 17 L = I1, I2
          R = -A(L,I) / A(I,I)
          A(L,I) = R
          DO 17 J = I1, M
17          A(L,J) = A(L,J) + R * A(I,J)
18    CONTINUE

      IF ( A(M,M) .NE. 0.0E+00 ) GO TO 19
      A(M,M) = EPS
C
C  THE VECTOR (1,1,...,1) IS STORED INTO THE RIGHT-HAND SIDE
C  VECTORS VECR(*,IVEC) AND VECR(*,IVEC-1), REPRESENTING THE
C  COMPLEX RIGHT HAND SIDE VECTOR.
C
19    DO 21 I = 1, N
        IF ( I .GT. M ) GO TO 20
        VECR(I,IVEC) = 1.0E+00
        VECR(I,IVEC-1) = 1.0E+00
        GO TO 21
20      VECR(I,IVEC) = 0.0E+00
        VECR(I,IVEC-1) = 0.0E+00
21    CONTINUE
C
C  THE INVERSE ITERATION IS PERFORMED ON THE MATRIX UNTIL THE
C  INFINITY NORM OF THE RIGHT-HAND SIDE IS GREATER
C  THAN THE BOUND DEFINED AS 0.01/(N*EX).
C
      BOUND = 0.01E+00 / ( EX * FLOAT ( N ) )
      NS = 0
      ITER = 1
      DO 22 I = 1, M
22      WORK(I) = H(I,I) - FKSI
C
C  THE SEQUENCE OF THE COMPLEX VECTORS Z(S) = P(S) + I * Q(S) AND
C  W(S+1) = U(S+1)+I*V(S+1) IS GIVEN BY THE RELATIONS
C  (A-FKSI-I*ETA)*I)*W(S+1) = Z(S) AND
C  Z(S+1) = W(S+1)/MAX(W(S+1)).
C  THE FINAL W(S) IS TAKEN AS THE COMPUTED EIGENVECTOR.
C
C  THE COMPUTATION OF THE RIGHT-HAND SIDE VECTOR
C  (A-FKSI*I)*P(S)-ETA*Q(S).  A IS AN UPPER-HESSENBERG MATRIX.
C
23    DO 27 I = 1, M
        D = WORK(I) * VECR(I,IVEC)
        IF ( I .EQ. 1 ) GO TO 24
        D = D + SUBDIA(I-1) * VECR(I-1,IVEC)
24      L = I + 1
        IF ( L .GT. M ) GO TO 26
        DO 25 K = L, M
25        D = D + H(I,K) * VECR(K,IVEC)
26      VECR(I,IVEC-1) = D - ETA * VECR(I,IVEC-1)
27    CONTINUE
C
C  GAUSSIAN ELIMINATION OF THE RIGHT-HAND SIDE VECTOR.
C
      K = M - 1
      DO 28 I = 1, K
        L = I + IWORK(I)
        R = VECR(L,IVEC-1)
        VECR(L,IVEC-1) = VECR(I,IVEC-1)
        VECR(I,IVEC-1) = R
        VECR(I+1,IVEC-1) = VECR(I+1,IVEC-1) + A(I+1,I) * R
        IF ( I .EQ. K ) GO TO 28
        VECR(I+2,IVEC-1) = VECR(I+2,IVEC-1) + A(I+2,I) * R
28    CONTINUE
C
C  THE COMPUTATION OF THE REAL PART U(S+1) OF THE COMPLEX
C  VECTOR W(S+1).  THE VECTOR U(S+1) IS OBTAINED AFTER THE
C  BACKSUBSTITUTION.
C
      DO 31 I = 1, M
        J = M - I + 1
        D = VECR(J,IVEC-1)
        IF ( J .EQ. M ) GO TO 30
        L = J + 1
        DO 29 K = L, M
          D1 = A(J,K)
29        D = D - D1 * VECR(K,IVEC-1)
30      VECR(J,IVEC-1) = D / A(J,J)
31    CONTINUE
C
C  THE COMPUTATION OF THE IMAGINARY PART V(S+1) OF THE VECTOR
C  W(S+1), WHERE V(S+1) = (P(S)-(A-FKSI*I)*U(S+1))/ETA.
C
      DO 35 I = 1, M
        D = WORK(I) * VECR(I,IVEC-1)
        IF ( I .EQ. 1 ) GO TO 32
        D = D + SUBDIA(I-1) * VECR(I-1,IVEC-1)
32      L = I + 1
        IF ( L .GT. M ) GO TO 34
        DO 33 K = L, M
33        D = D + H(I,K) * VECR(K,IVEC-1)
34      VECR(I,IVEC) = ( VECR(I,IVEC) - D ) / ETA
35    CONTINUE
C
C  THE COMPUTATION OF (INFINITY NORM OF W(S+1))**2.
C
      L = 1
      S = 0.0E+00
      DO 36 I = 1, M
        R = VECR(I,IVEC)**2 + VECR(I,IVEC-1)**2
        IF ( R .LE. S ) GO TO 36
        S = R
        L = I
36    CONTINUE
C
C  THE COMPUTATION OF THE VECTOR Z(S+1), WHERE Z(S+1) = W(S+1)/
C  (COMPONENT OF W(S+1) WITH THE LARGEST ABSOLUTE VALUE).
C
      U = VECR(L,IVEC-1)
      V = VECR(L,IVEC)
      DO 37 I = 1, M
        B = VECR(I,IVEC)
        R = VECR(I,IVEC-1)
        VECR(I,IVEC) = ( R * U + B * V ) / S
37      VECR(I,IVEC-1) = ( B * U - R * V ) / S
C
C  THE COMPUTATION OF THE RESIDUALS AND COMPARISON OF THE
C  RESIDUALS OF THE TWO SUCCESSIVE STEPS OF THE INVERSE
C  ITERATION.  IF THE INFINITY NORM OF THE RESIDUAL VECTOR IS
C  GREATER THAN THE INFINITY NORM OF THE PREVIOUS RESIDUAL
C  VECTOR, THE COMPUTED VECTOR OF THE PREVIOUS STEP IS TAKEN
C  AS THE COMPUTED APPROXIMATION TO THE EIGENVECTOR.
C
      B = 0.0E+00
      DO 41 I = 1, M
        R = WORK(I) * VECR(I,IVEC-1) - ETA * VECR(I,IVEC)
        U = WORK(I) * VECR(I,IVEC) + ETA * VECR(I,IVEC-1)
        IF ( I .EQ. 1 ) GO TO 38
        R = R + SUBDIA(I-1) * VECR(I-1,IVEC-1)
        U = U + SUBDIA(I-1) * VECR(I-1,IVEC)
38      L = I + 1
        IF ( L .GT. M ) GO TO 40
        DO 39 J = L, M
          R = R + H(I,J) * VECR(J,IVEC-1)
39        U = U + H(I,J) * VECR(J,IVEC)
40      U = R * R + U * U
        IF ( B .GE. U ) GO TO 41
        B = U
41    CONTINUE

      IF ( ITER .EQ. 1 ) GO TO 42
      IF ( PREVIS .LE. B ) GO TO 44
42    DO 43 I = 1, N
        WORK1(I) = VECR(I,IVEC)
43      WORK2(I) = VECR(I,IVEC-1)

      PREVIS = B
      IF ( NS .EQ. 1 ) GO TO 46
      IF ( ITER .GT. 6 ) GO TO 47
      ITER = ITER + 1
      IF ( BOUND .GT. SQRT ( S ) ) GO TO 23
      NS = 1
      GO TO 23

44    DO 45 I = 1, N
        VECR(I,IVEC) = WORK1(I)
45      VECR(I,IVEC-1) = WORK2(I)
46    INDIC(IVEC-1) = 2
      INDIC(IVEC) = 2

47    RETURN
      END
