dhseqr.f

Go to the documentation of this file.

      SUBROUTINE dhseqr( JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z,
     $                   LDZ, WORK, LWORK, INFO )
*
*  -- LAPACK 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 ..
      CHARACTER          COMPZ, JOB
      INTEGER            IHI, ILO, INFO, LDH, LDZ, LWORK, N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   H( LDH, * ), WI( * ), WORK( * ), WR( * ),
     $                   z( ldz, * )
*     ..
*
*  Purpose
*  =======
*
*  DHSEQR computes the eigenvalues of a real upper Hessenberg matrix H
*  and, optionally, the matrices T and Z from the Schur decomposition
*  H = Z T Z**T, where T is an upper quasi-triangular matrix (the Schur
*  form), and Z is the orthogonal matrix of Schur vectors.
*
*  Optionally Z may be postmultiplied into an input orthogonal matrix Q,
*  so that this routine can give the Schur factorization of a matrix A
*  which has been reduced to the Hessenberg form H by the orthogonal
*  matrix Q:  A = Q*H*Q**T = (QZ)*T*(QZ)**T.
*
*  Arguments
*  =========
*
*  JOB     (input) CHARACTER*1
*          = 'E':  compute eigenvalues only;
*          = 'S':  compute eigenvalues and the Schur form T.
*
*  COMPZ   (input) CHARACTER*1
*          = 'N':  no Schur vectors are computed;
*          = 'I':  Z is initialized to the unit matrix and the matrix Z
*                  of Schur vectors of H is returned;
*          = 'V':  Z must contain an orthogonal matrix Q on entry, and
*                  the product Q*Z is returned.
*
*  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 1:ILO-1 and IHI+1:N. ILO and IHI are normally
*          set by a previous call to DGEBAL, and then passed to SGEHRD
*          when the matrix output by DGEBAL is reduced to Hessenberg
*          form. Otherwise ILO and IHI should be set to 1 and N
*          respectively.
*          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
*
*  H       (input/output) DOUBLE PRECISION array, dimension (LDH,N)
*          On entry, the upper Hessenberg matrix H.
*          On exit, if JOB = 'S', H contains the upper quasi-triangular
*          matrix T from the Schur decomposition (the Schur form);
*          2-by-2 diagonal blocks (corresponding to complex conjugate
*          pairs of eigenvalues) are returned in standard form, with
*          H(i,i) = H(i+1,i+1) and H(i+1,i)*H(i,i+1) < 0. If JOB = 'E',
*          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. 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 JOB = 'S', 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).
*
*  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
*          If COMPZ = 'N': Z is not referenced.
*          If COMPZ = 'I': on entry, Z need not be set, and on exit, Z
*          contains the orthogonal matrix Z of the Schur vectors of H.
*          If COMPZ = 'V': on entry Z must contain an N-by-N matrix Q,
*          which is assumed to be equal to the unit matrix except for
*          the submatrix Z(ILO:IHI,ILO:IHI); on exit Z contains Q*Z.
*          Normally Q is the orthogonal matrix generated by DORGHR after
*          the call to DGEHRD which formed the Hessenberg matrix H.
*
*  LDZ     (input) INTEGER
*          The leading dimension of the array Z.
*          LDZ >= max(1,N) if COMPZ = 'I' or 'V'; LDZ >= 1 otherwise.
*
*  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
*  LWORK   (input) INTEGER
*          The dimension of the array WORK.  LWORK >= max(1,N).
*
*          If LWORK = -1, then a workspace query is assumed; the routine
*          only calculates the optimal size of the WORK array, returns
*          this value as the first entry of the WORK array, and no error
*          message related to LWORK is issued by XERBLA.
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value
*          > 0:  if INFO = i, DHSEQR failed to compute all of the
*                eigenvalues in a total of 30*(IHI-ILO+1) iterations;
*                elements 1:ilo-1 and i+1:n of WR and WI contain those
*                eigenvalues which have been successfully computed.
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, TWO
      parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
      DOUBLE PRECISION   CONST
      parameter( const = 1.5d+0 )
      INTEGER            NSMAX, LDS
      parameter( nsmax = 15, lds = nsmax )
*     ..
*     .. Local Scalars ..
      LOGICAL            INITZ, LQUERY, WANTT, WANTZ
      INTEGER            I, I1, I2, IERR, II, ITEMP, ITN, ITS, J, K, L,
     $                   maxb, nh, nr, ns, nv
      DOUBLE PRECISION   ABSW, OVFL, SMLNUM, TAU, TEMP, TST1, ULP, UNFL
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION   S( LDS, NSMAX ), V( NSMAX+1 ), VV( NSMAX+1 )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            IDAMAX, ILAENV
      DOUBLE PRECISION   DLAMCH, DLANHS, DLAPY2
      EXTERNAL           lsame, idamax, ilaenv, dlamch, dlanhs, dlapy2
