clahqr2.f

Go to the documentation of this file.

      SUBROUTINE clahqr2( 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            H( LDH, * ), W( * ), Z( LDZ, * )
*     ..
*
*  Purpose
*  =======
*
*  CLAHQR2 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.
*  This version of CLAHQR (not the standard LAPACK version) uses a
*    double-shift algorithm (like LAPACK's SLAHQR).
*  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).
*          CLAHQR 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 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 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 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, CLAHQR 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            ZERO
      parameter( zero = ( 0.0e+0, 0.0e+0 ) )
      REAL               RZERO, RONE
      parameter( rzero = 0.0e+0, rone = 1.0e+0 )
      REAL               DAT1, DAT2
      parameter( dat1 = 0.75e+0, dat2 = -0.4375e+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, I1, I2, ITN, ITS, J, K, L, M, NH, NR, NZ
      REAL               CS, OVFL, S, SMLNUM, TST1, ULP, UNFL
      COMPLEX            CDUM, H00, H10, H11, H12, H21, H22, H33, H33S,
     $                   h43h34, h44, h44s, sn, sum, t1, t2, t3, v1, v2,
     $                   v3
*     ..
*     .. Local Arrays ..
      REAL               RWORK( 1 )
      COMPLEX            V( 3 )
*     ..
*     .. External Functions ..
      REAL               SLAMCH, CLANHS
      EXTERNAL           slamch, clanhs
*     ..
*     .. External Subroutines ..
      EXTERNAL           slabad, ccopy, clanv2, clarfg, crot
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, real, conjg, aimag, max, min
*     ..
*     .. Statement Functions ..
      REAL               CABS1
*     ..
*     .. Statement Function definitions ..
      cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( 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 = slamch( 'Safe minimum' )
      ovfl = rone / unfl
      CALL slabad( unfl, ovfl )
      ulp = slamch( '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 = clanhs( '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( REAL( H( I,I-1 ) ) ) + ABS( REAL( 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 ccopy( nr, h( k, k-1 ), 1, v, 1 )
            CALL clarfg( 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 ) - conjg( 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 = conjg( t1 )*h( k, j ) +
     $                  conjg( t2 )*h( k+1, j ) +
     $                  conjg( 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*conjg( v2 )
                  h( j, k+2 ) = h( j, k+2 ) - sum*conjg( 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*conjg( v2 )
                     z( j, k+2 ) = z( j, k+2 ) - sum*conjg( 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 = conjg( t1 )*h( k, j ) +
     $                  conjg( 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*conjg( 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*conjg( 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 ) - CONJG( 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 )*CONJG( 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 clanv2( 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 crot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,
     $                    cs, sn )
            CALL crot( i-i1-1, h( i1, i-1 ), 1, h( i1, i ), 1, cs,
     $                 conjg( sn ) )
         END IF
         IF( wantz ) THEN
*
*           Apply the transformation to Z.
*
            CALL crot( nz, z( iloz, i-1 ), 1, z( iloz, i ), 1, cs,
     $                 conjg( 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 CLAHQR2
*
      END