subroutine symqr ( a, d, e, k0, n, na, eps, abscnv, vec, trd, 
     &  fail )

c*********************************************************************72
c
c  explanation of the parameters in the calling sequence:
c
c    a      a double dimensioned array.  if the matrix is not
c           initially tridiagonal, it is contained in the lower
c           triangle of a.  if eigenvectors are not requested,
c           the lower triangle of a is destroyed while the
c           elements above the diagonal are left undisturbed.
c           if eigenvectors are requested, they are returned in the
c           columns of a.
c
c    d      a singly subscripted array.  if the matrix is
c           initially tridiagonal, d contains its diagonal 
c           elements.  on return, d contains the eigenvalues of
c           the matrix.
c
c    e      a singly subscripted array.  if the matrix is
c           initially tridiagonal, e contains its off-diagonal
c           elements.  upon return, e(i) contains the number of
c           iterations required to compute the approximate
c           eigenvalue d(i).
c
c    k0     a real variable containing an initial origin shift to
c           be used until the computed shifts settle down.
c
c    n      an integer variable containing the first dimension
c           of the array a.
c
c    eps    a real variable containing a convergence tolerance.
c
c    abscnv a logical variable containing the value .true. if
c           the absolute convergence criterion is to be used
c           or the value .false. if the relative criterion
c           is to be used.
c
c    vec    a logical variable containing the value .true. if
c           eigenvectors are to be computed and returned in
c           the array a and otherwise containing the value
c           .false..
c
c    trd    a logical variable containing the value .true.
c           if the matrix is tridiagonal and located in the arrays
c           d and e and otherwise containing the value .false..
c
c    fail   an integer variable containing an error signal.
c           on return, the eigenvalues in d(fail+1), ..., d(n)
c           and their corresponding eigenvectors may be presumed
c           accurate.
c
      implicit none

      integer n
      integer na

      real a(na,n)
      logical abscnv
      real c
      real cb
      real cc
      real cd
      real con
      real d(n)
      real e(n)
      real eps
      integer fail
      integer i
      integer j
      real k
      real k0
      real k1
      real k2
      integer l
      integer l1
      integer ll
      integer ll1
      real max
      real ninf
      integer nl
      integer nm1
      integer nm2
      integer nnl
      real norm
      integer nu
      integer num1
      real p
      real q
      real r
      real s
      real s2
      logical shft
      real sum
      real sum1
      real temp
      real test
      real titter
      logical trd
      logical vec

      titter = 50.0e+00
      nm1 = n - 1
      nm2 = n - 2
      ninf = 0.0e+00
      fail = 0
c
c  signal error if n is not positive.
c
      if ( n .gt. 0 ) go to 1
      fail = -1
      return
c
c  special treatment for a matrix of order one.
c
1     if ( n .gt. 1 ) go to 5
      if ( .not. trd ) d(1) = a(1,1)
      if ( vec ) a(1,1) = 1.0e+00
      fail = 0
      return
c
c  if the matrix is tridiagonal, skip the reduction.c
c
5     if ( trd ) go to 100
      if ( n .eq. 2 ) go to 80
c
c  reduce the matrix to tridiagonal form by householder's method.
c
      do 70 l = 1, nm2

        l1 = l + 1
        d(l) = a(l,l)
        max = 0.0e+00

        do i = l1, n
          max = amax1 ( max, abs ( a(i,l) ) )
        end do

        if ( max .ne. 0.0 ) go to 13
        e(l) = 0.0e+00
        a(l,l) = 1.0e+00
        go to 70
13      sum = 0.0e+00

        do i = l1, n
          a(i,l) = a(i,l) / max
          sum = sum + a(i,l)**2
        end do

        s2 = sum
        s2 = sqrt ( s2 )
        if ( a(l1,l) .lt. 0.0 ) s2 = -s2
        e(l) = -s2 * max
        a(l1,l) = a(l1,l) + s2
        a(l,l) = s2 * a(l1,l)

        sum1 = 0.0e+00
        do i = l1, n
          sum = 0.0e+00
          do j = l1, i
            sum = sum + a(i,j) * a(j,l)
          end do

          do j = i + 1, n
            sum = sum + a(j,l) * a(j,i)
          end do
          e(i) = sum / a(l,l)
          sum1 = sum1 + a(i,l) * e(i)
        end do

        con = 0.5e+00 * sum1 / a(l,l)
        do i = l1, n
          e(i) = e(i) - con * a(i,l)
          do j = l1, i
            a(i,j) = a(i,j) - a(i,l) * e(j) - a(j,l) * e(i)
          end do
        end do

70    continue

80    d(nm1) = a(nm1,nm1)
      d(n) = a(n,n)
      e(nm1) = a(n,nm1)
c
c  if eigenvectors are required, initialize a.
c
100   if ( .not. vec ) go to 180
c
c  if the matrix was tridiagonal, set a equal to the identity matrix.
c
      if ( .not. trd .and. n .ne. 2 ) go to 130

      do i = 1, n
        do j = 1, n
          a(i,j) = 0.0e+00
        end do
        a(i,i) = 1.0e+00
      end do

      go to 180
