
Caveat receptor.  (Jack) dongarra@cs.utk.edu, (Eric Grosse) research!ehg
Careful! Anything free comes with no guarantee.
*
*PLEASE NOTE THAT netlib HAS MOVED, THE NEW ADDRESS IS netlib@ornl.gov.
*THE OLD ADDRESS, netlib@mcs.anl.gov, WILL BE TURNED OFF SOON.
*
*** from netlib, Mon Jun  4 12:50:06 EDT 1990 ***
      subroutine bqr(nm,n,mb,a,t,r,ierr,nv,rv)
c
      integer i,j,k,l,m,n,ii,ik,jk,jm,kj,kk,km,ll,mb,mk,mn,mz,
     x        m1,m2,m3,m4,ni,nm,nv,its,kj1,m21,m31,ierr,imult
      double precision a(nm,mb),rv(nv)
      double precision f,g,q,r,s,t,tst1,tst2,scale,pythag
c
c     this subroutine is a translation of the algol procedure bqr,
c     num. math. 16, 85-92(1970) by martin, reinsch, and wilkinson.
c     handbook for auto. comp., vol ii-linear algebra, 266-272(1971).
c
c     this subroutine finds the eigenvalue of smallest (usually)
c     magnitude of a real symmetric band matrix using the
c     qr algorithm with shifts of origin.  consecutive calls
c     can be made to find further eigenvalues.
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        mb is the (half) band width of the matrix, defined as the
c          number of adjacent diagonals, including the principal
c          diagonal, required to specify the non-zero portion of the
c          lower triangle of the matrix.
c
c        a contains the lower triangle of the symmetric band input
c          matrix stored as an n by mb array.  its lowest subdiagonal
c          is stored in the last n+1-mb positions of the first column,
c          its next subdiagonal in the last n+2-mb positions of the
c          second column, further subdiagonals similarly, and finally
c          its principal diagonal in the n positions of the last column.
c          contents of storages not part of the matrix are arbitrary.
c          on a subsequent call, its output contents from the previous
c          call should be passed.
c
c        t specifies the shift (of eigenvalues) applied to the diagonal
c          of a in forming the input matrix. what is actually determined
c          is the eigenvalue of a+ti (i is the identity matrix) nearest
c          to t.  on a subsequent call, the output value of t from the
c          previous call should be passed if the next nearest eigenvalue
c          is sought.
c
c        r should be specified as zero on the first call, and as its
c          output value from the previous call on a subsequent call.
c          it is used to determine when the last row and column of
c          the transformed band matrix can be regarded as negligible.
c
c        nv must be set to the dimension of the array parameter rv
c          as declared in the calling program dimension statement.
c
c     on output
c
c        a contains the transformed band matrix.  the matrix a+ti
c          derived from the output parameters is similar to the
c          input a+ti to within rounding errors.  its last row and
c          column are null (if ierr is zero).
c
c        t contains the computed eigenvalue of a+ti (if ierr is zero).
c
c        r contains the maximum of its input value and the norm of the
c          last column of the input matrix a.
c
c        ierr is set to
c          zero       for normal return,
c          n          if the eigenvalue has not been
c                     determined after 30 iterations.
c
c        rv is a temporary storage array of dimension at least
c          (2*mb**2+4*mb-3).  the first (3*mb-2) locations correspond
c          to the algol array b, the next (2*mb-1) locations correspond
c          to the algol array h, and the final (2*mb**2-mb) locations
c          correspond to the mb by (2*mb-1) algol array u.
c
c     note. for a subsequent call, n should be replaced by n-1, but
c     mb should not be altered even when it exceeds the current n.
c
c     calls pythag for  dsqrt(a*a + b*b) .
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
      m1 = min0(mb,n)
      m = m1 - 1
      m2 = m + m
      m21 = m2 + 1
      m3 = m21 + m
      m31 = m3 + 1
      m4 = m31 + m2
      mn = m + n
      mz = mb - m1
      its = 0
c     .......... test for convergence ..........
   40 g = a(n,mb)
      if (m .eq. 0) go to 360
      f = 0.0d0
c
      do 50 k = 1, m
         mk = k + mz
         f = f + dabs(a(n,mk))
   50 continue
c
      if (its .eq. 0 .and. f .gt. r) r = f
      tst1 = r
      tst2 = tst1 + f
      if (tst2 .le. tst1) go to 360
      if (its .eq. 30) go to 1000
      its = its + 1
