c-----------------------------------------------------------------------
c\BeginDoc
c
c\Name: dlaqrb
c
c\Description:
c  Compute the eigenvalues and the Schur decomposition of an upper 
c  Hessenberg submatrix in rows and columns ILO to IHI.  Only the
c  last component of the Schur vectors are computed.
c
c  This is mostly a modification of the LAPACK routine dlahqr.
c  
c\Usage:
c  call dlaqrb
c     ( WANTT, N, ILO, IHI, H, LDH, WR, WI,  Z, INFO )
c
c\Arguments
c  WANTT   Logical variable.  (INPUT)
c          = .TRUE. : the full Schur form T is required;
c          = .FALSE.: only eigenvalues are required.
c
c  N       Integer.  (INPUT)
c          The order of the matrix H.  N >= 0.
c
c  ILO     Integer.  (INPUT)
c  IHI     Integer.  (INPUT)
c          It is assumed that H is already upper quasi-triangular in
c          rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
c          ILO = 1). SLAQRB works primarily with the Hessenberg
c          submatrix in rows and columns ILO to IHI, but applies
c          transformations to all of H if WANTT is .TRUE..
c          1 <= ILO <= max(1,IHI); IHI <= N.
c
c  H       Double precision array, dimension (LDH,N).  (INPUT/OUTPUT)
c          On entry, the upper Hessenberg matrix H.
c          On exit, if WANTT is .TRUE., H is upper quasi-triangular in
c          rows and columns ILO:IHI, with any 2-by-2 diagonal blocks in
c          standard form. If WANTT is .FALSE., the contents of H are
c          unspecified on exit.
c
c  LDH     Integer.  (INPUT)
c          The leading dimension of the array H. LDH >= max(1,N).
c
c  WR      Double precision array, dimension (N).  (OUTPUT)
c  WI      Double precision array, dimension (N).  (OUTPUT)
c          The real and imaginary parts, respectively, of the computed
c          eigenvalues ILO to IHI are stored in the corresponding
c          elements of WR and WI. If two eigenvalues are computed as a
c          complex conjugate pair, they are stored in consecutive
c          elements of WR and WI, say the i-th and (i+1)th, with
c          WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the
c          eigenvalues are stored in the same order as on the diagonal
c          of the Schur form returned in H, with WR(i) = H(i,i), and, if
c          H(i:i+1,i:i+1) is a 2-by-2 diagonal block,
c          WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).
c
c  Z       Double precision array, dimension (N).  (OUTPUT)
c          On exit Z contains the last components of the Schur vectors.
c
c  INFO    Integer.  (OUPUT)
c          = 0: successful exit
c          > 0: SLAQRB failed to compute all the eigenvalues ILO to IHI
c               in a total of 30*(IHI-ILO+1) iterations; if INFO = i,
c               elements i+1:ihi of WR and WI contain those eigenvalues
c               which have been successfully computed.
c
c\Remarks
c  1. None.
c
c-----------------------------------------------------------------------
c
c\BeginLib
c
c\Local variables:
c     xxxxxx  real
c
c\Routines called:
c     dlabad  LAPACK routine that computes machine constants.
c     dlamch  LAPACK routine that determines machine constants.
c     dlanhs  LAPACK routine that computes various norms of a matrix.
c     dlanv2  LAPACK routine that computes the Schur factorization of
c             2 by 2 nonsymmetric matrix in standard form.
c     dlarfg  LAPACK Householder reflection construction routine.
c     dcopy   Level 1 BLAS that copies one vector to another.
c     drot    Level 1 BLAS that applies a rotation to a 2 by 2 matrix.

