c    @(#)hqr.f	1.1
c
      subroutine hqr(nm,n,low,igh,h,wr,wi,ierr)
      integer i,j,k,l,m,n,en,ll,mm,na,nm,igh,its,low,mp2,enm2,ierr
      real h(nm,n),wr(n),wi(n)
      real p,q,r,s,t,w,x,y,zz,norm,machep
      logical notlas
      machep = 2.**(-47)
      ierr = 0
      norm = 0.0
      k = 1
      do 50 i = 1, n
         do 40 j = k, n
   40    norm = norm + abs(h(i,j))
         k = i
         if (i .ge. low .and. i .le. igh) go to 50
         wr(i) = h(i,i)
         wi(i) = 0.0
   50 continue
      en = igh
      t = 0.0
   60 if (en .lt. low) go to 1001
      its = 0
      na = en - 1
      enm2 = na - 1
   70 do 80 ll = low, en
         l = en + low - ll
         if (l .eq. low) go to 100
         s = abs(h(l-1,l-1)) + abs(h(l,l))
         if (s .eq. 0.0) s = norm
         if (abs(h(l,l-1)) .le. machep * s) go to 100
   80 continue
  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 (its .eq. 30) go to 1000
      if (its .ne. 10 .and. its .ne. 20) go to 130
      t = t + x
      do 120 i = low, en
  120 h(i,i) = h(i,i) - x
      s = abs(h(en,na)) + abs(h(na,enm2))
      x = 0.75 * s
      y = x
      w = -0.4375 * s * s
  130 its = its + 1
      do 140 mm = l, enm2
         m = enm2 + l - mm
         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
         if (abs(h(m,m-1)) * (abs(q) + abs(r)) .le. machep * abs(p)
     1    * (abs(h(m-1,m-1)) + abs(zz) + abs(h(m+1,m+1)))) go to 150
  140 continue
  150 mp2 = m + 2
      do 160 i = mp2, en
         h(i,i-2) = 0.0
         if (i .eq. mp2) go to 160
         h(i,i-3) = 0.0
  160 continue
      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.0
         if (notlas) r = h(k+2,k-1)
         x = abs(p) + abs(q) + abs(r)
         if (x .eq. 0.0) 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
         do 210 j = k, en
            p = h(k,j) + q * h(k+1,j)
            if (.not. notlas) go to 200
            p = p + r * h(k+2,j)
            h(k+2,j) = h(k+2,j) - p * zz
  200       h(k+1,j) = h(k+1,j) - p * y
            h(k,j) = h(k,j) - p * x
  210    continue
         j = min0(en,k+3)
         do 230 i = l, j
            p = x * h(i,k) + y * h(i,k+1)
            if (.not. notlas) go to 220
            p = p + zz * h(i,k+2)
            h(i,k+2) = h(i,k+2) - p * r
  220       h(i,k+1) = h(i,k+1) - p * q
            h(i,k) = h(i,k) - p
  230    continue
  260 continue
      go to 70
  270 wr(en) = x + t
      wi(en) = 0.0
      en = na
      go to 60
  280 p = (y - x) / 2.0
      q = p * p + w
      zz = sqrt(abs(q))
      x = x + t
      if (q .lt. 0.0) go to 320
      zz = p + sign(zz,p)
      wr(na) = x + zz
      wr(en) = wr(na)
      if (zz .ne. 0.0) wr(en) = x - w / zz
      wi(na) = 0.0
      wi(en) = 0.0
      go to 330
  320 wr(na) = x + p
      wr(en) = x + p
      wi(na) = zz
      wi(en) = -zz
  330 en = enm2
      go to 60
 1000 ierr = en
 1001 return
      end
