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 k1 = 1, m
          do l1 = k1, m
            a(k1,l1) = veci(k1,l1)
          end do
          if ( k1 .ne. 1 ) then
            a(k1,k1-1) = subdia(k1-1)
          end if
        end do
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 i = 1, n
        do j = i, n
          a(i,j) = 0.0e+00
          a(j,i) = 0.0e+00
        end do
        a(i,i) = 1.0e+00
      end do

      do k = 1, n - 2
        l = k + 1
        do j = 2, n
          d1 = 0.0e+00
          do i = l, n
            d2 = veci(i,k)
            d1 = d1 + d2 * a(j,i)
          end do
          do i = l, n
            a(j,i) = a(j,i) - real ( veci(i,k) * d1 )
          end do
        end do
      end do
c
c  the computation of the eigenvectors of the original non-
c  scaled matrix.
c
      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      continue

        do j = 1, n
          d1 = 0.0e+00
          d2 = 0.0e+00
          do k = 1, n
            d3 = a(j,k)
            d1 = d1 + d3 * vecr(k,i)
            if ( l .ne. 0 ) then
              d2 = d2 + d3 * vecr(k,i-1)
            end if
          end do
          work(j) = real ( d1 / prfact(j) )
          if ( l .ne. 0 ) then
            subdia(j) = real ( d2 / prfact(j) )
          end if
        end do
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 m = 1, n
          d1 = d1 + work(m)**2
        end do
        d1 = sqrt ( d1 )

        do m = 1, n
          veci(m,i) = 0.0e+00
          vecr(m,i) = real ( work(m) / d1 )
        end do

        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 .lt. r1 ) then
            r = r1
            l = j
          end if
22      continue

        d3 = work(l)
        r1 = subdia(l)
        do 23 j = 1, n
          d1 = work(j)
          d2 = subdia(j)
          vecr(j,i) = real ( ( d1 * d3 + d2 * r1 ) / r )
          veci(j,i) = real ( ( 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,*)
      real bound1
      real bound2
      double precision column
      real enorm
      double precision factor
      double precision fnorm
      real h(nm,*)
      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 i = 1, n
        do j = 1, n
          h(i,j) = a(i,j)
        end do
        prfact(i) = 1.0e+00
      end do

      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 j = 1, n
          if ( i .ne. j ) then
            a(i,j) = real ( a(i,j) * factor )
            a(j,i) = real ( a(j,i) / factor )
          end if
        end do

        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 i = 1, n
        do j = 1, n
          q = a(i,j)
          fnorm = fnorm + q * q
        end do
      end do
      fnorm = dsqrt ( fnorm )

      do i = 1, n
        do j = 1, n
          a(i,j) = a(i,j) / fnorm
        end do
      end do

      enorm = real ( fnorm )
      go to 13

11    do i = 1, n
        do j = 1, n
          a(i,j) = h(i,j)
        end do
      end do

      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,*)
      real eps
      real ex
      real evi(nm)
      real evr(nm)
      real h(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

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 i = l, n
          h(i,k) = a(i,k)
          s = s + abs ( a(i,k) )
        end do

        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 i = l, n
          sr = a(i,k)
          sr = sr / s
          a(i,k) = real ( sr )
          sr2 = sr2 + sr * sr
        end do

        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 i = l, n
          h(i,k) = h(i,k) / x
          subdia(i) = a(i,k) / sr2
        end do
c
c  premultiplication by the matrix pr.
c
        do j = l, n
          sr = 0.0e+00
          do i = l, n
            sr = sr + a(i,k) * a(i,j)
          end do
          do i = l, n
            a(i,j) = a(i,j) - subdia(i) * sr
          end do
        end do
c
c  postmultiplication by the matrix pr.
c
        do j = 1, n
          sr = 0.0e+00
          do i = l, n
            sr = sr + a(j,i) * a(i,k)
          end do
          do i = l, n
            a(j,i) = a(j,i) - subdia(i) * sr
          end do
        end do

12    continue

      do k = 1, m
        a(k+1,k) = subdia(k)
      end do
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 k = 1, n
        indic(k) = 0
        if ( k .ne. n ) eps = eps + subdia(k)**2
        do i = k, n
          h(k,i) = a(k,i)
          eps = eps + a(k,i)**2
        end do
      end do

      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 i = 1, n
        indic(i) = 1
        evr(i) = a(i,i)
        evi(i) = 0.0e+00
      end do

      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
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) = real ( a(i+1,j) - s * x )
          if ( i + 2 .gt. m ) go to 28
          a(i+2,j) = real ( 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) = real ( a(j,i) - s )
          a(j,i+1) = real ( a(j,i+1) - s * x )
          if ( i + 2 .le. m ) then
            a(j,i+2) = real ( a(j,i+2) - s * y )
          end if
31      continue

        if ( i + 3 .gt. m ) go to 33
        s = -a(i+3,i+2) * y * sr
32      a(i+3,i) = real ( s )
        a(i+3,i+1) = real ( s * x )
        a(i+3,i+2) = real ( 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 = real ( 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    continue

      t = real ( dsqrt ( -s ) )
      evr(k) = r
      evi(k) = t
      evr(m) = r
      evi(m) = -t
      m = m - 2
      go to 18

37    continue

      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 k = l, m
          sr = work(k)
          s = s - sr * a(j,k)
        end do
14      work(j) = real ( 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 i = 1, m
        work(i) = work(i) / r
      end do
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,*)
      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,*)
      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,*)
      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 k = i, j
            d = d + h(i,k) * h(k,j)
          end do
4         a(i,j) = real ( 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 j = 1, i1
          a(j,i) = a(j,i) + r * h(j,i+1)
        end do
        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
        do k = l, m
          d = d + h(i,k) * vecr(k,ivec)
        end do
        vecr(i,ivec-1) = real ( 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 k = l, m
          d1 = a(j,k)
          d = d - d1 * vecr(k,ivec-1)
        end do
30      vecr(j,ivec-1) = real ( 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
        do k = l, m
          d = d + h(i,k) * vecr(k,ivec-1)
        end do
        vecr(i,ivec) = real ( ( 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 i = 1, m
        r = vecr(i,ivec)**2 + vecr(i,ivec-1)**2
        if ( r .gt. s ) then
          s = r
          l = i
        end if
      end do
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 i = 1, n
        work1(i) = vecr(i,ivec)
        work2(i) = vecr(i,ivec-1)
      end do

      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