@process directive('"    (')
      subroutine hqr(nm,n,low,igh,h,wr,wi,ierr)
c
      integer i,j,k,l,m,n,en,na,nm,igh,itn,its,low,enm2,ierr
      real h(nm,n),wr(n),wi(n)
      real p,q,r,s,t,w,x,y,zz,norm,tst1,tst2,foo
      logical notlas
c
c     this subroutine is a translation of the algol procedure hqr,
c     num. math. 14, 219-231(1970) by martin, peters, and wilkinson.
c     handbook for auto. comp., vol.ii-linear algebra, 359-371(1971).
c
c     this subroutine finds the eigenvalues of a real
c     upper hessenberg matrix by the qr method.
c
c     on input
c
c        nm must be set to the row dimension of two-dimensional
c          array parameters as declared in the calling program
c          dimension statement.
c
c        n is the order of the matrix.
c
c        low and igh are integers determined by the balancing
c          subroutine  balanc.  if  balanc  has not been used,
c          set low=1, igh=n.
c
c        h contains the upper hessenberg matrix.  information about
c          the transformations used in the reduction to hessenberg
c          form by  elmhes  or  orthes, if performed, is stored
c          in the remaining triangle under the hessenberg matrix.
c
c     on output
c
c        h has been destroyed.  therefore, it must be saved
c          before calling  hqr  if subsequent calculation and
c          back transformation of eigenvectors is to be performed.
c
c        wr and wi contain the real and imaginary parts,
c          respectively, of the eigenvalues.  the eigenvalues
c          are unordered except that complex conjugate pairs
c          of values appear consecutively with the eigenvalue
c          having the positive imaginary part first.  if an
c          error exit is made, the eigenvalues should be correct
c          for indices ierr+1,...,n.
c
c        ierr is set to
c          zero       for normal return,
c          j          if the limit of 30*n iterations is exhausted
c                     while the j-th eigenvalue is being sought.
c
c     Questions and comments should be directed to Alan K. Cline,
c     Pleasant Valley Software, 8603 Altus Cove, Austin, TX 78759.
c     Electronic mail to cline@cs.utexas.edu.
c
c     this version dated january 1989. (for the IBM 3090vf)
c
c     ------------------------------------------------------------------
c
      call xuflow(0)
      ierr = 0
      norm = 0.0e0
c     .......... store roots isolated by balanc
c                and compute matrix norm ..........
      do 50 j = 1, n
         do 50 i = 1,min0(n,j+1)
   50       norm = norm + abs(h(i,j))
c
c"    ( prefer vector
      do 52 i = 1, low-1
         wr(i) = h(i,i)
         wi(i) = 0.0e0
   52 continue
c
c"    ( prefer vector
      do 54 i = igh+1,n
         wr(i) = h(i,i)
         wi(i) = 0.0e0
   54 continue
c
      en = igh
      t = 0.0e0
      itn = 30*n
c     .......... search for next eigenvalues ..........
   60 if (en .lt. low) go to 1001
      its = 0
      na = en - 1
      enm2 = na - 1
c     .......... look for single small sub-diagonal element
c                for l=en step -1 until low do -- ..........
c     .......... changed lower bound to "low+1" and removed
c                test "if (l.eq.low) goto 100" [JMM] -- ..........
   70 do 80 l = en, low+1, -1
         s = abs(h(l-1,l-1)) + abs(h(l,l))
         if (s .eq. 0.0e0) s = norm
         tst1 = s
         tst2 = tst1 + abs(h(l,l-1))
         if (tst2 .eq. tst1) go to 100
   80 continue
c     .......... form shift ..........
  100 x = h(en,en)
      if (l .eq. en) go to 270
      y = h(na,na)
      w = h(en,na) * h(na,en)
      if (l .eq. na) go to 280
      if (itn .eq. 0) go to 1000
      if (its .ne. 10 .and. its .ne. 20) go to 130
c     .......... form exceptional shift ..........
      t = t + x
c
c"    ( prefer vector
      do 120 i = low, en
  120 h(i,i) = h(i,i) - x
c
      s = abs(h(en,na)) + abs(h(na,enm2))
      x = 0.75e0 * s
      y = x
      w = -0.4375e0 * s * s
  130 its = its + 1
      itn = itn - 1
