clahqr()

subroutine clahqr	(	logical	wantt,
		logical	wantz,
		integer	n,
		integer	ilo,
		integer	ihi,
		complex, dimension( ldh, * )	h,
		integer	ldh,
		complex, dimension( * )	w,
		integer	iloz,
		integer	ihiz,
		complex, dimension( ldz, * )	z,
		integer	ldz,
		integer	info )

CLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.

Download CLAHQR + dependencies [TGZ] [ZIP] [TXT]

Purpose:

!>
!>    CLAHQR 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.
!> 

Parameters

[in]	WANTT	!> WANTT is LOGICAL !> = .TRUE. : the full Schur form T is required; !> = .FALSE.: only eigenvalues are required. !>
[in]	WANTZ	!> WANTZ is LOGICAL !> = .TRUE. : the matrix of Schur vectors Z is required; !> = .FALSE.: Schur vectors are not required. !>
[in]	N	!> N is INTEGER !> The order of the matrix H. N >= 0. !>
[in]	ILO	!> ILO is INTEGER !>
[in]	IHI	!> 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). !> 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. !>
[in,out]	H	!> H is COMPLEX 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. !>
[in]	LDH	!> LDH is INTEGER !> The leading dimension of the array H. LDH >= max(1,N). !>
[out]	W	!> W is 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). !>
[in]	ILOZ	!> ILOZ is INTEGER !>
[in]	IHIZ	!> IHIZ is INTEGER !> Specify the rows of Z to which transformations must be !> applied if WANTZ is .TRUE.. !> 1 <= ILOZ <= ILO; IHI <= IHIZ <= N. !>
[in,out]	Z	!> Z is COMPLEX 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. !>
[in]	LDZ	!> LDZ is INTEGER !> The leading dimension of the array Z. LDZ >= max(1,N). !>
[out]	INFO	!> INFO is INTEGER !> = 0: successful exit !> > 0: if INFO = i, CLAHQR 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.) !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors:

!>
!>     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 CLAHQR 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).
!> 

Definition at line 191 of file clahqr.f.

      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            H( LDH, * ), W( * ), Z( LDZ, * )
*     ..
*
*  =========================================================
*
*     .. Parameters ..
      COMPLEX            ZERO, ONE
      parameter( zero = ( 0.0e0, 0.0e0 ),
     $                   one = ( 1.0e0, 0.0e0 ) )
      REAL               RZERO, RONE, HALF
      parameter( rzero = 0.0e0, rone = 1.0e0, half = 0.5e0 )
      REAL               DAT1
      parameter( dat1 = 3.0e0 / 4.0e0 )
      INTEGER            KEXSH
      parameter( kexsh = 10 )
*     ..
*     .. Local Scalars ..
      COMPLEX            CDUM, H11, H11S, H22, SC, SUM, T, T1, TEMP, U,
     $                   V2, X, Y
      REAL               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            V( 2 )
*     ..
*     .. External Functions ..
      COMPLEX            CLADIV
      REAL               SLAMCH
      EXTERNAL           cladiv, slamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           ccopy, clarfg, cscal
*     ..
*     .. Statement Functions ..
      REAL               CABS1
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, aimag, conjg, max, min, real, sqrt
*     ..
*     .. 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
*
*     ==== 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( aimag( 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 = conjg( sc ) / abs( sc )
            h( i, i-1 ) = abs( h( i, i-1 ) )
            CALL cscal( jhi-i+1, sc, h( i, i ), ldh )
            CALL cscal( min( jhi, i+1 )-jlo+1, conjg( sc ), h( jlo,
     $                  i ),
     $                  1 )
            IF( wantz )
     $         CALL cscal( ihiz-iloz+1, conjg( 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 = slamch( 'SAFE MINIMUM' )
      safmax = rone / safmin
      ulp = slamch( 'PRECISION' )
      smlnum = safmin*( real( 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( real( h( k-1, k-2 ) ) )
               IF( k+1.LE.ihi )
     $            tst = tst + abs( real( 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( real( 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( real( h( i, i-1 ) ) )
            t = s + h( i, i )
         ELSE IF( mod(kdefl,kexsh).EQ.0 ) THEN
*
*           Exceptional shift.
*
            s = dat1*abs( real( 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( real( x / sx )*real( y )+aimag( x / sx )*
     $                aimag( y ).LT.rzero )y = -y
               END IF
               t = t - u*cladiv( 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 = real( h( m+1, m ) )
            s = cabs1( h11s ) + abs( h21 )
            h11s = h11s / s
            h21 = h21 / s
            v( 1 ) = h11s
            v( 2 ) = h21
            h10 = real( 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 = real( 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 CLARFG, and hence
*           after the call T2 ( = T1*V(2) ) is also real.
*
            IF( k.GT.m )
     $         CALL ccopy( 2, h( k, k-1 ), 1, v, 1 )
            CALL clarfg( 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 = real( 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 = conjg( 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*conjg( 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*conjg( 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 )*conjg( 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 cscal( i2-j, temp, h( j, j+1 ), ldh )
                     CALL cscal( j-i1, conjg( temp ), h( i1, j ), 1 )
                     IF( wantz ) THEN
                        CALL cscal( nz, conjg( 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( aimag( temp ).NE.rzero ) THEN
            rtemp = abs( temp )
            h( i, i-1 ) = rtemp
            temp = temp / rtemp
            IF( i2.GT.i )
     $         CALL cscal( i2-i, conjg( temp ), h( i, i+1 ), ldh )
            CALL cscal( i-i1, temp, h( i1, i ), 1 )
            IF( wantz ) THEN
               CALL cscal( 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 CLAHQR
*

Here is the call graph for this function:

Here is the caller graph for this function: