LAPACK: SRC/clahqr.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
00002      $                   IHIZ, Z, LDZ, INFO )
00003 *
00004 *  -- LAPACK auxiliary routine (version 3.2) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
00010       LOGICAL            WANTT, WANTZ
00011 *     ..
00012 *     .. Array Arguments ..
00013       COMPLEX            H( LDH, * ), W( * ), Z( LDZ, * )
00014 *     ..
00015 *
00016 *     Purpose
00017 *     =======
00018 *
00019 *     CLAHQR is an auxiliary routine called by CHSEQR to update the
00020 *     eigenvalues and Schur decomposition already computed by CHSEQR, by
00021 *     dealing with the Hessenberg submatrix in rows and columns ILO to
00022 *     IHI.
00023 *
00024 *     Arguments
00025 *     =========
00026 *
00027 *     WANTT   (input) LOGICAL
00028 *          = .TRUE. : the full Schur form T is required;
00029 *          = .FALSE.: only eigenvalues are required.
00030 *
00031 *     WANTZ   (input) LOGICAL
00032 *          = .TRUE. : the matrix of Schur vectors Z is required;
00033 *          = .FALSE.: Schur vectors are not required.
00034 *
00035 *     N       (input) INTEGER
00036 *          The order of the matrix H.  N >= 0.
00037 *
00038 *     ILO     (input) INTEGER
00039 *     IHI     (input) INTEGER
00040 *          It is assumed that H is already upper triangular in rows and
00041 *          columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1).
00042 *          CLAHQR works primarily with the Hessenberg submatrix in rows
00043 *          and columns ILO to IHI, but applies transformations to all of
00044 *          H if WANTT is .TRUE..
00045 *          1 <= ILO <= max(1,IHI); IHI <= N.
00046 *
00047 *     H       (input/output) COMPLEX array, dimension (LDH,N)
00048 *          On entry, the upper Hessenberg matrix H.
00049 *          On exit, if INFO is zero and if WANTT is .TRUE., then H
00050 *          is upper triangular in rows and columns ILO:IHI.  If INFO
00051 *          is zero and if WANTT is .FALSE., then the contents of H
00052 *          are unspecified on exit.  The output state of H in case
00053 *          INF is positive is below under the description of INFO.
00054 *
00055 *     LDH     (input) INTEGER
00056 *          The leading dimension of the array H. LDH >= max(1,N).
00057 *
00058 *     W       (output) COMPLEX array, dimension (N)
00059 *          The computed eigenvalues ILO to IHI are stored in the
00060 *          corresponding elements of W. If WANTT is .TRUE., the
00061 *          eigenvalues are stored in the same order as on the diagonal
00062 *          of the Schur form returned in H, with W(i) = H(i,i).
00063 *
00064 *     ILOZ    (input) INTEGER
00065 *     IHIZ    (input) INTEGER
00066 *          Specify the rows of Z to which transformations must be
00067 *          applied if WANTZ is .TRUE..
00068 *          1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
00069 *
00070 *     Z       (input/output) COMPLEX array, dimension (LDZ,N)
00071 *          If WANTZ is .TRUE., on entry Z must contain the current
00072 *          matrix Z of transformations accumulated by CHSEQR, and on
00073 *          exit Z has been updated; transformations are applied only to
00074 *          the submatrix Z(ILOZ:IHIZ,ILO:IHI).
00075 *          If WANTZ is .FALSE., Z is not referenced.
00076 *
00077 *     LDZ     (input) INTEGER
00078 *          The leading dimension of the array Z. LDZ >= max(1,N).
00079 *
00080 *     INFO    (output) INTEGER
00081 *           =   0: successful exit
00082 *          .GT. 0: if INFO = i, CLAHQR failed to compute all the
00083 *                  eigenvalues ILO to IHI in a total of 30 iterations
00084 *                  per eigenvalue; elements i+1:ihi of W contain
00085 *                  those eigenvalues which have been successfully
00086 *                  computed.
00087 *
00088 *                  If INFO .GT. 0 and WANTT is .FALSE., then on exit,
00089 *                  the remaining unconverged eigenvalues are the
00090 *                  eigenvalues of the upper Hessenberg matrix
00091 *                  rows and columns ILO thorugh INFO of the final,
00092 *                  output value of H.
00093 *
00094 *                  If INFO .GT. 0 and WANTT is .TRUE., then on exit
00095 *          (*)       (initial value of H)*U  = U*(final value of H)
00096 *                  where U is an orthognal matrix.    The final
00097 *                  value of H is upper Hessenberg and triangular in
00098 *                  rows and columns INFO+1 through IHI.
00099 *
00100 *                  If INFO .GT. 0 and WANTZ is .TRUE., then on exit
00101 *                      (final value of Z)  = (initial value of Z)*U
00102 *                  where U is the orthogonal matrix in (*)
00103 *                  (regardless of the value of WANTT.)
00104 *
00105 *     Further Details
00106 *     ===============
00107 *
00108 *     02-96 Based on modifications by
00109 *     David Day, Sandia National Laboratory, USA
00110 *
00111 *     12-04 Further modifications by
00112 *     Ralph Byers, University of Kansas, USA
00113 *     This is a modified version of CLAHQR from LAPACK version 3.0.
00114 *     It is (1) more robust against overflow and underflow and
00115 *     (2) adopts the more conservative Ahues & Tisseur stopping
00116 *     criterion (LAWN 122, 1997).
00117 *
00118 *     =========================================================
00119 *
00120 *     .. Parameters ..
00121       INTEGER            ITMAX
00122       PARAMETER          ( ITMAX = 30 )
00123       COMPLEX            ZERO, ONE
00124       PARAMETER          ( ZERO = ( 0.0e0, 0.0e0 ),
00125      $                   ONE = ( 1.0e0, 0.0e0 ) )
00126       REAL               RZERO, RONE, HALF
00127       PARAMETER          ( RZERO = 0.0e0, RONE = 1.0e0, HALF = 0.5e0 )
00128       REAL               DAT1
00129       PARAMETER          ( DAT1 = 3.0e0 / 4.0e0 )
00130 *     ..
00131 *     .. Local Scalars ..
00132       COMPLEX            CDUM, H11, H11S, H22, SC, SUM, T, T1, TEMP, U,
00133      $                   V2, X, Y
00134       REAL               AA, AB, BA, BB, H10, H21, RTEMP, S, SAFMAX,
00135      $                   SAFMIN, SMLNUM, SX, T2, TST, ULP
00136       INTEGER            I, I1, I2, ITS, J, JHI, JLO, K, L, M, NH, NZ
00137 *     ..
00138 *     .. Local Arrays ..
00139       COMPLEX            V( 2 )
00140 *     ..
00141 *     .. External Functions ..
00142       COMPLEX            CLADIV
00143       REAL               SLAMCH
00144       EXTERNAL           CLADIV, SLAMCH
00145 *     ..
00146 *     .. External Subroutines ..
00147       EXTERNAL           CCOPY, CLARFG, CSCAL, SLABAD
00148 *     ..
00149 *     .. Statement Functions ..
00150       REAL               CABS1
00151 *     ..
00152 *     .. Intrinsic Functions ..
00153       INTRINSIC          ABS, AIMAG, CONJG, MAX, MIN, REAL, SQRT
00154 *     ..
00155 *     .. Statement Function definitions ..
00156       CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
00157 *     ..
00158 *     .. Executable Statements ..
00159 *
00160       INFO = 0
00161 *
00162 *     Quick return if possible
00163 *
00164       IF( N.EQ.0 )
00165      $   RETURN
00166       IF( ILO.EQ.IHI ) THEN
00167          W( ILO ) = H( ILO, ILO )
00168          RETURN
00169       END IF
00170 *
00171 *     ==== clear out the trash ====
00172       DO 10 J = ILO, IHI - 3
00173          H( J+2, J ) = ZERO
00174          H( J+3, J ) = ZERO
00175    10 CONTINUE
00176       IF( ILO.LE.IHI-2 )
00177      $   H( IHI, IHI-2 ) = ZERO
00178 *     ==== ensure that subdiagonal entries are real ====
00179       IF( WANTT ) THEN
00180          JLO = 1
00181          JHI = N
00182       ELSE
00183          JLO = ILO
00184          JHI = IHI
00185       END IF
00186       DO 20 I = ILO + 1, IHI
00187          IF( AIMAG( H( I, I-1 ) ).NE.RZERO ) THEN
00188 *           ==== The following redundant normalization
00189 *           .    avoids problems with both gradual and
00190 *           .    sudden underflow in ABS(H(I,I-1)) ====
00191             SC = H( I, I-1 ) / CABS1( H( I, I-1 ) )
00192             SC = CONJG( SC ) / ABS( SC )
00193             H( I, I-1 ) = ABS( H( I, I-1 ) )
00194             CALL CSCAL( JHI-I+1, SC, H( I, I ), LDH )
00195             CALL CSCAL( MIN( JHI, I+1 )-JLO+1, CONJG( SC ), H( JLO, I ),
00196      $                  1 )
00197             IF( WANTZ )
00198      $         CALL CSCAL( IHIZ-ILOZ+1, CONJG( SC ), Z( ILOZ, I ), 1 )
00199          END IF
00200    20 CONTINUE
00201 *
00202       NH = IHI - ILO + 1
00203       NZ = IHIZ - ILOZ + 1
00204 *
00205 *     Set machine-dependent constants for the stopping criterion.
00206 *
00207       SAFMIN = SLAMCH( 'SAFE MINIMUM' )
00208       SAFMAX = RONE / SAFMIN
00209       CALL SLABAD( SAFMIN, SAFMAX )
00210       ULP = SLAMCH( 'PRECISION' )
00211       SMLNUM = SAFMIN*( REAL( NH ) / ULP )
00212 *
00213 *     I1 and I2 are the indices of the first row and last column of H
00214 *     to which transformations must be applied. If eigenvalues only are
00215 *     being computed, I1 and I2 are set inside the main loop.
00216 *
00217       IF( WANTT ) THEN
00218          I1 = 1
00219          I2 = N
00220       END IF
00221 *
00222 *     The main loop begins here. I is the loop index and decreases from
00223 *     IHI to ILO in steps of 1. Each iteration of the loop works
00224 *     with the active submatrix in rows and columns L to I.
00225 *     Eigenvalues I+1 to IHI have already converged. Either L = ILO, or
00226 *     H(L,L-1) is negligible so that the matrix splits.
00227 *
00228       I = IHI
00229    30 CONTINUE
00230       IF( I.LT.ILO )
00231      $   GO TO 150
00232 *
00233 *     Perform QR iterations on rows and columns ILO to I until a
00234 *     submatrix of order 1 splits off at the bottom because a
00235 *     subdiagonal element has become negligible.
00236 *
00237       L = ILO
00238       DO 130 ITS = 0, ITMAX
00239 *
00240 *        Look for a single small subdiagonal element.
00241 *
00242          DO 40 K = I, L + 1, -1
00243             IF( CABS1( H( K, K-1 ) ).LE.SMLNUM )
00244      $         GO TO 50
00245             TST = CABS1( H( K-1, K-1 ) ) + CABS1( H( K, K ) )
00246             IF( TST.EQ.ZERO ) THEN
00247                IF( K-2.GE.ILO )
00248      $            TST = TST + ABS( REAL( H( K-1, K-2 ) ) )
00249                IF( K+1.LE.IHI )
00250      $            TST = TST + ABS( REAL( H( K+1, K ) ) )
00251             END IF
00252 *           ==== The following is a conservative small subdiagonal
00253 *           .    deflation criterion due to Ahues & Tisseur (LAWN 122,
00254 *           .    1997). It has better mathematical foundation and
00255 *           .    improves accuracy in some examples.  ====
00256             IF( ABS( REAL( H( K, K-1 ) ) ).LE.ULP*TST ) THEN
00257                AB = MAX( CABS1( H( K, K-1 ) ), CABS1( H( K-1, K ) ) )
00258                BA = MIN( CABS1( H( K, K-1 ) ), CABS1( H( K-1, K ) ) )
00259                AA = MAX( CABS1( H( K, K ) ),
00260      $              CABS1( H( K-1, K-1 )-H( K, K ) ) )
00261                BB = MIN( CABS1( H( K, K ) ),
00262      $              CABS1( H( K-1, K-1 )-H( K, K ) ) )
00263                S = AA + AB
00264                IF( BA*( AB / S ).LE.MAX( SMLNUM,
00265      $             ULP*( BB*( AA / S ) ) ) )GO TO 50
00266             END IF
00267    40    CONTINUE
00268    50    CONTINUE
00269          L = K
00270          IF( L.GT.ILO ) THEN
00271 *
00272 *           H(L,L-1) is negligible
00273 *
00274             H( L, L-1 ) = ZERO
00275          END IF
00276 *
00277 *        Exit from loop if a submatrix of order 1 has split off.
00278 *
00279          IF( L.GE.I )
00280      $      GO TO 140
00281 *
00282 *        Now the active submatrix is in rows and columns L to I. If
00283 *        eigenvalues only are being computed, only the active submatrix
00284 *        need be transformed.
00285 *
00286          IF( .NOT.WANTT ) THEN
00287             I1 = L
00288             I2 = I
00289          END IF
00290 *
00291          IF( ITS.EQ.10 ) THEN
00292 *
00293 *           Exceptional shift.
00294 *
00295             S = DAT1*ABS( REAL( H( L+1, L ) ) )
00296             T = S + H( L, L )
00297          ELSE IF( ITS.EQ.20 ) THEN
00298 *
00299 *           Exceptional shift.
00300 *
00301             S = DAT1*ABS( REAL( H( I, I-1 ) ) )
00302             T = S + H( I, I )
00303          ELSE
00304 *
00305 *           Wilkinson's shift.
00306 *
00307             T = H( I, I )
00308             U = SQRT( H( I-1, I ) )*SQRT( H( I, I-1 ) )
00309             S = CABS1( U )
00310             IF( S.NE.RZERO ) THEN
00311                X = HALF*( H( I-1, I-1 )-T )
00312                SX = CABS1( X )
00313                S = MAX( S, CABS1( X ) )
00314                Y = S*SQRT( ( X / S )**2+( U / S )**2 )
00315                IF( SX.GT.RZERO ) THEN
00316                   IF( REAL( X / SX )*REAL( Y )+AIMAG( X / SX )*
00317      $                AIMAG( Y ).LT.RZERO )Y = -Y
00318                END IF
00319                T = T - U*CLADIV( U, ( X+Y ) )
00320             END IF
00321          END IF
00322 *
00323 *        Look for two consecutive small subdiagonal elements.
00324 *
00325          DO 60 M = I - 1, L + 1, -1
00326 *
00327 *           Determine the effect of starting the single-shift QR
00328 *           iteration at row M, and see if this would make H(M,M-1)
00329 *           negligible.
00330 *
00331             H11 = H( M, M )
00332             H22 = H( M+1, M+1 )
00333             H11S = H11 - T
00334             H21 = REAL( H( M+1, M ) )
00335             S = CABS1( H11S ) + ABS( H21 )
00336             H11S = H11S / S
00337             H21 = H21 / S
00338             V( 1 ) = H11S
00339             V( 2 ) = H21
00340             H10 = REAL( H( M, M-1 ) )
00341             IF( ABS( H10 )*ABS( H21 ).LE.ULP*
00342      $          ( CABS1( H11S )*( CABS1( H11 )+CABS1( H22 ) ) ) )
00343      $          GO TO 70
00344    60    CONTINUE
00345          H11 = H( L, L )
00346          H22 = H( L+1, L+1 )
00347          H11S = H11 - T
00348          H21 = REAL( H( L+1, L ) )
00349          S = CABS1( H11S ) + ABS( H21 )
00350          H11S = H11S / S
00351          H21 = H21 / S
00352          V( 1 ) = H11S
00353          V( 2 ) = H21
00354    70    CONTINUE
00355 *
00356 *        Single-shift QR step
00357 *
00358          DO 120 K = M, I - 1
00359 *
00360 *           The first iteration of this loop determines a reflection G
00361 *           from the vector V and applies it from left and right to H,
00362 *           thus creating a nonzero bulge below the subdiagonal.
00363 *
00364 *           Each subsequent iteration determines a reflection G to
00365 *           restore the Hessenberg form in the (K-1)th column, and thus
00366 *           chases the bulge one step toward the bottom of the active
00367 *           submatrix.
00368 *
00369 *           V(2) is always real before the call to CLARFG, and hence
00370 *           after the call T2 ( = T1*V(2) ) is also real.
00371 *
00372             IF( K.GT.M )
00373      $         CALL CCOPY( 2, H( K, K-1 ), 1, V, 1 )
00374             CALL CLARFG( 2, V( 1 ), V( 2 ), 1, T1 )
00375             IF( K.GT.M ) THEN
00376                H( K, K-1 ) = V( 1 )
00377                H( K+1, K-1 ) = ZERO
00378             END IF
00379             V2 = V( 2 )
00380             T2 = REAL( T1*V2 )
00381 *
00382 *           Apply G from the left to transform the rows of the matrix
00383 *           in columns K to I2.
00384 *
00385             DO 80 J = K, I2
00386                SUM = CONJG( T1 )*H( K, J ) + T2*H( K+1, J )
00387                H( K, J ) = H( K, J ) - SUM
00388                H( K+1, J ) = H( K+1, J ) - SUM*V2
00389    80       CONTINUE
00390 *
00391 *           Apply G from the right to transform the columns of the
00392 *           matrix in rows I1 to min(K+2,I).
00393 *
00394             DO 90 J = I1, MIN( K+2, I )
00395                SUM = T1*H( J, K ) + T2*H( J, K+1 )
00396                H( J, K ) = H( J, K ) - SUM
00397                H( J, K+1 ) = H( J, K+1 ) - SUM*CONJG( V2 )
00398    90       CONTINUE
00399 *
00400             IF( WANTZ ) THEN
00401 *
00402 *              Accumulate transformations in the matrix Z
00403 *
00404                DO 100 J = ILOZ, IHIZ
00405                   SUM = T1*Z( J, K ) + T2*Z( J, K+1 )
00406                   Z( J, K ) = Z( J, K ) - SUM
00407                   Z( J, K+1 ) = Z( J, K+1 ) - SUM*CONJG( V2 )
00408   100          CONTINUE
00409             END IF
00410 *
00411             IF( K.EQ.M .AND. M.GT.L ) THEN
00412 *
00413 *              If the QR step was started at row M > L because two
00414 *              consecutive small subdiagonals were found, then extra
00415 *              scaling must be performed to ensure that H(M,M-1) remains
00416 *              real.
00417 *
00418                TEMP = ONE - T1
00419                TEMP = TEMP / ABS( TEMP )
00420                H( M+1, M ) = H( M+1, M )*CONJG( TEMP )
00421                IF( M+2.LE.I )
00422      $            H( M+2, M+1 ) = H( M+2, M+1 )*TEMP
00423                DO 110 J = M, I
00424                   IF( J.NE.M+1 ) THEN
00425                      IF( I2.GT.J )
00426      $                  CALL CSCAL( I2-J, TEMP, H( J, J+1 ), LDH )
00427                      CALL CSCAL( J-I1, CONJG( TEMP ), H( I1, J ), 1 )
00428                      IF( WANTZ ) THEN
00429                         CALL CSCAL( NZ, CONJG( TEMP ), Z( ILOZ, J ), 1 )
00430                      END IF
00431                   END IF
00432   110          CONTINUE
00433             END IF
00434   120    CONTINUE
00435 *
00436 *        Ensure that H(I,I-1) is real.
00437 *
00438          TEMP = H( I, I-1 )
00439          IF( AIMAG( TEMP ).NE.RZERO ) THEN
00440             RTEMP = ABS( TEMP )
00441             H( I, I-1 ) = RTEMP
00442             TEMP = TEMP / RTEMP
00443             IF( I2.GT.I )
00444      $         CALL CSCAL( I2-I, CONJG( TEMP ), H( I, I+1 ), LDH )
00445             CALL CSCAL( I-I1, TEMP, H( I1, I ), 1 )
00446             IF( WANTZ ) THEN
00447                CALL CSCAL( NZ, TEMP, Z( ILOZ, I ), 1 )
00448             END IF
00449          END IF
00450 *
00451   130 CONTINUE
00452 *
00453 *     Failure to converge in remaining number of iterations
00454 *
00455       INFO = I
00456       RETURN
00457 *
00458   140 CONTINUE
00459 *
00460 *     H(I,I-1) is negligible: one eigenvalue has converged.
00461 *
00462       W( I ) = H( I, I )
00463 *
00464 *     return to start of the main loop with new value of I.
00465 *
00466       I = L - 1
00467       GO TO 30
00468 *
00469   150 CONTINUE
00470       RETURN
00471 *
00472 *     End of CLAHQR
00473 *
00474       END