dlahqr.f

*> \brief \b DLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DLAHQR + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlahqr.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlahqr.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlahqr.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
*                          ILOZ, IHIZ, Z, LDZ, INFO )
*
*       .. Scalar Arguments ..
*       INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
*       LOGICAL            WANTT, WANTZ
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION   H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*>    DLAHQR is an auxiliary routine called by DHSEQR to update the
*>    eigenvalues and Schur decomposition already computed by DHSEQR, by
*>    dealing with the Hessenberg submatrix in rows and columns ILO to
*>    IHI.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] WANTT
*> \verbatim
*>          WANTT is LOGICAL
*>          = .TRUE. : the full Schur form T is required;
*>          = .FALSE.: only eigenvalues are required.
*> \endverbatim
*>
*> \param[in] WANTZ
*> \verbatim
*>          WANTZ is LOGICAL
*>          = .TRUE. : the matrix of Schur vectors Z is required;
*>          = .FALSE.: Schur vectors are not required.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix H.  N >= 0.
*> \endverbatim
*>
*> \param[in] ILO
*> \verbatim
*>          ILO is INTEGER
*> \endverbatim
*>
*> \param[in] IHI
*> \verbatim
*>          IHI is INTEGER
*>          It is assumed that H is already upper quasi-triangular in
*>          rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
*>          ILO = 1). DLAHQR works primarily with the Hessenberg
*>          submatrix in rows and columns ILO to IHI, but applies
*>          transformations to all of H if WANTT is .TRUE..
*>          1 <= ILO <= max(1,IHI); IHI <= N.
*> \endverbatim
*>
*> \param[in,out] H
*> \verbatim
*>          H is DOUBLE PRECISION array, dimension (LDH,N)
*>          On entry, the upper Hessenberg matrix H.
*>          On exit, if INFO is zero and if WANTT is .TRUE., H is upper
*>          quasi-triangular in rows and columns ILO:IHI, with any
*>          2-by-2 diagonal blocks in standard form. If INFO is zero
*>          and WANTT is .FALSE., the contents of H are unspecified on
*>          exit.  The output state of H if INFO is nonzero is given
*>          below under the description of INFO.
*> \endverbatim
*>
*> \param[in] LDH
*> \verbatim
*>          LDH is INTEGER
*>          The leading dimension of the array H. LDH >= max(1,N).
*> \endverbatim
*>
*> \param[out] WR
*> \verbatim
*>          WR is DOUBLE PRECISION array, dimension (N)
*> \endverbatim
*>
*> \param[out] WI
*> \verbatim
*>          WI is DOUBLE PRECISION array, dimension (N)
*>          The real and imaginary parts, respectively, of the computed
*>          eigenvalues ILO to IHI are stored in the corresponding
*>          elements of WR and WI. If two eigenvalues are computed as a
*>          complex conjugate pair, they are stored in consecutive
*>          elements of WR and WI, say the i-th and (i+1)th, with
*>          WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the
*>          eigenvalues are stored in the same order as on the diagonal
*>          of the Schur form returned in H, with WR(i) = H(i,i), and, if
*>          H(i:i+1,i:i+1) is a 2-by-2 diagonal block,
*>          WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).
*> \endverbatim
*>
*> \param[in] ILOZ
*> \verbatim
*>          ILOZ is INTEGER
*> \endverbatim
*>
*> \param[in] IHIZ
*> \verbatim
*>          IHIZ is INTEGER
*>          Specify the rows of Z to which transformations must be
*>          applied if WANTZ is .TRUE..
*>          1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
*> \endverbatim
*>
*> \param[in,out] Z
*> \verbatim
*>          Z is DOUBLE PRECISION array, dimension (LDZ,N)
*>          If WANTZ is .TRUE., on entry Z must contain the current
*>          matrix Z of transformations accumulated by DHSEQR, and on
*>          exit Z has been updated; transformations are applied only to
*>          the submatrix Z(ILOZ:IHIZ,ILO:IHI).
*>          If WANTZ is .FALSE., Z is not referenced.
*> \endverbatim
*>
*> \param[in] LDZ
*> \verbatim
*>          LDZ is INTEGER
*>          The leading dimension of the array Z. LDZ >= max(1,N).
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>           = 0:  successful exit
*>           > 0:  If INFO = i, DLAHQR failed to compute all the
*>                  eigenvalues ILO to IHI in a total of 30 iterations
*>                  per eigenvalue; elements i+1:ihi of WR and WI
*>                  contain those eigenvalues which have been
*>                  successfully computed.
*>
*>                  If INFO > 0 and WANTT is .FALSE., then on exit,
*>                  the remaining unconverged eigenvalues are the
*>                  eigenvalues of the upper Hessenberg matrix rows
*>                  and columns ILO through INFO of the final, output
*>                  value of H.
*>
*>                  If INFO > 0 and WANTT is .TRUE., then on exit
*>          (*)       (initial value of H)*U  = U*(final value of H)
*>                  where U is an orthogonal matrix.    The final
*>                  value of H is upper Hessenberg and triangular in
*>                  rows and columns INFO+1 through IHI.
*>
*>                  If INFO > 0 and WANTZ is .TRUE., then on exit
*>                      (final value of Z)  = (initial value of Z)*U
*>                  where U is the orthogonal matrix in (*)
*>                  (regardless of the value of WANTT.)
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup lahqr
*
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>     02-96 Based on modifications by
*>     David Day, Sandia National Laboratory, USA
*>
*>     12-04 Further modifications by
*>     Ralph Byers, University of Kansas, USA
*>     This is a modified version of DLAHQR from LAPACK version 3.0.
*>     It is (1) more robust against overflow and underflow and
*>     (2) adopts the more conservative Ahues & Tisseur stopping
*>     criterion (LAWN 122, 1997).
*> \endverbatim
*>
*  =====================================================================
      SUBROUTINE dlahqr( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
     $                   ILOZ, IHIZ, Z, LDZ, INFO )
      IMPLICIT NONE
