dlahqr.f

Go to the documentation of this file.

      SUBROUTINE dlahqr( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
     $                   ILOZ, IHIZ, Z, LDZ, INFO )
*
*  -- LAPACK auxiliary routine (version 3.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     June 30, 1999
*
*     .. Scalar Arguments ..
      LOGICAL            WANTT, WANTZ
      INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )
*     ..
*
*  Purpose
*  =======
*
*  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.
*
*  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 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.
*
*  H       (input/output) DOUBLE PRECISION array, dimension (LDH,N)
*          On entry, the upper Hessenberg matrix H.
*          On exit, 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 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).
*
*  WR      (output) DOUBLE PRECISION array, dimension (N)
*  WI      (output) 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).
*
*  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) 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.
*
*  LDZ     (input) INTEGER
*          The leading dimension of the array Z. LDZ >= max(1,N).
*
*  INFO    (output) INTEGER
*          = 0: successful exit
*          > 0: DLAHQR failed to compute all the eigenvalues ILO to IHI
*               in a total of 30*(IHI-ILO+1) iterations; if INFO = i,
*               elements i+1:ihi of WR and WI contain those eigenvalues
*               which have been successfully computed.
*
*  Further Details
*  ===============
*
*  2-96 Based on modifications by
*     David Day, Sandia National Laboratory, USA
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, HALF
      parameter( zero = 0.0d+0, one = 1.0d+0, half = 0.5d0 )
      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   AVE, CS, DISC, H00, H10, H11, H12, H21, H22,
     $                   h33, h33s, h43h34, h44, h44s, ovfl, s, smlnum,
     $                   sn, sum, t1, t2, t3, tst1, ulp, unfl, v1, v2,
     $                   v3
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION   V( 3 ), WORK( 1 )
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH, DLANHS
      EXTERNAL           dlamch, dlanhs
*     ..
*     .. External Subroutines ..
      EXTERNAL           dcopy, dlanv2, dlarfg, drot
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min, sign, 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
*
      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 = one / 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 = abs( h( k-1, k-1 ) ) + abs( h( k, k ) )
            IF( tst1.EQ.zero )
     $         tst1 = dlanhs( '1', i-l+1, h( l, l ), ldh, work )
            IF( abs( 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( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
            h44 = dat1*s + h( i, i )
            h33 = h44
            h43h34 = dat2*s*s
         ELSE
*
*           Prepare to use Francis' double shift
*           (i.e. 2nd degree generalized Rayleigh quotient)
*
            h44 = h( i, i )
            h33 = h( i-1, i-1 )
            h43h34 = h( i, i-1 )*h( i-1, i )
            s = h( i-1, i-2 )*h( i-1, i-2 )
            disc = ( h33-h44 )*half
            disc = disc*disc + h43h34
            IF( disc.GT.zero ) THEN
*
*              Real roots: use Wilkinson's shift twice
*
               disc = sqrt( disc )
               ave = half*( h33+h44 )
               IF( abs( h33 )-abs( h44 ).GT.zero ) THEN
                  h33 = h33*h44 - h43h34
                  h44 = h33 / ( sign( disc, ave )+ave )
               ELSE
                  h44 = sign( disc, ave ) + ave
               END IF
               h33 = h44
               h43h34 = zero
            END IF
         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 = abs( v1 ) + abs( 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 = abs( v1 )*( abs( h00 )+abs( h11 )+abs( h22 ) )
            IF( abs( h10 )*( abs( v2 )+abs( 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 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
               h( k, k-1 ) = -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 = 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
   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 = 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
   70          CONTINUE
*
               IF( wantz ) THEN
*
*                 Accumulate transformations in the matrix Z
*
                  DO 80 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
   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 = 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
   90          CONTINUE
*
*              Apply G from the right to transform the columns of the
*              matrix in rows I1 to min(K+3,I).
*
               DO 100 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
  100          CONTINUE
*
               IF( wantz ) THEN
*
*                 Accumulate transformations in the matrix Z
*
                  DO 110 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
  110             CONTINUE
               END IF
            END IF
  120    CONTINUE
*
  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.
*
         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
*
*     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 DLAHQR
*
      END