c     .......... form shift from bottom 2 by 2 minor ..........
      if (f .gt. 0.25d0 * r .and. its .lt. 5) go to 90
      f = a(n,mb-1)
      if (f .eq. 0.0d0) go to 70
      q = (a(n-1,mb) - g) / (2.0d0 * f)
      s = pythag(q,1.0d0)
      g = g - f / (q + dsign(s,q))
   70 t = t + g
c
      do 80 i = 1, n
   80 a(i,mb) = a(i,mb) - g
c
   90 do 100 k = m31, m4
  100 rv(k) = 0.0d0
c
      do 350 ii = 1, mn
         i = ii - m
         ni = n - ii
         if (ni .lt. 0) go to 230
c     .......... form column of shifted matrix a-g*i ..........
         l = max0(1,2-i)
c
         do 110 k = 1, m3
  110    rv(k) = 0.0d0
c
         do 120 k = l, m1
            km = k + m
            mk = k + mz
            rv(km) = a(ii,mk)
  120    continue
c
         ll = min0(m,ni)
         if (ll .eq. 0) go to 135
c
         do 130 k = 1, ll
            km = k + m21
            ik = ii + k
            mk = mb - k
            rv(km) = a(ik,mk)
  130    continue
c     .......... pre-multiply with householder reflections ..........
  135    ll = m2
         imult = 0
c     .......... multiplication procedure ..........
  140    kj = m4 - m1
c
         do 170 j = 1, ll
            kj = kj + m1
            jm = j + m3
            if (rv(jm) .eq. 0.0d0) go to 170
            f = 0.0d0
c
            do 150 k = 1, m1
               kj = kj + 1
               jk = j + k - 1
               f = f + rv(kj) * rv(jk)
  150       continue
c
            f = f / rv(jm)
            kj = kj - m1
c
            do 160 k = 1, m1
               kj = kj + 1
               jk = j + k - 1
               rv(jk) = rv(jk) - rv(kj) * f
  160       continue
c
            kj = kj - m1
  170    continue
c
         if (imult .ne. 0) go to 280
c     .......... householder reflection ..........
         f = rv(m21)
         s = 0.0d0
         rv(m4) = 0.0d0
         scale = 0.0d0
c
         do 180 k = m21, m3
  180    scale = scale + dabs(rv(k))
c
         if (scale .eq. 0.0d0) go to 210
c
         do 190 k = m21, m3
  190    s = s + (rv(k)/scale)**2
c
         s = scale * scale * s
         g = -dsign(dsqrt(s),f)
         rv(m21) = g
         rv(m4) = s - f * g
         kj = m4 + m2 * m1 + 1
         rv(kj) = f - g
c
         do 200 k = 2, m1
            kj = kj + 1
            km = k + m2
            rv(kj) = rv(km)
  200    continue
c     .......... save column of triangular factor r ..........
  210    do 220 k = l, m1
            km = k + m
            mk = k + mz
            a(ii,mk) = rv(km)
  220    continue
c
  230    l = max0(1,m1+1-i)
         if (i .le. 0) go to 300
c     .......... perform additional steps ..........
         do 240 k = 1, m21
  240    rv(k) = 0.0d0
c
         ll = min0(m1,ni+m1)
c     .......... get row of triangular factor r ..........
         do 250 kk = 1, ll
            k = kk - 1
            km = k + m1
            ik = i + k
            mk = mb - k
            rv(km) = a(ik,mk)
  250    continue
c     .......... post-multiply with householder reflections ..........
         ll = m1
         imult = 1
         go to 140
c     .......... store column of new a matrix ..........
  280    do 290 k = l, m1
            mk = k + mz
            a(i,mk) = rv(k)
  290    continue
c     .......... update householder reflections ..........
  300    if (l .gt. 1) l = l - 1
         kj1 = m4 + l * m1
c
         do 320 j = l, m2
            jm = j + m3
            rv(jm) = rv(jm+1)
c
            do 320 k = 1, m1
               kj1 = kj1 + 1
               kj = kj1 - m1
               rv(kj) = rv(kj1)
  320    continue
c
  350 continue
c
      go to 40
c     .......... convergence ..........
  360 t = t + g
c
      do 380 i = 1, n
  380 a(i,mb) = a(i,mb) - g
c
      do 400 k = 1, m1
         mk = k + mz
         a(n,mk) = 0.0d0
  400 continue
c
      go to 1001
c     .......... set error -- no convergence to
c                eigenvalue after 30 iterations ..........
 1000 ierr = n
 1001 return
      end