*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
      LOGICAL            WANTT, WANTZ
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   H( ldh, * ), WI( * ), WR( * ), Z( ldz, * )
*     ..
*
*  =========================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, TWO
      parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
      DOUBLE PRECISION   DAT1, DAT2
      parameter( dat1 = 3.0d0 / 4.0d0, dat2 = -0.4375d0 )
      INTEGER            KEXSH
      parameter( kexsh = 10 )
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION   AA, AB, BA, BB, CS, DET, H11, H12, H21, H21S,
     $                   h22, rt1i, rt1r, rt2i, rt2r, rtdisc, s, safmax,
     $                   safmin, smlnum, sn, sum, t1, t2, t3, tr, tst,
     $                   ulp, v2, v3
      INTEGER            I, I1, I2, ITS, ITMAX, J, K, L, M, NH, NR, NZ,
     $                   kdefl 
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION   V( 3 )
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           dlamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           dcopy, dlabad, dlanv2, dlarfg, drot
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, max, min, sqrt
*     ..
*     .. Executable Statements ..
*
      info = 0
*
*     Quick return if possible
*
      IF( n.EQ.0 )
     $   RETURN
      IF( ilo.EQ.ihi ) THEN
         wr( ilo ) = h( ilo, ilo )
         wi( ilo ) = zero
         RETURN
      END IF
*
*     ==== clear out the trash ====
      DO 10 j = ilo, ihi - 3
         h( j+2, j ) = zero
         h( j+3, j ) = zero
   10 CONTINUE
      IF( ilo.LE.ihi-2 )
     $   h( ihi, ihi-2 ) = zero
*
      nh = ihi - ilo + 1
      nz = ihiz - iloz + 1
*
*     Set machine-dependent constants for the stopping criterion.
*
      safmin = dlamch( 'SAFE MINIMUM' )
      safmax = one / safmin
      CALL dlabad( safmin, safmax )
      ulp = dlamch( 'PRECISION' )
      smlnum = safmin*( dble( nh ) / ulp )
*
*     I1 and I2 are the indices of the first row and last column of H
*     to which transformations must be applied. If eigenvalues only are
*     being computed, I1 and I2 are set inside the main loop.
*
      IF( wantt ) THEN
         i1 = 1
         i2 = n
      END IF
*
*     ITMAX is the total number of QR iterations allowed.
*
      itmax = 30 * max( 10, nh )
*
*     KDEFL counts the number of iterations since a deflation
*
      kdefl = 0
*
*     The main loop begins here. I is the loop index and decreases from
*     IHI to ILO in steps of 1 or 2. Each iteration of the loop works
*     with the active submatrix in rows and columns L to I.
*     Eigenvalues I+1 to IHI have already converged. Either L = ILO or
*     H(L,L-1) is negligible so that the matrix splits.
*
      i = ihi
   20 CONTINUE
      l = ilo
      IF( i.LT.ilo )
     $   GO TO 160
*
*     Perform QR iterations on rows and columns ILO to I until a
*     submatrix of order 1 or 2 splits off at the bottom because a
*     subdiagonal element has become negligible.
*
      DO 140 its = 0, itmax
