@process directive('"    (')
      subroutine corth(nm,n,low,igh,ar,ai,ortr,orti)
c
      integer i,j,m,n,mp,nm,igh,low
      real ar(nm,n),ai(nm,n),ortr(igh),orti(igh)
      real f,g,h,fi,fr,scale
c
c     this subroutine is a translation of a complex analogue of
c     the algol procedure orthes, num. math. 12, 349-368(1968)
c     by martin and wilkinson.
c     handbook for auto. comp., vol.ii-linear algebra, 339-358(1971).
c
c     given a complex general matrix, this subroutine
c     reduces a submatrix situated in rows and columns
c     low through igh to upper hessenberg form by
c     unitary similarity transformations.
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  cbal.  if  cbal  has not been used,
c          set low=1, igh=n.
c
c        ar and ai contain the real and imaginary parts,
c          respectively, of the complex input matrix.
c
c     on output
c
c        ar and ai contain the real and imaginary parts,
c          respectively, of the hessenberg matrix.  information
c          about the unitary transformations used in the reduction
c          is stored in the remaining triangles under the
c          hessenberg matrix.
c
c        ortr and orti contain further information about the
c          transformations.  only elements low through igh are used.
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)
      do 180 m = low+1, igh-1
         h = 0.0e0
         ortr(m) = 0.0e0
         orti(m) = 0.0e0
         scale = 0.0e0
c     .......... scale column (algol tol then not needed) ..........
         do 90 i = m, igh
   90    scale = scale + abs(ar(i,m-1)) + abs(ai(i,m-1))
c
         if (scale .eq. 0.0e0) go to 180
         mp = m + igh
c     .......... for i=igh step -1 until m do -- ..........
         do 100 i = m, igh
            ortr(i) = ar(i,m-1) / scale
            orti(i) = ai(i,m-1) / scale
            h = h + ortr(i) * ortr(i) + orti(i) * orti(i)
  100    continue
c
         g = sqrt(h)
c--         f = pythag(ortr(m),orti(m))
         if (abs(ortr(m)).gt.abs(orti(m))) then
            f = abs(ortr(m))*sqrt(1e0+(orti(m)/ortr(m))**2)
         else if (orti(m).ne.0e0) then
            f = abs(orti(m))*sqrt((ortr(m)/orti(m))**2+1e0)
         else
            f = abs(ortr(m))
         endif
         if (f .eq. 0.0e0) go to 103
         h = h + f * g
         g = g / f
         ortr(m) = (1.0e0 + g) * ortr(m)
         orti(m) = (1.0e0 + g) * orti(m)
         go to 105
c
  103    ortr(m) = g
         ar(m,m-1) = scale
c     .......... form (i-(u*ut)/h) * a ..........
c"    ( prefer vector
  105    do 130 j = m, n
            fr = 0.0e0
            fi = 0.0e0
c     .......... for i=igh step -1 until m do -- ..........
            do 110 i = m, igh
               fr = fr + ortr(i) * ar(i,j) + orti(i) * ai(i,j)
               fi = fi + ortr(i) * ai(i,j) - orti(i) * ar(i,j)
  110       continue
c
            fr = -(fr / h)
            fi = -(fi / h)
c
            do 120 i = m, igh
               ar(i,j) = ar(i,j) + fr * ortr(i) - fi * orti(i)
               ai(i,j) = ai(i,j) + fr * orti(i) + fi * ortr(i)
  120       continue
c
  130    continue
c     .......... form (i-(u*ut)/h)*a*(i-(u*ut)/h) ..........
         do 160 i = 1, igh
            fr = 0.0e0
            fi = 0.0e0
c     .......... for j=igh step -1 until m do -- ..........
            do 140 j = m, igh
               fr = fr + ortr(j) * ar(i,j) - orti(j) * ai(i,j)
               fi = fi + ortr(j) * ai(i,j) + orti(j) * ar(i,j)
  140       continue
c
            fr = -(fr / h)
            fi = -(fi / h)
c
            do 150 j = m, igh
               ar(i,j) = ar(i,j) + fr * ortr(j) + fi * orti(j)
               ai(i,j) = ai(i,j) - fr * orti(j) + fi * ortr(j)
  150       continue
c
  160    continue
c
         ortr(m) = scale * ortr(m)
         orti(m) = scale * orti(m)
         ar(m,m-1) = -g * ar(m,m-1)
         ai(m,m-1) = -g * ai(m,m-1)
  180 continue
c
      return
      end