*     ..
*     .. External Subroutines ..
      EXTERNAL           dcopy, dgemv, dlacpy, dlahqr, dlarfg, dlarfx,
     $                   dlaset, dscal, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min
*     ..
*     .. Executable Statements ..
*
*     Decode and test the input parameters
*
      wantt = lsame( job, 'S' )
      initz = lsame( compz, 'I' )
      wantz = initz .OR. lsame( compz, 'V' )
*
      info = 0
      work( 1 ) = max( 1, n )
      lquery = ( lwork.EQ.-1 )
      IF( .NOT.lsame( job, 'E' ) .AND. .NOT.wantt ) THEN
         info = -1
      ELSE IF( .NOT.lsame( compz, 'N' ) .AND. .NOT.wantz ) THEN
         info = -2
      ELSE IF( n.LT.0 ) THEN
         info = -3
      ELSE IF( ilo.LT.1 .OR. ilo.GT.max( 1, n ) ) THEN
         info = -4
      ELSE IF( ihi.LT.min( ilo, n ) .OR. ihi.GT.n ) THEN
         info = -5
      ELSE IF( ldh.LT.max( 1, n ) ) THEN
         info = -7
      ELSE IF( ldz.LT.1 .OR. wantz .AND. ldz.LT.max( 1, n ) ) THEN
         info = -11
      ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
         info = -13
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'DHSEQR', -info )
         RETURN
      ELSE IF( lquery ) THEN
         RETURN
      END IF
*
*     Initialize Z, if necessary
*
      IF( initz )
     $   CALL dlaset( 'Full', n, n, zero, one, z, ldz )
*
*     Store the eigenvalues isolated by DGEBAL.
*
      DO 10 i = 1, ilo - 1
         wr( i ) = h( i, i )
         wi( i ) = zero
   10 CONTINUE
      DO 20 i = ihi + 1, n
         wr( i ) = h( i, i )
         wi( i ) = zero
   20 CONTINUE
*
*     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
*
*     Set rows and columns ILO to IHI to zero below the first
*     subdiagonal.
*
      DO 40 j = ilo, ihi - 2
         DO 30 i = j + 2, n
            h( i, j ) = zero
   30    CONTINUE
   40 CONTINUE
      nh = ihi - ilo + 1