*
*        Look for a single small subdiagonal element.
*
         DO 30 k = i, l + 1, -1
            IF( abs( h( k, k-1 ) ).LE.smlnum )
     $         GO TO 40
            tst = abs( h( k-1, k-1 ) ) + abs( h( k, k ) )
            IF( tst.EQ.zero ) THEN
               IF( k-2.GE.ilo )
     $            tst = tst + abs( h( k-1, k-2 ) )
               IF( k+1.LE.ihi )
     $            tst = tst + abs( h( k+1, k ) )
            END IF
*           ==== The following is a conservative small subdiagonal
*           .    deflation  criterion due to Ahues & Tisseur (LAWN 122,
*           .    1997). It has better mathematical foundation and
*           .    improves accuracy in some cases.  ====
            IF( abs( h( k, k-1 ) ).LE.ulp*tst ) THEN
               ab = max( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) )
               ba = min( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) )
               aa = max( abs( h( k, k ) ),
     $              abs( h( k-1, k-1 )-h( k, k ) ) )
               bb = min( abs( h( k, k ) ),
     $              abs( h( k-1, k-1 )-h( k, k ) ) )
               s = aa + ab
               IF( ba*( ab / s ).LE.max( smlnum,
     $             ulp*( bb*( aa / s ) ) ) )GO TO 40
            END IF
   30    CONTINUE
   40    CONTINUE
         l = k
         IF( l.GT.ilo ) THEN
*
*           H(L,L-1) is negligible
*
            h( l, l-1 ) = zero
         END IF
*
*        Exit from loop if a submatrix of order 1 or 2 has split off.
*
         IF( l.GE.i-1 )
     $      GO TO 150
         kdefl = kdefl + 1
*
*        Now the active submatrix is in rows and columns L to I. If
*        eigenvalues only are being computed, only the active submatrix
*        need be transformed.
*
         IF( .NOT.wantt ) THEN
            i1 = l
            i2 = i
         END IF
*
         IF( mod(kdefl,2*kexsh).EQ.0 ) THEN
*
*           Exceptional shift.
*
            s = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
            h11 = dat1*s + h( i, i )
            h12 = dat2*s
            h21 = s
            h22 = h11
         ELSE IF( mod(kdefl,kexsh).EQ.0 ) THEN
*
*           Exceptional shift.
*
            s = abs( h( l+1, l ) ) + abs( h( l+2, l+1 ) )
            h11 = dat1*s + h( l, l )
            h12 = dat2*s
            h21 = s
            h22 = h11
         ELSE
*
*           Prepare to use Francis' double shift
*           (i.e. 2nd degree generalized Rayleigh quotient)
*
            h11 = h( i-1, i-1 )
            h21 = h( i, i-1 )
            h12 = h( i-1, i )
            h22 = h( i, i )
         END IF
         s = abs( h11 ) + abs( h12 ) + abs( h21 ) + abs( h22 )
         IF( s.EQ.zero ) THEN
            rt1r = zero
            rt1i = zero
            rt2r = zero
            rt2i = zero
         ELSE
            h11 = h11 / s
            h21 = h21 / s
            h12 = h12 / s
            h22 = h22 / s
            tr = ( h11+h22 ) / two
            det = ( h11-tr )*( h22-tr ) - h12*h21
            rtdisc = sqrt( abs( det ) )
            IF( det.GE.zero ) THEN
*
*              ==== complex conjugate shifts ====
*
               rt1r = tr*s
               rt2r = rt1r
               rt1i = rtdisc*s
               rt2i = -rt1i
            ELSE
*
*              ==== real shifts (use only one of them)  ====
*
               rt1r = tr + rtdisc
               rt2r = tr - rtdisc
               IF( abs( rt1r-h22 ).LE.abs( rt2r-h22 ) ) THEN
                  rt1r = rt1r*s
                  rt2r = rt1r
               ELSE
                  rt2r = rt2r*s
                  rt1r = rt2r
               END IF
               rt1i = zero
               rt2i = zero
            END IF
         END IF
*
*        Look for two consecutive small subdiagonal elements.
*
         DO 50 m = i - 2, l, -1