c
c\Author
c     Danny Sorensen               Phuong Vu
c     Richard Lehoucq              CRPC / Rice University
c     Dept. of Computational &     Houston, Texas 
c     Applied Mathematics
c     Rice University           
c     Houston, Texas            
c
c\Revision history:
c     xx/xx/92: Version ' 2.4'
c               Modified from the LAPACK routine dlahqr so that only the
c               last component of the Schur vectors are computed.
c
c\SCCS Information: @(#) 
c FILE: laqrb.F   SID: 2.2   DATE OF SID: 8/27/96   RELEASE: 2
c
c\Remarks
c     1. None
c
c\EndLib
c
c-----------------------------------------------------------------------
c
      subroutine dlaqrb ( wantt, n, ilo, ihi, h, ldh, wr, wi,
     &                    z, info )
c
c     %------------------%
c     | Scalar Arguments |
c     %------------------%
c
      logical    wantt
      integer    ihi, ilo, info, ldh, n
c
c     %-----------------%
c     | Array Arguments |
c     %-----------------%
c
      Double precision
     &           h( ldh, * ), wi( * ), wr( * ), z( * )
c
c     %------------%
c     | Parameters |
c     %------------%
c
      Double precision
     &           zero, one, dat1, dat2
      parameter (zero = 0.0D+0, one = 1.0D+0, dat1 = 7.5D-1, 
     &           dat2 = -4.375D-1)
c
c     %------------------------%
c     | Local Scalars & Arrays |
c     %------------------------%
c
      integer    i, i1, i2, itn, its, j, k, l, m, nh, nr
      Double precision
     &           cs, h00, h10, h11, h12, h21, h22, h33, h33s,
     &           h43h34, h44, h44s, ovfl, s, smlnum, sn, sum,
     &           t1, t2, t3, tst1, ulp, unfl, v1, v2, v3
      Double precision
     &           v( 3 ), work( 1 )
c
c     %--------------------%
c     | External Functions |
c     %--------------------%
c
      Double precision
     &           dlamch, dlanhs
      external   dlamch, dlanhs
c
c     %----------------------%
c     | External Subroutines |
c     %----------------------%
c
      external   dcopy, dlabad, dlanv2, dlarfg, drot
c
c     %-----------------------%
c     | Executable Statements |
c     %-----------------------%
c
      info = 0
c
c     %--------------------------%
c     | Quick return if possible |
c     %--------------------------%
c
      if( n.eq.0 )
     &   return
      if( ilo.eq.ihi ) then
         wr( ilo ) = h( ilo, ilo )
         wi( ilo ) = zero
         return
      end if
c 
c     %---------------------------------------------%
c     | Initialize the vector of last components of |
c     | the Schur vectors for accumulation.         |
c     %---------------------------------------------%
c
      do 5 j = 1, n-1
         z(j) = zero
  5   continue 
      z(n) = one
c 
      nh = ihi - ilo + 1
c
c     %-------------------------------------------------------------%
c     | Set machine-dependent constants for the stopping criterion. |
c     | If norm(H) <= sqrt(OVFL), overflow should not occur.        |
c     %-------------------------------------------------------------%
c
      unfl = dlamch( 'safe minimum' )
      ovfl = one / unfl
      call dlabad( unfl, ovfl )
      ulp = dlamch( 'precision' )
      smlnum = unfl*( nh / ulp )
c
c     %---------------------------------------------------------------%
c     | I1 and I2 are the indices of the first row and last column    |
c     | of H to which transformations must be applied. If eigenvalues |
c     | only are computed, I1 and I2 are set inside the main loop.    |
c     | Zero out H(J+2,J) = ZERO for J=1:N if WANTT = .TRUE.          |
c     | else H(J+2,J) for J=ILO:IHI-ILO-1 if WANTT = .FALSE.          |
c     %---------------------------------------------------------------%
c
      if( wantt ) then
         i1 = 1
         i2 = n
         do 8 i=1,i2-2
            h(i1+i+1,i) = zero
 8       continue
      else
         do 9 i=1, ihi-ilo-1
            h(ilo+i+1,ilo+i-1) = zero
 9       continue
      end if
c 
c     %---------------------------------------------------%
c     | ITN is the total number of QR iterations allowed. |
c     %---------------------------------------------------%
c
      itn = 30*nh
c 
c     ------------------------------------------------------------------
c     The main loop begins here. I is the loop index and decreases from
c     IHI to ILO in steps of 1 or 2. Each iteration of the loop works
c     with the active submatrix in rows and columns L to I.
c     Eigenvalues I+1 to IHI have already converged. Either L = ILO or
c     H(L,L-1) is negligible so that the matrix splits.
c     ------------------------------------------------------------------
c 
      i = ihi
   10 continue
      l = ilo
      if( i.lt.ilo )
     &   go to 150
 
c     %--------------------------------------------------------------%
c     | Perform QR iterations on rows and columns ILO to I until a   |
c     | submatrix of order 1 or 2 splits off at the bottom because a |
c     | subdiagonal element has become negligible.                   |
c     %--------------------------------------------------------------%
 
      do 130 its = 0, itn
c
c        %----------------------------------------------%
c        | Look for a single small subdiagonal element. |
c        %----------------------------------------------%
c
         do 20 k = i, l + 1, -1
            tst1 = abs( h( k-1, k-1 ) ) + abs( h( k, k ) )
            if( tst1.eq.zero )
     &         tst1 = dlanhs( '1', i-l+1, h( l, l ), ldh, work )
            if( abs( h( k, k-1 ) ).le.max( ulp*tst1, smlnum ) )
     &         go to 30
   20    continue
   30    continue
         l = k
         if( l.gt.ilo ) then
c
c           %------------------------%
c           | H(L,L-1) is negligible |
c           %------------------------%
c
            h( l, l-1 ) = zero
         end if
c
c        %-------------------------------------------------------------%
c        | Exit from loop if a submatrix of order 1 or 2 has split off |
c        %-------------------------------------------------------------%
c
         if( l.ge.i-1 )
     &      go to 140
c
c        %---------------------------------------------------------%
c        | Now the active submatrix is in rows and columns L to I. |
c        | If eigenvalues only are being computed, only the active |
c        | submatrix need be transformed.                          |
c        %---------------------------------------------------------%
c
         if( .not.wantt ) then
            i1 = l
            i2 = i
         end if
c 
         if( its.eq.10 .or. its.eq.20 ) then
c
c           %-------------------%
c           | Exceptional shift |
c           %-------------------%
c
            s = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
            h44 = dat1*s
            h33 = h44
            h43h34 = dat2*s*s
c
         else
c
c           %-----------------------------------------%
c           | Prepare to use Wilkinson's double shift |
c           %-----------------------------------------%
c
            h44 = h( i, i )
            h33 = h( i-1, i-1 )
            h43h34 = h( i, i-1 )*h( i-1, i )
         end if
c
c        %-----------------------------------------------------%
c        | Look for two consecutive small subdiagonal elements |
c        %-----------------------------------------------------%
c
         do 40 m = i - 2, l, -1
c
c           %---------------------------------------------------------%
c           | Determine the effect of starting the double-shift QR    |
c           | iteration at row M, and see if this would make H(M,M-1) |
c           | negligible.                                             |
c           %---------------------------------------------------------%
c
            h11 = h( m, m )
            h22 = h( m+1, m+1 )
            h21 = h( m+1, m )
            h12 = h( m, m+1 )
            h44s = h44 - h11
            h33s = h33 - h11
            v1 = ( h33s*h44s-h43h34 ) / h21 + h12
            v2 = h22 - h11 - h33s - h44s
            v3 = h( m+2, m+1 )
            s = abs( v1 ) + abs( v2 ) + abs( v3 )
            v1 = v1 / s
            v2 = v2 / s
            v3 = v3 / s
            v( 1 ) = v1
            v( 2 ) = v2
            v( 3 ) = v3
            if( m.eq.l )
     &         go to 50
            h00 = h( m-1, m-1 )
            h10 = h( m, m-1 )
            tst1 = abs( v1 )*( abs( h00 )+abs( h11 )+abs( h22 ) )
            if( abs( h10 )*( abs( v2 )+abs( v3 ) ).le.ulp*tst1 )
     &         go to 50
   40    continue
   50    continue
c
c        %----------------------%
c        | Double-shift QR step |
c        %----------------------%
c
         do 120 k = m, i - 1
c 
c           ------------------------------------------------------------
c           The first iteration of this loop determines a reflection G
c           from the vector V and applies it from left and right to H,
c           thus creating a nonzero bulge below the subdiagonal.
c
c           Each subsequent iteration determines a reflection G to
c           restore the Hessenberg form in the (K-1)th column, and thus
c           chases the bulge one step toward the bottom of the active
c           submatrix. NR is the order of G.
c           ------------------------------------------------------------
c 
            nr = min( 3, i-k+1 )
            if( k.gt.m )
     &         call dcopy( nr, h( k, k-1 ), 1, v, 1 )
            call dlarfg( nr, v( 1 ), v( 2 ), 1, t1 )
            if( k.gt.m ) then
               h( k, k-1 ) = v( 1 )
               h( k+1, k-1 ) = zero
               if( k.lt.i-1 )
     &            h( k+2, k-1 ) = zero
            else if( m.gt.l ) then
               h( k, k-1 ) = -h( k, k-1 )
            end if
            v2 = v( 2 )
            t2 = t1*v2
            if( nr.eq.3 ) then
               v3 = v( 3 )
               t3 = t1*v3
c
c              %------------------------------------------------%
c              | Apply G from the left to transform the rows of |
c              | the matrix in columns K to I2.                 |
c              %------------------------------------------------%
c
               do 60 j = k, i2
                  sum = h( k, j ) + v2*h( k+1, j ) + v3*h( k+2, j )
                  h( k, j ) = h( k, j ) - sum*t1
                  h( k+1, j ) = h( k+1, j ) - sum*t2
                  h( k+2, j ) = h( k+2, j ) - sum*t3
   60          continue
c
c              %----------------------------------------------------%
c              | Apply G from the right to transform the columns of |
c              | the matrix in rows I1 to min(K+3,I).               |
c              %----------------------------------------------------%
c
               do 70 j = i1, min( k+3, i )
                  sum = h( j, k ) + v2*h( j, k+1 ) + v3*h( j, k+2 )
                  h( j, k ) = h( j, k ) - sum*t1
                  h( j, k+1 ) = h( j, k+1 ) - sum*t2
                  h( j, k+2 ) = h( j, k+2 ) - sum*t3
   70          continue
c
c              %----------------------------------%
c              | Accumulate transformations for Z |
c              %----------------------------------%
c
               sum      = z( k ) + v2*z( k+1 ) + v3*z( k+2 )
               z( k )   = z( k ) - sum*t1
               z( k+1 ) = z( k+1 ) - sum*t2
               z( k+2 ) = z( k+2 ) - sum*t3
 
            else if( nr.eq.2 ) then
c
c              %------------------------------------------------%
c              | Apply G from the left to transform the rows of |
c              | the matrix in columns K to I2.                 |
c              %------------------------------------------------%
c
               do 90 j = k, i2
                  sum = h( k, j ) + v2*h( k+1, j )
                  h( k, j ) = h( k, j ) - sum*t1
                  h( k+1, j ) = h( k+1, j ) - sum*t2
   90          continue
c
c              %----------------------------------------------------%
c              | Apply G from the right to transform the columns of |
c              | the matrix in rows I1 to min(K+3,I).               |
c              %----------------------------------------------------%
c
               do 100 j = i1, i
                  sum = h( j, k ) + v2*h( j, k+1 )
                  h( j, k ) = h( j, k ) - sum*t1
                  h( j, k+1 ) = h( j, k+1 ) - sum*t2
  100          continue
c
c              %----------------------------------%
c              | Accumulate transformations for Z |
c              %----------------------------------%
c
               sum      = z( k ) + v2*z( k+1 )
               z( k )   = z( k ) - sum*t1
               z( k+1 ) = z( k+1 ) - sum*t2
            end if
  120    continue
 
  130 continue
c
c     %-------------------------------------------------------%
c     | Failure to converge in remaining number of iterations |
c     %-------------------------------------------------------%
c
      info = i
      return
 
  140 continue
 
      if( l.eq.i ) then
c
c        %------------------------------------------------------%
c        | H(I,I-1) is negligible: one eigenvalue has converged |
c        %------------------------------------------------------%
c
         wr( i ) = h( i, i )
         wi( i ) = zero

      else if( l.eq.i-1 ) then
c
c        %--------------------------------------------------------%
c        | H(I-1,I-2) is negligible;                              |
c        | a pair of eigenvalues have converged.                  |
c        |                                                        |
c        | Transform the 2-by-2 submatrix to standard Schur form, |
c        | and compute and store the eigenvalues.                 |
c        %--------------------------------------------------------%
c
         call dlanv2( h( i-1, i-1 ), h( i-1, i ), h( i, i-1 ),
     &                h( i, i ), wr( i-1 ), wi( i-1 ), wr( i ), wi( i ),
     &                cs, sn )
 
         if( wantt ) then
c
c           %-----------------------------------------------------%
c           | Apply the transformation to the rest of H and to Z, |
c           | as required.                                        |
c           %-----------------------------------------------------%
c
            if( i2.gt.i )
     &         call drot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,
     &                    cs, sn )
            call drot( i-i1-1, h( i1, i-1 ), 1, h( i1, i ), 1, cs, sn )
            sum      = cs*z( i-1 ) + sn*z( i )
            z( i )   = cs*z( i )   - sn*z( i-1 )
            z( i-1 ) = sum
         end if
      end if
c
c     %---------------------------------------------------------%
c     | Decrement number of remaining iterations, and return to |
c     | start of the main loop with new value of I.             |
c     %---------------------------------------------------------%
c
      itn = itn - its
      i = l - 1
      go to 10
 
  150 continue
      return
c
c     %---------------%
c     | End of dlaqrb |
c     %---------------%
c
      end