c     .......... look for two consecutive small
c                sub-diagonal elements.
c                for m=en-2 step -1 until l do -- ..........
      do 140 m = en-2, l, -1
         zz = h(m,m)
         r = x - zz
         s = y - zz
         p = (r * s - w) / h(m+1,m) + h(m,m+1)
         q = h(m+1,m+1) - zz - r - s
         r = h(m+2,m+1)
         s = abs(p) + abs(q) + abs(r)
         p = p / s
         q = q / s
         r = r / s
         if (m .eq. l) go to 150
         tst1 = abs(p)*(abs(h(m-1,m-1)) + abs(zz) + abs(h(m+1,m+1)))
         tst2 = tst1 + abs(h(m,m-1))*(abs(q) + abs(r))
         if (tst2 .eq. tst1) go to 150
  140 continue
c
  150 h(m+2,m) = 0.0e0
c"    ( prefer vector
      do 160 i = m+3, en
         h(i,i-2) = 0.0e0
         h(i,i-3) = 0.0e0
  160 continue
c     .......... double qr step involving rows l to en and
c                columns m to en ..........
      do 260 k = m, na
         notlas = k .ne. na
         if (k .eq. m) go to 170
         p = h(k,k-1)
         q = h(k+1,k-1)
         r = 0.0e0
         if (notlas) r = h(k+2,k-1)
         x = abs(p) + abs(q) + abs(r)
         if (x .eq. 0.0e0) go to 260
         p = p / x
         q = q / x
         r = r / x
  170    s = sign(sqrt(p*p+q*q+r*r),p)
         if (k .eq. m) go to 180
         h(k,k-1) = -s * x
         go to 190
  180    if (l .ne. m) h(k,k-1) = -h(k,k-1)
  190    p = p + s
         x = p / s
         y = q / s
         zz = r / s
         q = q / p
         r = r / p
         if (notlas) go to 225
c     .......... row modification ..........
c"    ( prefer vector
         do 200 j = k, en
            foo = h(k,j) + q * h(k+1,j)
            h(k,j) = h(k,j) - foo * x
            h(k+1,j) = h(k+1,j) - foo * y
  200    continue
c
         j = min0(en,k+3)
c     .......... column modification ..........
         do 210 i = l, j
            foo = x * h(i,k) + y * h(i,k+1)
            h(i,k) = h(i,k) - foo
            h(i,k+1) = h(i,k+1) - foo * q
  210    continue
         go to 255
  225    continue
c     .......... row modification ..........
c"    ( prefer vector
         do 230 j = k, en
            foo = h(k,j) + q * h(k+1,j) + r * h(k+2,j)
            h(k,j) = h(k,j) - foo * x
            h(k+1,j) = h(k+1,j) - foo * y
            h(k+2,j) = h(k+2,j) - foo * zz
  230    continue
c
         j = min0(en,k+3)
c     .......... column modification ..........
         do 240 i = l, j
            foo = x * h(i,k) + y * h(i,k+1) + zz * h(i,k+2)
            h(i,k) = h(i,k) - foo
            h(i,k+1) = h(i,k+1) - foo * q
            h(i,k+2) = h(i,k+2) - foo * r
  240    continue
  255    continue
c
  260 continue
c
      go to 70
c     .......... one root found ..........
  270 wr(en) = x + t
      wi(en) = 0.0e0
      en = na
      go to 60
c     .......... two roots found ..........
  280 p = (y - x) / 2.0e0
      q = p * p + w
      zz = sqrt(abs(q))
      x = x + t
      if (q .lt. 0.0e0) go to 320
c     .......... real pair ..........
      zz = p + sign(zz,p)
      wr(na) = x + zz
      wr(en) = wr(na)
      if (zz .ne. 0.0e0) wr(en) = x - w / zz
      wi(na) = 0.0e0
      wi(en) = 0.0e0
      go to 330
c     .......... complex pair ..........
  320 wr(na) = x + p
      wr(en) = x + p
      wi(na) = zz
      wi(en) = -zz
  330 en = enm2
      go to 60
c     .......... set error -- all eigenvalues have not
c                converged after 30*n iterations ..........
 1000 ierr = en
 1001 return
      end