c
c  if the matrix was not tridiagonal, multiply out the
c  transformations obtained in the householder reduction.c
c
130   a(n,n) = 1.0e+00
      a(nm1,nm1) = 1.0e+00
      a(nm1,n) = 0.0e+00

      do l = 1, nm2

        ll = nm2 - l + 1
        ll1 = ll + 1

        do i = ll1, n
          sum = 0.0e+00
          do j = ll1, n
            sum = sum + a(j,ll) * a(j,i)
          end do
          a(ll,i) = sum / a(ll,ll)
        end do

        do i = ll1, n
          do j = ll1, n
            a(i,j) = a(i,j) - a(i,ll) * a(ll,j)
          end do
        end do

        do i = ll1, n
          a(i,ll) = 0.0e+00
          a(ll,i) = 0.0e+00
        end do

        a(ll,ll) = 1.0e+00

      end do
c
c  if an absolute convergence criterion is requested,
c  ( abscnv = .true. ), compute the infinity norm of the matrix.
c
180   if ( .not. abscnv ) go to 200
      ninf = amax1 ( 
     &  ninf, abs ( d(1) ) + abs ( e(1) ) + abs ( e(i-1) ) )
      if ( n .eq. 2 ) go to 200

      do i = 2, nm1
        ninf = amax1 ( ninf, 
     &    abs ( d(i) ) + abs ( e(i) ) + abs ( e(i-1) ) )
      end do
c
c  start the qr iteration.
c
200   nu = n
      num1 = n - 1
      shft = .false.
      k1 = k0
      test = ninf * eps
      e(n) = 0.0e+00
c
c  check for convergence and locate the submatrix in which the
c  qr step is to be performed.
c
210   do 220 nnl = 1, num1
        nl = num1 - nnl + 1
        if ( .not. abscnv ) 
     &    test = eps * amin1 ( abs ( d(nl) ), abs ( d(nl+1) ) )
        if ( abs ( e(nl) ) .le. test ) go to 230
220   continue

      go to 240
230   e(nl) = 0.0e+00
      nl = nl + 1
      if ( nl .ne. nu ) go to 240
      if ( num1 .eq. 1 ) return
      nu = num1
      num1 = nu - 1
      go to 210

240   e(nu) = e(nu) + 1.0e+00
      if ( e(nu) .le. titter ) go to 250
      fail = nu
      return
c
c  calculate the shift.
c
250   cb = ( d(num1) - d(nu) ) / 2.0e+00
      max = amax1 ( abs ( cb ), abs ( e(num1) ) )
      cb = cb / max
      cc = ( e(num1) / max )**2
      cd = sqrt ( cb**2 + cc )
      if ( cb .ne. 0.0e+00 ) cd = sign ( cd, cb )
      k2 = d(nu) - max * cc / ( cb + cd )
      if ( shft ) go to 270
      if ( abs ( k2 - k1 ) .lt. 0.5e+00 * abs ( k2 ) ) go to 260
      k1 = k2
      k = k0
      go to 300
260   shft = .true.
270   k = k2
c
c  perform one qr step with shift k on rows and columns
c  nl through nu.
c
300   p = d(nl) - k
      q = e(nl)
      call sincos ( p, q, c, s, norm )

      do 380 i = nl, num1
c
c  if required, rotate the eigenvectors.
c
        if ( .not. vec ) go to 330

        do j = 1, n
          temp = c * a(j,i) + s * a(j,i+1)
          a(j,i+1) = -s * a(j,i) + c*a(j,i+1)
          a(j,i) = temp
        end do
c
c  perform the similarity transformation and calculate the next
c  rotation.
c
330     d(i) = c * d(i) + s * e(i)
        temp = c * e(i) + s * d(i+1)
        d(i+1) = -s * e(i) + c * d(i+1)
        e(i) = -s * k
        d(i) = c * d(i) + s * temp
        if ( i .eq. num1 ) go to 380
        if ( abs ( s ) .gt. abs ( c ) ) go to 350
        r = s / c
        d(i+1) = -s * e(i) + c * d(i+1)
        p = d(i+1) - k
        q = c * e(i+1)
        call sincos ( p, q, c, s, norm )
        e(i) = r * norm
        e(i+1) = q
        go to 380
350     p = c * e(i) + s * d(i+1)
        q = s * e(i+1)
        d(i+1) = c * p / s + k
        e(i+1) = c * e(i+1)
        call sincos ( p, q, c, s, norm )
        e(i) = norm

380   continue

      temp = c * e(num1) + s * d(nu)
      d(nu) = -s * e(num1) + c * d(nu)
      e(num1) = temp
      go to 210

      end
      subroutine sincos ( p, q, c, s, norm )

c*********************************************************************72
c
c  sincos calculates the rotation corresponding to the vector (p,q).
c
      implicit none

      real c
      real norm
      real p
      real pp
      real q
      real qq
      real s

      pp = abs ( p )
      qq = abs ( q )
      if ( qq .gt. pp ) go to 510
      norm = pp * sqrt ( 1.0e+00 + ( qq / pp )**2 )
      go to 520
510   if ( qq .eq. 0.0e+00 ) go to 530
      norm = qq * sqrt ( 1.0e+00 + ( pp / qq )**2 )
520   c = p / norm
      s = q / norm
      return
530   c = 1.0e+00
      s = 0.0e+00
      norm = 0.0e+00
      return
      end