*           Determine the effect of starting the double-shift QR
*           iteration at row M, and see if this would make H(M,M-1)
*           negligible.  (The following uses scaling to avoid
*           overflows and most underflows.)
*
            h21s = h( m+1, m )
            s = abs( h( m, m )-rt2r ) + abs( rt2i ) + abs( h21s )
            h21s = h( m+1, m ) / s
            v( 1 ) = h21s*h( m, m+1 ) + ( h( m, m )-rt1r )*
     $               ( ( h( m, m )-rt2r ) / s ) - rt1i*( rt2i / s )
            v( 2 ) = h21s*( h( m, m )+h( m+1, m+1 )-rt1r-rt2r )
            v( 3 ) = h21s*h( m+2, m+1 )
            s = abs( v( 1 ) ) + abs( v( 2 ) ) + abs( v( 3 ) )
            v( 1 ) = v( 1 ) / s
            v( 2 ) = v( 2 ) / s
            v( 3 ) = v( 3 ) / s
            IF( m.EQ.l )
     $         GO TO 60
            IF( abs( h( m, m-1 ) )*( abs( v( 2 ) )+abs( v( 3 ) ) ).LE.
     $          ulp*abs( v( 1 ) )*( abs( h( m-1, m-1 ) )+abs( h( m,
     $          m ) )+abs( h( m+1, m+1 ) ) ) )GO TO 60
   50    CONTINUE
   60    CONTINUE
*
*        Double-shift QR step
*
         DO 130 k = m, i - 1
*
*           The first iteration of this loop determines a reflection G
*           from the vector V and applies it from left and right to H,
*           thus creating a nonzero bulge below the subdiagonal.
*
*           Each subsequent iteration determines a reflection G to
*           restore the Hessenberg form in the (K-1)th column, and thus
*           chases the bulge one step toward the bottom of the active
*           submatrix. NR is the order of G.
*
            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
*               ==== Use the following instead of
*               .    H( K, K-1 ) = -H( K, K-1 ) to
*               .    avoid a bug when v(2) and v(3)
*               .    underflow. ====
               h( k, k-1 ) = h( k, k-1 )*( one-t1 )
            END IF
            v2 = v( 2 )
            t2 = t1*v2
            IF( nr.EQ.3 ) THEN
               v3 = v( 3 )
               t3 = t1*v3
*
*              Apply G from the left to transform the rows of the matrix
*              in columns K to I2.
*
               DO 70 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
   70          CONTINUE
*
*              Apply G from the right to transform the columns of the
*              matrix in rows I1 to min(K+3,I).
*
               DO 80 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
   80          CONTINUE
*
               IF( wantz ) THEN
*
*                 Accumulate transformations in the matrix Z
*
                  DO 90 j = iloz, ihiz
                     sum = z( j, k ) + v2*z( j, k+1 ) + v3*z( j, k+2 )
                     z( j, k ) = z( j, k ) - sum*t1
                     z( j, k+1 ) = z( j, k+1 ) - sum*t2
                     z( j, k+2 ) = z( j, k+2 ) - sum*t3
   90             CONTINUE
               END IF
            ELSE IF( nr.EQ.2 ) THEN
*
*              Apply G from the left to transform the rows of the matrix
*              in columns K to I2.
*
               DO 100 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
  100          CONTINUE
*
*              Apply G from the right to transform the columns of the
*              matrix in rows I1 to min(K+3,I).
*
               DO 110 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
  110          CONTINUE
*
               IF( wantz ) THEN
*
*                 Accumulate transformations in the matrix Z
*
                  DO 120 j = iloz, ihiz
                     sum = z( j, k ) + v2*z( j, k+1 )
                     z( j, k ) = z( j, k ) - sum*t1
                     z( j, k+1 ) = z( j, k+1 ) - sum*t2
  120             CONTINUE
               END IF
            END IF
  130    CONTINUE
*
  140 CONTINUE
*
*     Failure to converge in remaining number of iterations
*
      info = i
      RETURN
*
  150 CONTINUE
*
      IF( l.EQ.i ) THEN
*
*        H(I,I-1) is negligible: one eigenvalue has converged.
*
         wr( i ) = h( i, i )
         wi( i ) = zero
      ELSE IF( l.EQ.i-1 ) THEN
*
*        H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
*
*        Transform the 2-by-2 submatrix to standard Schur form,
*        and compute and store the eigenvalues.
*
         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
*
*           Apply the transformation to the rest of H.
*
            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 )
         END IF
         IF( wantz ) THEN
*
*           Apply the transformation to Z.
*
            CALL drot( nz, z( iloz, i-1 ), 1, z( iloz, i ), 1, cs, sn )
         END IF
      END IF
*     reset deflation counter
      kdefl = 0
*
*     return to start of the main loop with new value of I.
*
      i = l - 1
      GO TO 20
*
  160 CONTINUE
      RETURN
*
*     End of DLAHQR
*
      END