zlahqr.f

*> \brief \b ZLAHQR 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 ZLAHQR + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlahqr.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlahqr.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlahqr.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
*                          IHIZ, Z, LDZ, INFO )
*
*       .. Scalar Arguments ..
*       INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
*       LOGICAL            WANTT, WANTZ
*       ..
*       .. Array Arguments ..
*       COMPLEX*16         H( LDH, * ), W( * ), Z( LDZ, * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*>    ZLAHQR is an auxiliary routine called by CHSEQR to update the
*>    eigenvalues and Schur decomposition already computed by CHSEQR, 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 triangular in rows and
*>          columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1).
*>          ZLAHQR 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 COMPLEX*16 array, dimension (LDH,N)
*>          On entry, the upper Hessenberg matrix H.
*>          On exit, if INFO is zero and if WANTT is .TRUE., then H
*>          is upper triangular in rows and columns ILO:IHI.  If INFO
*>          is zero and if WANTT is .FALSE., then the contents of H
*>          are unspecified on exit.  The output state of H in case
*>          INF is positive is 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] W
*> \verbatim
*>          W is COMPLEX*16 array, dimension (N)
*>          The computed eigenvalues ILO to IHI are stored in the
*>          corresponding elements of W. If WANTT is .TRUE., the
*>          eigenvalues are stored in the same order as on the diagonal
*>          of the Schur form returned in H, with W(i) = H(i,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 COMPLEX*16 array, dimension (LDZ,N)
*>          If WANTZ is .TRUE., on entry Z must contain the current
*>          matrix Z of transformations accumulated by CHSEQR, 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, ZLAHQR failed to compute all the
*>                  eigenvalues ILO to IHI in a total of 30 iterations
*>                  per eigenvalue; elements i+1:ihi of W 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 Contributors:
*  ==================
*>
*> \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 ZLAHQR 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 zlahqr( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, 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 ..
      COMPLEX*16         H( ldh, * ), W( * ), Z( ldz, * )
*     ..
*
*  =========================================================
*
*     .. Parameters ..
      COMPLEX*16         ZERO, ONE
      parameter( zero = ( 0.0d0, 0.0d0 ),
     $                   one = ( 1.0d0, 0.0d0 ) )
      DOUBLE PRECISION   RZERO, RONE, HALF
      parameter( rzero = 0.0d0, rone = 1.0d0, half = 0.5d0 )
      DOUBLE PRECISION   DAT1
      parameter( dat1 = 3.0d0 / 4.0d0 )
      INTEGER            KEXSH
      parameter( kexsh = 10 )
*     ..
*     .. Local Scalars ..
      COMPLEX*16         CDUM, H11, H11S, H22, SC, SUM, T, T1, TEMP, U,
     $                   v2, x, y
      DOUBLE PRECISION   AA, AB, BA, BB, H10, H21, RTEMP, S, SAFMAX,
     $                   safmin, smlnum, sx, t2, tst, ulp
      INTEGER            I, I1, I2, ITS, ITMAX, J, JHI, JLO, K, L, M,
     $                   nh, nz, kdefl
*     ..
*     .. Local Arrays ..
      COMPLEX*16         V( 2 )
*     ..
*     .. External Functions ..
      COMPLEX*16         ZLADIV
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           zladiv, dlamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           dlabad, zcopy, zlarfg, zscal
*     ..
*     .. Statement Functions ..
      DOUBLE PRECISION   CABS1
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, dconjg, dimag, max, min, sqrt
*     ..
*     .. Statement Function definitions ..
      cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
*     ..
*     .. Executable Statements ..
*
      info = 0
*
*     Quick return if possible
*
      IF( n.EQ.0 )
     $   RETURN
      IF( ilo.EQ.ihi ) THEN
         w( ilo ) = h( ilo, ilo )
         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
*     ==== ensure that subdiagonal entries are real ====
      IF( wantt ) THEN
         jlo = 1
         jhi = n
      ELSE
         jlo = ilo
         jhi = ihi
      END IF
      DO 20 i = ilo + 1, ihi
         IF( dimag( h( i, i-1 ) ).NE.rzero ) THEN
*           ==== The following redundant normalization
*           .    avoids problems with both gradual and
*           .    sudden underflow in ABS(H(I,I-1)) ====
            sc = h( i, i-1 ) / cabs1( h( i, i-1 ) )
            sc = dconjg( sc ) / abs( sc )
            h( i, i-1 ) = abs( h( i, i-1 ) )
            CALL zscal( jhi-i+1, sc, h( i, i ), ldh )
            CALL zscal( min( jhi, i+1 )-jlo+1, dconjg( sc ),
     $                  h( jlo, i ), 1 )
            IF( wantz )
     $         CALL zscal( ihiz-iloz+1, dconjg( sc ), z( iloz, i ), 1 )
         END IF
   20 CONTINUE
*
      nh = ihi - ilo + 1
      nz = ihiz - iloz + 1
*
*     Set machine-dependent constants for the stopping criterion.
*
      safmin = dlamch( 'SAFE MINIMUM' )
      safmax = rone / 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. 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
   30 CONTINUE
      IF( i.LT.ilo )
     $   GO TO 150
*
*     Perform QR iterations on rows and columns ILO to I until a
*     submatrix of order 1 splits off at the bottom because a
*     subdiagonal element has become negligible.
*
      l = ilo
      DO 130 its = 0, itmax
*
*        Look for a single small subdiagonal element.
*
         DO 40 k = i, l + 1, -1
            IF( cabs1( h( k, k-1 ) ).LE.smlnum )
     $         GO TO 50
            tst = cabs1( h( k-1, k-1 ) ) + cabs1( h( k, k ) )
            IF( tst.EQ.zero ) THEN
               IF( k-2.GE.ilo )
     $            tst = tst + abs( dble( h( k-1, k-2 ) ) )
               IF( k+1.LE.ihi )
     $            tst = tst + abs( dble( 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 examples.  ====
            IF( abs( dble( h( k, k-1 ) ) ).LE.ulp*tst ) THEN
               ab = max( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
               ba = min( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) )
               aa = max( cabs1( h( k, k ) ),
     $              cabs1( h( k-1, k-1 )-h( k, k ) ) )
               bb = min( cabs1( h( k, k ) ),
     $              cabs1( h( k-1, k-1 )-h( k, k ) ) )
               s = aa + ab
               IF( ba*( ab / s ).LE.max( smlnum,
     $             ulp*( bb*( aa / s ) ) ) )GO TO 50
            END IF
   40    CONTINUE
   50    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 has split off.
*
         IF( l.GE.i )
     $      GO TO 140
         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 = dat1*abs( dble( h( i, i-1 ) ) )
            t = s + h( i, i )
         ELSE IF( mod(kdefl,kexsh).EQ.0 ) THEN
*
*           Exceptional shift.
*
            s = dat1*abs( dble( h( l+1, l ) ) )
            t = s + h( l, l )
         ELSE
*
*           Wilkinson's shift.
*
            t = h( i, i )
            u = sqrt( h( i-1, i ) )*sqrt( h( i, i-1 ) )
            s = cabs1( u )
            IF( s.NE.rzero ) THEN
               x = half*( h( i-1, i-1 )-t )
               sx = cabs1( x )
               s = max( s, cabs1( x ) )
               y = s*sqrt( ( x / s )**2+( u / s )**2 )
               IF( sx.GT.rzero ) THEN
                  IF( dble( x / sx )*dble( y )+dimag( x / sx )*
     $                dimag( y ).LT.rzero )y = -y
               END IF
               t = t - u*zladiv( u, ( x+y ) )
            END IF
         END IF
*
*        Look for two consecutive small subdiagonal elements.
*
         DO 60 m = i - 1, l + 1, -1
*
*           Determine the effect of starting the single-shift QR
*           iteration at row M, and see if this would make H(M,M-1)
*           negligible.
*
            h11 = h( m, m )
            h22 = h( m+1, m+1 )
            h11s = h11 - t
            h21 = dble( h( m+1, m ) )
            s = cabs1( h11s ) + abs( h21 )
            h11s = h11s / s
            h21 = h21 / s
            v( 1 ) = h11s
            v( 2 ) = h21
            h10 = dble( h( m, m-1 ) )
            IF( abs( h10 )*abs( h21 ).LE.ulp*
     $          ( cabs1( h11s )*( cabs1( h11 )+cabs1( h22 ) ) ) )
     $          GO TO 70
   60    CONTINUE
         h11 = h( l, l )
         h22 = h( l+1, l+1 )
         h11s = h11 - t
         h21 = dble( h( l+1, l ) )
         s = cabs1( h11s ) + abs( h21 )
         h11s = h11s / s
         h21 = h21 / s
         v( 1 ) = h11s
         v( 2 ) = h21
   70    CONTINUE
*
*        Single-shift QR step
*
         DO 120 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.
*
*           V(2) is always real before the call to ZLARFG, and hence
*           after the call T2 ( = T1*V(2) ) is also real.
*
            IF( k.GT.m )
     $         CALL zcopy( 2, h( k, k-1 ), 1, v, 1 )
            CALL zlarfg( 2, v( 1 ), v( 2 ), 1, t1 )
            IF( k.GT.m ) THEN
               h( k, k-1 ) = v( 1 )
               h( k+1, k-1 ) = zero
            END IF
            v2 = v( 2 )
            t2 = dble( t1*v2 )
*
*           Apply G from the left to transform the rows of the matrix
*           in columns K to I2.
*
            DO 80 j = k, i2
               sum = dconjg( t1 )*h( k, j ) + t2*h( k+1, j )
               h( k, j ) = h( k, j ) - sum
               h( k+1, j ) = h( k+1, j ) - sum*v2
   80       CONTINUE
*
*           Apply G from the right to transform the columns of the
*           matrix in rows I1 to min(K+2,I).
*
            DO 90 j = i1, min( k+2, i )
               sum = t1*h( j, k ) + t2*h( j, k+1 )
               h( j, k ) = h( j, k ) - sum
               h( j, k+1 ) = h( j, k+1 ) - sum*dconjg( v2 )
   90       CONTINUE
*
            IF( wantz ) THEN
*
*              Accumulate transformations in the matrix Z
*
               DO 100 j = iloz, ihiz
                  sum = t1*z( j, k ) + t2*z( j, k+1 )
                  z( j, k ) = z( j, k ) - sum
                  z( j, k+1 ) = z( j, k+1 ) - sum*dconjg( v2 )
  100          CONTINUE
            END IF
*
            IF( k.EQ.m .AND. m.GT.l ) THEN
*
*              If the QR step was started at row M > L because two
*              consecutive small subdiagonals were found, then extra
*              scaling must be performed to ensure that H(M,M-1) remains
*              real.
*
               temp = one - t1
               temp = temp / abs( temp )
               h( m+1, m ) = h( m+1, m )*dconjg( temp )
               IF( m+2.LE.i )
     $            h( m+2, m+1 ) = h( m+2, m+1 )*temp
               DO 110 j = m, i
                  IF( j.NE.m+1 ) THEN
                     IF( i2.GT.j )
     $                  CALL zscal( i2-j, temp, h( j, j+1 ), ldh )
                     CALL zscal( j-i1, dconjg( temp ), h( i1, j ), 1 )
                     IF( wantz ) THEN
                        CALL zscal( nz, dconjg( temp ), z( iloz, j ),
     $                              1 )
                     END IF
                  END IF
  110          CONTINUE
            END IF
  120    CONTINUE
*
*        Ensure that H(I,I-1) is real.
*
         temp = h( i, i-1 )
         IF( dimag( temp ).NE.rzero ) THEN
            rtemp = abs( temp )
            h( i, i-1 ) = rtemp
            temp = temp / rtemp
            IF( i2.GT.i )
     $         CALL zscal( i2-i, dconjg( temp ), h( i, i+1 ), ldh )
            CALL zscal( i-i1, temp, h( i1, i ), 1 )
            IF( wantz ) THEN
               CALL zscal( nz, temp, z( iloz, i ), 1 )
            END IF
         END IF
*
  130 CONTINUE
*
*     Failure to converge in remaining number of iterations
*
      info = i
      RETURN
*
  140 CONTINUE
*
*     H(I,I-1) is negligible: one eigenvalue has converged.
*
      w( i ) = h( i, i )
*     reset deflation counter
      kdefl = 0
*
*     return to start of the main loop with new value of I.
*
      i = l - 1
      GO TO 30
*
  150 CONTINUE
      RETURN
*
*     End of ZLAHQR
*
      END