zlahqr2.f

Go to the documentation of this file.

      SUBROUTINE zlahqr2( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
     $                    IHIZ, Z, LDZ, INFO )
*
*  -- ScaLAPACK routine (version 1.7) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     June 22, 2000
*
*     .. Scalar Arguments ..
      LOGICAL            WANTT, WANTZ
      INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
*     ..
*     .. Array Arguments ..
      COMPLEX*16         H( LDH, * ), W( * ), Z( LDZ, * )
*     ..
*
*  Purpose
*  =======
*
*  ZLAHQR2 is an auxiliary routine called by ZHSEQR to update the
*    eigenvalues and Schur decomposition already computed by ZHSEQR, by
*    dealing with the Hessenberg submatrix in rows and columns ILO to IHI.
*  This version of ZLAHQR (not the standard LAPACK version) uses a
*    double-shift algorithm (like LAPACK's DLAHQR).
*  Unlike the standard LAPACK convention, this does not assume the
*    subdiagonal is real, nor does it work to preserve this quality if
*    given.
*
*  Arguments
*  =========
*
*  WANTT   (input) LOGICAL
*          = .TRUE. : the full Schur form T is required;
*          = .FALSE.: only eigenvalues are required.
*
*  WANTZ   (input) LOGICAL
*          = .TRUE. : the matrix of Schur vectors Z is required;
*          = .FALSE.: Schur vectors are not required.
*
*  N       (input) INTEGER
*          The order of the matrix H.  N >= 0.
*
*  ILO     (input) INTEGER
*  IHI     (input) 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.
*
*  H       (input/output) COMPLEX*16 array, dimension (LDH,N)
*          On entry, the upper Hessenberg matrix H.
*          On exit, if WANTT is .TRUE., H is upper triangular in rows
*          and columns ILO:IHI.  If WANTT is .FALSE., the contents of H
*          are unspecified on exit.
*
*  LDH     (input) INTEGER
*          The leading dimension of the array H. LDH >= max(1,N).
*
*  W       (output) 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).
*
*  ILOZ    (input) INTEGER
*  IHIZ    (input) INTEGER
*          Specify the rows of Z to which transformations must be
*          applied if WANTZ is .TRUE..
*          1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
*
*  Z       (input/output) COMPLEX*16 array, dimension (LDZ,N)
*          If WANTZ is .TRUE., on entry Z must contain the current
*          matrix Z of transformations, 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.
*
*  LDZ     (input) INTEGER
*          The leading dimension of the array Z. LDZ >= max(1,N).
*
*  INFO    (output) INTEGER
*          = 0: successful exit
*          > 0: if INFO = i, ZLAHQR failed to compute all the
*               eigenvalues ILO to IHI in a total of 30*(IHI-ILO+1)
*               iterations; elements i+1:ihi of W contain those
*               eigenvalues which have been successfully computed.
*
*  Further Details
*  ===============
*
*  Modified by Mark R. Fahey, June, 2000
*
*  =====================================================================
*
*     .. Parameters ..
      COMPLEX*16         ZERO
      parameter( zero = ( 0.0d+0, 0.0d+0 ) )
      DOUBLE PRECISION   RZERO, RONE
      parameter( rzero = 0.0d+0, rone = 1.0d+0 )
      DOUBLE PRECISION   DAT1, DAT2
      parameter( dat1 = 0.75d+0, dat2 = -0.4375d+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, I1, I2, ITN, ITS, J, K, L, M, NH, NR, NZ
      DOUBLE PRECISION   CS, OVFL, S, SMLNUM, TST1, ULP, UNFL
      COMPLEX*16         CDUM, H00, H10, H11, H12, H21, H22, H33, H33S,
     $                   h43h34, h44, h44s, sn, sum, t1, t2, t3, v1, v2,
     $                   v3
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION   RWORK( 1 )
      COMPLEX*16         V( 3 )
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH, ZLANHS
      EXTERNAL           dlamch, zlanhs
*     ..
*     .. External Subroutines ..
      EXTERNAL           dlabad, zcopy, zlanv2, zlarfg, zrot
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, dconjg, dimag, max, min
*     ..
*     .. Statement Functions ..
      DOUBLE PRECISION   CABS1
*     ..
*     .. 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
*
      nh = ihi - ilo + 1
      nz = ihiz - iloz + 1
*
*     Set machine-dependent constants for the stopping criterion.
*     If norm(H) <= sqrt(OVFL), overflow should not occur.
*
      unfl = dlamch( 'Safe minimum' )
      ovfl = rone / unfl
      CALL dlabad( unfl, ovfl )
      ulp = dlamch( 'Precision' )
      smlnum = unfl*( 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
*
*     ITN is the total number of QR iterations allowed.
*
      itn = 30*nh
*
*     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
   10 CONTINUE
      l = ilo
      IF( i.LT.ilo )
     $   GO TO 150
*
*     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 130 its = 0, itn
*
*        Look for a single small subdiagonal element.
*
         DO 20 k = i, l + 1, -1
            tst1 = cabs1( h( k-1, k-1 ) ) + cabs1( h( k, k ) )
            IF( tst1.EQ.rzero )
     $         tst1 = zlanhs( '1', i-l+1, h( l, l ), ldh, rwork )
            IF( cabs1( h( k, k-1 ) ).LE.max( ulp*tst1, smlnum ) )
     $         GO TO 30
   20    CONTINUE
   30    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 140
*
*        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( its.EQ.10 .OR. its.EQ.20 ) THEN
*
*           Exceptional shift.
*
*            S = ABS( DBLE( H( I,I-1 ) ) ) + ABS( DBLE( H( I-1,I-2 ) ) )
            s = cabs1( h( i, i-1 ) ) + cabs1( h( i-1, i-2 ) )
            h44 = dat1*s
            h33 = h44
            h43h34 = dat2*s*s
         ELSE
*
*           Prepare to use Wilkinson's shift.
*
            h44 = h( i, i )
            h33 = h( i-1, i-1 )
            h43h34 = h( i, i-1 )*h( i-1, i )
         END IF
*
*        Look for two consecutive small subdiagonal elements.
*
         DO 40 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.
*
            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 = cabs1( v1 ) + cabs1( 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 = cabs1( v1 )*( cabs1( h00 )+cabs1( h11 )+
     $             cabs1( h22 ) )
            IF( cabs1( h10 )*( cabs1( v2 )+cabs1( v3 ) ).LE.ulp*tst1 )
     $         GO TO 50
   40    CONTINUE
   50    CONTINUE
*
*        Double-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.  NR is the order of G
*
            nr = min( 3, i-k+1 )
            IF( k.GT.m )
     $         CALL zcopy( nr, h( k, k-1 ), 1, v, 1 )
            CALL zlarfg( 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
*              The real double-shift code uses H( K, K-1 ) = -H( K, K-1 )
*              instead of the following.
               h( k, k-1 ) = h( k, k-1 ) - dconjg( t1 )*h( k, k-1 )
            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 60 j = k, i2
                  sum = dconjg( t1 )*h( k, j ) +
     $                  dconjg( t2 )*h( k+1, j ) +
     $                  dconjg( t3 )*h( k+2, j )
                  h( k, j ) = h( k, j ) - sum
                  h( k+1, j ) = h( k+1, j ) - sum*v2
                  h( k+2, j ) = h( k+2, j ) - sum*v3
   60          CONTINUE
*
*              Apply G from the right to transform the columns of the
*              matrix in rows I1 to min(K+3,I).
*
               DO 70 j = i1, min( k+3, i )
                  sum = t1*h( j, k ) + t2*h( j, k+1 ) + t3*h( j, k+2 )
                  h( j, k ) = h( j, k ) - sum
                  h( j, k+1 ) = h( j, k+1 ) - sum*dconjg( v2 )
                  h( j, k+2 ) = h( j, k+2 ) - sum*dconjg( v3 )
   70          CONTINUE
*
               IF( wantz ) THEN
*
*              Accumulate transformations in the matrix Z
*
                  DO 80 j = iloz, ihiz
                     sum = t1*z( j, k ) + t2*z( j, k+1 ) +
     $                     t3*z( j, k+2 )
                     z( j, k ) = z( j, k ) - sum
                     z( j, k+1 ) = z( j, k+1 ) - sum*dconjg( v2 )
                     z( j, k+2 ) = z( j, k+2 ) - sum*dconjg( v3 )
   80             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 90 j = k, i2
                  sum = dconjg( t1 )*h( k, j ) +
     $                  dconjg( t2 )*h( k+1, j )
                  h( k, j ) = h( k, j ) - sum
                  h( k+1, j ) = h( k+1, j ) - sum*v2
   90          CONTINUE
*
*              Apply G from the right to transform the columns of the
*              matrix in rows I1 to min(K+2,I).
*
               DO 100 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 )
  100          CONTINUE
*
               IF( wantz ) THEN
*
*                 Accumulate transformations in the matrix Z
*
                  DO 110 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 )
  110             CONTINUE
               END IF
            END IF
*
*           Since at the start of the QR step we have for M > L
*              H( K, K-1 ) = H( K, K-1 ) - DCONJG( T1 )*H( K, K-1 )
*           then we don't need to do the following
*           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 H(M,M-1)
*              must also be updated by a factor of (1-T1).
*              TEMP = ONE - T1
*              H( m, m-1 ) = H( m, m-1 )*DCONJG( TEMP )
*           END IF
  120    CONTINUE
*
*        Ensure that H(I,I-1) is real.
*
  130 CONTINUE
*
*     Failure to converge in remaining number of iterations
*
      info = i
      RETURN
*
  140 CONTINUE
*
      IF( l.EQ.i ) THEN
*
*        H(I,I-1) is negligible: one eigenvalue has converged.
*
         w( i ) = h( i, i )
*
      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 zlanv2( h( i-1, i-1 ), h( i-1, i ), h( i, i-1 ),
     $                h( i, i ), w( i-1 ), w( i ), cs, sn )
*
         IF( wantt ) THEN
*
*           Apply the transformation to the rest of H.
*
            IF( i2.GT.i )
     $         CALL zrot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,
     $                    cs, sn )
            CALL zrot( i-i1-1, h( i1, i-1 ), 1, h( i1, i ), 1, cs,
     $                 dconjg( sn ) )
         END IF
         IF( wantz ) THEN
*
*           Apply the transformation to Z.
*
            CALL zrot( nz, z( iloz, i-1 ), 1, z( iloz, i ), 1, cs,
     $                 dconjg( sn ) )
         END IF
*
      END IF
*
*     Decrement number of remaining iterations, and return to start of
*     the main loop with new value of I.
*
      itn = itn - its
      i = l - 1
      GO TO 10
*
  150 CONTINUE
      RETURN
*
*     End of ZLAHQR2
*
      END