subroutine ratqr(n,eps1,d,e,e2,m,w,ind,bd,type,idef,ierr)
c
      integer i,j,k,m,n,ii,jj,k1,idef,ierr,jdef
      double precision d(n),e(n),e2(n),w(n),bd(n)
      double precision f,p,q,r,s,ep,qp,err,tot,eps1,delta,epslon
      integer ind(n)
      logical type
c
c     this subroutine is a translation of the algol procedure ratqr,
c     num. math. 11, 264-272(1968) by reinsch and bauer.
c     handbook for auto. comp., vol.ii-linear algebra, 257-265(1971).
c
c     this subroutine finds the algebraically smallest or largest
c     eigenvalues of a symmetric tridiagonal matrix by the
c     rational qr method with newton corrections.
c
c     on input
c
c        n is the order of the matrix.
c
c        eps1 is a theoretical absolute error tolerance for the
c          computed eigenvalues.  if the input eps1 is non-positive,
c          or indeed smaller than its default value, it is reset
c          at each iteration to the respective default value,
c          namely, the product of the relative machine precision
c          and the magnitude of the current eigenvalue iterate.
c          the theoretical absolute error in the k-th eigenvalue
c          is usually not greater than k times eps1.
c
c        d contains the diagonal elements of the input matrix.
c
c        e contains the subdiagonal elements of the input matrix
c          in its last n-1 positions.  e(1) is arbitrary.
c
c        e2 contains the squares of the corresponding elements of e.
c          e2(1) is arbitrary.
c
c        m is the number of eigenvalues to be found.
c
c        idef should be set to 1 if the input matrix is known to be
c          positive definite, to -1 if the input matrix is known to
c          be negative definite, and to 0 otherwise.
c
c        type should be set to .true. if the smallest eigenvalues
c          are to be found, and to .false. if the largest eigenvalues
c          are to be found.
c
c     on output
c
c        eps1 is unaltered unless it has been reset to its
c          (last) default value.
c
c        d and e are unaltered (unless w overwrites d).
c
c        elements of e2, corresponding to elements of e regarded
c          as negligible, have been replaced by zero causing the
c          matrix to split into a direct sum of submatrices.
c          e2(1) is set to 0.0d0 if the smallest eigenvalues have been
c          found, and to 2.0d0 if the largest eigenvalues have been
c          found.  e2 is otherwise unaltered (unless overwritten by bd).
c
c        w contains the m algebraically smallest eigenvalues in
c          ascending order, or the m largest eigenvalues in
c          descending order.  if an error exit is made because of
c          an incorrect specification of idef, no eigenvalues
c          are found.  if the newton iterates for a particular
c          eigenvalue are not monotone, the best estimate obtained
c          is returned and ierr is set.  w may coincide with d.
c
c        ind contains in its first m positions the submatrix indices
c          associated with the corresponding eigenvalues in w --
c          1 for eigenvalues belonging to the first submatrix from
c          the top, 2 for those belonging to the second submatrix, etc..
c
c        bd contains refined bounds for the theoretical errors of the
c          corresponding eigenvalues in w.  these bounds are usually
c          within the tolerance specified by eps1.  bd may coincide
c          with e2.
c
c        ierr is set to
c          zero       for normal return,
c          6*n+1      if  idef  is set to 1 and  type  to .true.
c                     when the matrix is not positive definite, or
c                     if  idef  is set to -1 and  type  to .false.
c                     when the matrix is not negative definite,
c          5*n+k      if successive iterates to the k-th eigenvalue
c                     are not monotone increasing, where k refers
c                     to the last such occurrence.
c
c     note that subroutine tridib is generally faster and more
c     accurate than ratqr if the eigenvalues are clustered.
c
c     questions and comments should be directed to burton s. garbow,
c     mathematics and computer science div, argonne national laboratory
c
c     this version dated august 1983.
c
c     ------------------------------------------------------------------
c
      ierr = 0
      jdef = idef