*
*     Determine the order of the multi-shift QR algorithm to be used.
*
      ns = ilaenv( 4, 'DHSEQR', job // compz, n, ilo, ihi, -1 )
      maxb = ilaenv( 8, 'DHSEQR', job // compz, n, ilo, ihi, -1 )
      IF( ns.LE.2 .OR. ns.GT.nh .OR. maxb.GE.nh ) THEN
*
*        Use the standard double-shift algorithm
*
         CALL dlahqr( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,
     $                ihi, z, ldz, info )
         RETURN
      END IF
      maxb = max( 3, maxb )
      ns = min( ns, maxb, nsmax )
*
*     Now 2 < NS <= MAXB < NH.
*
*     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 multiple-shift 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 at most MAXB. 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
   50 CONTINUE
      l = ilo
      IF( i.LT.ilo )
     $   GO TO 170
*
*     Perform multiple-shift QR iterations on rows and columns ILO to I
*     until a submatrix of order at most MAXB splits off at the bottom
*     because a subdiagonal element has become negligible.
*
      DO 150 its = 0, itn
*
*        Look for a single small subdiagonal element.
*
         DO 60 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 70
   60    CONTINUE
   70    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 <= MAXB has split off.
*
         IF( l.GE.i-maxb+1 )
     $      GO TO 160
*
*        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.20 .OR. its.EQ.30 ) THEN
*
*           Exceptional shifts.
*
            DO 80 ii = i - ns + 1, i
               wr( ii ) = const*( abs( h( ii, ii-1 ) )+
     $                    abs( h( ii, ii ) ) )
               wi( ii ) = zero
   80       CONTINUE
         ELSE
*
*           Use eigenvalues of trailing submatrix of order NS as shifts.
*
            CALL dlacpy( 'Full', ns, ns, h( i-ns+1, i-ns+1 ), ldh, s,
     $                   lds )
            CALL dlahqr( .false., .false., ns, 1, ns, s, lds,
     $                   wr( i-ns+1 ), wi( i-ns+1 ), 1, ns, z, ldz,
     $                   ierr )
            IF( ierr.GT.0 ) THEN
*
*              If DLAHQR failed to compute all NS eigenvalues, use the
*              unconverged diagonal elements as the remaining shifts.
*
               DO 90 ii = 1, ierr
                  wr( i-ns+ii ) = s( ii, ii )
                  wi( i-ns+ii ) = zero
   90          CONTINUE
            END IF
         END IF
*
*        Form the first column of (G-w(1)) (G-w(2)) . . . (G-w(ns))
*        where G is the Hessenberg submatrix H(L:I,L:I) and w is
*        the vector of shifts (stored in WR and WI). The result is
*        stored in the local array V.
*
         v( 1 ) = one
         DO 100 ii = 2, ns + 1
            v( ii ) = zero
  100    CONTINUE
         nv = 1
         DO 120 j = i - ns + 1, i
            IF( wi( j ).GE.zero ) THEN
               IF( wi( j ).EQ.zero ) THEN
*
*                 real shift
*
                  CALL dcopy( nv+1, v, 1, vv, 1 )
                  CALL dgemv( 'No transpose', nv+1, nv, one, h( l, l ),
     $                        ldh, vv, 1, -wr( j ), v, 1 )
                  nv = nv + 1
               ELSE IF( wi( j ).GT.zero ) THEN
*
*                 complex conjugate pair of shifts
*
                  CALL dcopy( nv+1, v, 1, vv, 1 )
                  CALL dgemv( 'No transpose', nv+1, nv, one, h( l, l ),
     $                        ldh, v, 1, -two*wr( j ), vv, 1 )
                  itemp = idamax( nv+1, vv, 1 )
                  temp = one / max( abs( vv( itemp ) ), smlnum )
                  CALL dscal( nv+1, temp, vv, 1 )
                  absw = dlapy2( wr( j ), wi( j ) )
                  temp = ( temp*absw )*absw
                  CALL dgemv( 'No transpose', nv+2, nv+1, one,
     $                        h( l, l ), ldh, vv, 1, temp, v, 1 )
                  nv = nv + 2
               END IF
*
*              Scale V(1:NV) so that max(abs(V(i))) = 1. If V is zero,
*              reset it to the unit vector.
*
               itemp = idamax( nv, v, 1 )
               temp = abs( v( itemp ) )
               IF( temp.EQ.zero ) THEN
                  v( 1 ) = one
                  DO 110 ii = 2, nv
                     v( ii ) = zero
  110             CONTINUE
               ELSE
                  temp = max( temp, smlnum )
                  CALL dscal( nv, one / temp, v, 1 )
               END IF
            END IF
  120    CONTINUE
*
*        Multiple-shift QR step
*
         DO 140 k = l, 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( ns+1, i-k+1 )
            IF( k.GT.l )
     $         CALL dcopy( nr, h( k, k-1 ), 1, v, 1 )
            CALL dlarfg( nr, v( 1 ), v( 2 ), 1, tau )
            IF( k.GT.l ) THEN
               h( k, k-1 ) = v( 1 )
               DO 130 ii = k + 1, i
                  h( ii, k-1 ) = zero
  130          CONTINUE
            END IF
            v( 1 ) = one
*
*           Apply G from the left to transform the rows of the matrix in
*           columns K to I2.
*
            CALL dlarfx( 'Left', nr, i2-k+1, v, tau, h( k, k ), ldh,
     $                   work )
*
*           Apply G from the right to transform the columns of the
*           matrix in rows I1 to min(K+NR,I).
*
            CALL dlarfx( 'Right', min( k+nr, i )-i1+1, nr, v, tau,
     $                   h( i1, k ), ldh, work )
*
            IF( wantz ) THEN
*
*              Accumulate transformations in the matrix Z
*
               CALL dlarfx( 'Right', nh, nr, v, tau, z( ilo, k ), ldz,
     $                      work )
            END IF
  140    CONTINUE
*
  150 CONTINUE
*
*     Failure to converge in remaining number of iterations
*
      info = i
      RETURN
*
  160 CONTINUE
*
*     A submatrix of order <= MAXB in rows and columns L to I has split
*     off. Use the double-shift QR algorithm to handle it.
*
      CALL dlahqr( wantt, wantz, n, l, i, h, ldh, wr, wi, ilo, ihi, z,
     $             ldz, info )
      IF( info.GT.0 )
     $   RETURN
*
*     Decrement number of remaining iterations, and return to start of
*     the main loop with a new value of I.
*
      itn = itn - its
      i = l - 1
      GO TO 50
*
  170 CONTINUE
      work( 1 ) = max( 1, n )
      RETURN
*
*     End of DHSEQR
*
      END