c     .......... copy d array into w ..........
      do 20 i = 1, n
   20 w(i) = d(i)
c
      if (type) go to 40
      j = 1
      go to 400
   40 err = 0.0d0
      s = 0.0d0
c     .......... look for small sub-diagonal entries and define
c                initial shift from lower gerschgorin bound.
c                copy e2 array into bd ..........
      tot = w(1)
      q = 0.0d0
      j = 0
c
      do 100 i = 1, n
         p = q
         if (i .eq. 1) go to 60
         if (p .gt. epslon(dabs(d(i)) + dabs(d(i-1)))) go to 80
   60    e2(i) = 0.0d0
   80    bd(i) = e2(i)
c     .......... count also if element of e2 has underflowed ..........
         if (e2(i) .eq. 0.0d0) j = j + 1
         ind(i) = j
         q = 0.0d0
         if (i .ne. n) q = dabs(e(i+1))
         tot = dmin1(w(i)-p-q,tot)
  100 continue
c
      if (jdef .eq. 1 .and. tot .lt. 0.0d0) go to 140
c
      do 110 i = 1, n
  110 w(i) = w(i) - tot
c
      go to 160
  140 tot = 0.0d0
c
  160 do 360 k = 1, m
c     .......... next qr transformation ..........
  180    tot = tot + s
         delta = w(n) - s
         i = n
         f = dabs(epslon(tot))
         if (eps1 .lt. f) eps1 = f
         if (delta .gt. eps1) go to 190
         if (delta .lt. (-eps1)) go to 1000
         go to 300
c     .......... replace small sub-diagonal squares by zero
c                to reduce the incidence of underflows ..........
  190    if (k .eq. n) go to 210
         k1 = k + 1
         do 200 j = k1, n
            if (bd(j) .le. (epslon(w(j)+w(j-1))) ** 2) bd(j) = 0.0d0
  200    continue
c
  210    f = bd(n) / delta
         qp = delta + f
         p = 1.0d0
         if (k .eq. n) go to 260
         k1 = n - k
c     .......... for i=n-1 step -1 until k do -- ..........
         do 240 ii = 1, k1
            i = n - ii
            q = w(i) - s - f
            r = q / qp
            p = p * r + 1.0d0
            ep = f * r
            w(i+1) = qp + ep
            delta = q - ep
            if (delta .gt. eps1) go to 220
            if (delta .lt. (-eps1)) go to 1000
            go to 300
  220       f = bd(i) / q
            qp = delta + f
            bd(i+1) = qp * ep
  240    continue
c
  260    w(k) = qp
         s = qp / p
         if (tot + s .gt. tot) go to 180
c     .......... set error -- irregular end of iteration.
c                deflate minimum diagonal element ..........
         ierr = 5 * n + k
         s = 0.0d0
         delta = qp
c
         do 280 j = k, n
            if (w(j) .gt. delta) go to 280
            i = j
            delta = w(j)
  280    continue
c     .......... convergence ..........
  300    if (i .lt. n) bd(i+1) = bd(i) * f / qp
         ii = ind(i)
         if (i .eq. k) go to 340
         k1 = i - k
c     .......... for j=i-1 step -1 until k do -- ..........
         do 320 jj = 1, k1
            j = i - jj
            w(j+1) = w(j) - s
            bd(j+1) = bd(j)
            ind(j+1) = ind(j)
  320    continue
c
  340    w(k) = tot
         err = err + dabs(delta)
         bd(k) = err
         ind(k) = ii
  360 continue
c
      if (type) go to 1001
      f = bd(1)
      e2(1) = 2.0d0
      bd(1) = f
      j = 2
c     .......... negate elements of w for largest values ..........
  400 do 500 i = 1, n
  500 w(i) = -w(i)
c
      jdef = -jdef
      go to (40,1001), j
c     .......... set error -- idef specified incorrectly ..........
 1000 ierr = 6 * n + 1
 1001 return
      end