LAPACK: SRC/dlahqr.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
00002      $                   ILOZ, 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       DOUBLE PRECISION   H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )
00014 *     ..
00015 *
00016 *     Purpose
00017 *     =======
00018 *
00019 *     DLAHQR is an auxiliary routine called by DHSEQR to update the
00020 *     eigenvalues and Schur decomposition already computed by DHSEQR, 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 quasi-triangular in
00041 *          rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
00042 *          ILO = 1). DLAHQR works primarily with the Hessenberg
00043 *          submatrix in rows and columns ILO to IHI, but applies
00044 *          transformations to all of H if WANTT is .TRUE..
00045 *          1 <= ILO <= max(1,IHI); IHI <= N.
00046 *
00047 *     H       (input/output) DOUBLE PRECISION array, dimension (LDH,N)
00048 *          On entry, the upper Hessenberg matrix H.
00049 *          On exit, if INFO is zero and if WANTT is .TRUE., H is upper
00050 *          quasi-triangular in rows and columns ILO:IHI, with any
00051 *          2-by-2 diagonal blocks in standard form. If INFO is zero
00052 *          and WANTT is .FALSE., the contents of H are unspecified on
00053 *          exit.  The output state of H if INFO is nonzero is given
00054 *          below under the description of INFO.
00055 *
00056 *     LDH     (input) INTEGER
00057 *          The leading dimension of the array H. LDH >= max(1,N).
00058 *
00059 *     WR      (output) DOUBLE PRECISION array, dimension (N)
00060 *     WI      (output) DOUBLE PRECISION array, dimension (N)
00061 *          The real and imaginary parts, respectively, of the computed
00062 *          eigenvalues ILO to IHI are stored in the corresponding
00063 *          elements of WR and WI. If two eigenvalues are computed as a
00064 *          complex conjugate pair, they are stored in consecutive
00065 *          elements of WR and WI, say the i-th and (i+1)th, with
00066 *          WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the
00067 *          eigenvalues are stored in the same order as on the diagonal
00068 *          of the Schur form returned in H, with WR(i) = H(i,i), and, if
00069 *          H(i:i+1,i:i+1) is a 2-by-2 diagonal block,
00070 *          WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).
00071 *
00072 *     ILOZ    (input) INTEGER
00073 *     IHIZ    (input) INTEGER
00074 *          Specify the rows of Z to which transformations must be
00075 *          applied if WANTZ is .TRUE..
00076 *          1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
00077 *
00078 *     Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
00079 *          If WANTZ is .TRUE., on entry Z must contain the current
00080 *          matrix Z of transformations accumulated by DHSEQR, and on
00081 *          exit Z has been updated; transformations are applied only to
00082 *          the submatrix Z(ILOZ:IHIZ,ILO:IHI).
00083 *          If WANTZ is .FALSE., Z is not referenced.
00084 *
00085 *     LDZ     (input) INTEGER
00086 *          The leading dimension of the array Z. LDZ >= max(1,N).
00087 *
00088 *     INFO    (output) INTEGER
00089 *           =   0: successful exit
00090 *          .GT. 0: If INFO = i, DLAHQR failed to compute all the
00091 *                  eigenvalues ILO to IHI in a total of 30 iterations
00092 *                  per eigenvalue; elements i+1:ihi of WR and WI
00093 *                  contain those eigenvalues which have been
00094 *                  successfully computed.
00095 *
00096 *                  If INFO .GT. 0 and WANTT is .FALSE., then on exit,
00097 *                  the remaining unconverged eigenvalues are the
00098 *                  eigenvalues of the upper Hessenberg matrix rows
00099 *                  and columns ILO thorugh INFO of the final, output
00100 *                  value of H.
00101 *
00102 *                  If INFO .GT. 0 and WANTT is .TRUE., then on exit
00103 *          (*)       (initial value of H)*U  = U*(final value of H)
00104 *                  where U is an orthognal matrix.    The final
00105 *                  value of H is upper Hessenberg and triangular in
00106 *                  rows and columns INFO+1 through IHI.
00107 *
00108 *                  If INFO .GT. 0 and WANTZ is .TRUE., then on exit
00109 *                      (final value of Z)  = (initial value of Z)*U
00110 *                  where U is the orthogonal matrix in (*)
00111 *                  (regardless of the value of WANTT.)
00112 *
00113 *     Further Details
00114 *     ===============
00115 *
00116 *     02-96 Based on modifications by
00117 *     David Day, Sandia National Laboratory, USA
00118 *
00119 *     12-04 Further modifications by
00120 *     Ralph Byers, University of Kansas, USA
00121 *     This is a modified version of DLAHQR from LAPACK version 3.0.
00122 *     It is (1) more robust against overflow and underflow and
00123 *     (2) adopts the more conservative Ahues & Tisseur stopping
00124 *     criterion (LAWN 122, 1997).
00125 *
00126 *     =========================================================
00127 *
00128 *     .. Parameters ..
00129       INTEGER            ITMAX
00130       PARAMETER          ( ITMAX = 30 )
00131       DOUBLE PRECISION   ZERO, ONE, TWO
00132       PARAMETER          ( ZERO = 0.0d0, ONE = 1.0d0, TWO = 2.0d0 )
00133       DOUBLE PRECISION   DAT1, DAT2
00134       PARAMETER          ( DAT1 = 3.0d0 / 4.0d0, DAT2 = -0.4375d0 )
00135 *     ..
00136 *     .. Local Scalars ..
00137       DOUBLE PRECISION   AA, AB, BA, BB, CS, DET, H11, H12, H21, H21S,
00138      $                   H22, RT1I, RT1R, RT2I, RT2R, RTDISC, S, SAFMAX,
00139      $                   SAFMIN, SMLNUM, SN, SUM, T1, T2, T3, TR, TST,
00140      $                   ULP, V2, V3
00141       INTEGER            I, I1, I2, ITS, J, K, L, M, NH, NR, NZ
00142 *     ..
00143 *     .. Local Arrays ..
00144       DOUBLE PRECISION   V( 3 )
00145 *     ..
00146 *     .. External Functions ..
00147       DOUBLE PRECISION   DLAMCH
00148       EXTERNAL           DLAMCH
00149 *     ..
00150 *     .. External Subroutines ..
00151       EXTERNAL           DCOPY, DLABAD, DLANV2, DLARFG, DROT
00152 *     ..
00153 *     .. Intrinsic Functions ..
00154       INTRINSIC          ABS, DBLE, MAX, MIN, SQRT
00155 *     ..
00156 *     .. Executable Statements ..
00157 *
00158       INFO = 0
00159 *
00160 *     Quick return if possible
00161 *
00162       IF( N.EQ.0 )
00163      $   RETURN
00164       IF( ILO.EQ.IHI ) THEN
00165          WR( ILO ) = H( ILO, ILO )
00166          WI( ILO ) = ZERO
00167          RETURN
00168       END IF
00169 *
00170 *     ==== clear out the trash ====
00171       DO 10 J = ILO, IHI - 3
00172          H( J+2, J ) = ZERO
00173          H( J+3, J ) = ZERO
00174    10 CONTINUE
00175       IF( ILO.LE.IHI-2 )
00176      $   H( IHI, IHI-2 ) = ZERO
00177 *
00178       NH = IHI - ILO + 1
00179       NZ = IHIZ - ILOZ + 1
00180 *
00181 *     Set machine-dependent constants for the stopping criterion.
00182 *
00183       SAFMIN = DLAMCH( 'SAFE MINIMUM' )
00184       SAFMAX = ONE / SAFMIN
00185       CALL DLABAD( SAFMIN, SAFMAX )
00186       ULP = DLAMCH( 'PRECISION' )
00187       SMLNUM = SAFMIN*( DBLE( NH ) / ULP )
00188 *
00189 *     I1 and I2 are the indices of the first row and last column of H
00190 *     to which transformations must be applied. If eigenvalues only are
00191 *     being computed, I1 and I2 are set inside the main loop.
00192 *
00193       IF( WANTT ) THEN
00194          I1 = 1
00195          I2 = N
00196       END IF
00197 *
00198 *     The main loop begins here. I is the loop index and decreases from
00199 *     IHI to ILO in steps of 1 or 2. Each iteration of the loop works
00200 *     with the active submatrix in rows and columns L to I.
00201 *     Eigenvalues I+1 to IHI have already converged. Either L = ILO or
00202 *     H(L,L-1) is negligible so that the matrix splits.
00203 *
00204       I = IHI
00205    20 CONTINUE
00206       L = ILO
00207       IF( I.LT.ILO )
00208      $   GO TO 160
00209 *
00210 *     Perform QR iterations on rows and columns ILO to I until a
00211 *     submatrix of order 1 or 2 splits off at the bottom because a
00212 *     subdiagonal element has become negligible.
00213 *
00214       DO 140 ITS = 0, ITMAX
00215 *
00216 *        Look for a single small subdiagonal element.
00217 *
00218          DO 30 K = I, L + 1, -1
00219             IF( ABS( H( K, K-1 ) ).LE.SMLNUM )
00220      $         GO TO 40
00221             TST = ABS( H( K-1, K-1 ) ) + ABS( H( K, K ) )
00222             IF( TST.EQ.ZERO ) THEN
00223                IF( K-2.GE.ILO )
00224      $            TST = TST + ABS( H( K-1, K-2 ) )
00225                IF( K+1.LE.IHI )
00226      $            TST = TST + ABS( H( K+1, K ) )
00227             END IF
00228 *           ==== The following is a conservative small subdiagonal
00229 *           .    deflation  criterion due to Ahues & Tisseur (LAWN 122,
00230 *           .    1997). It has better mathematical foundation and
00231 *           .    improves accuracy in some cases.  ====
00232             IF( ABS( H( K, K-1 ) ).LE.ULP*TST ) THEN
00233                AB = MAX( ABS( H( K, K-1 ) ), ABS( H( K-1, K ) ) )
00234                BA = MIN( ABS( H( K, K-1 ) ), ABS( H( K-1, K ) ) )
00235                AA = MAX( ABS( H( K, K ) ),
00236      $              ABS( H( K-1, K-1 )-H( K, K ) ) )
00237                BB = MIN( ABS( H( K, K ) ),
00238      $              ABS( H( K-1, K-1 )-H( K, K ) ) )
00239                S = AA + AB
00240                IF( BA*( AB / S ).LE.MAX( SMLNUM,
00241      $             ULP*( BB*( AA / S ) ) ) )GO TO 40
00242             END IF
00243    30    CONTINUE
00244    40    CONTINUE
00245          L = K
00246          IF( L.GT.ILO ) THEN
00247 *
00248 *           H(L,L-1) is negligible
00249 *
00250             H( L, L-1 ) = ZERO
00251          END IF
00252 *
00253 *        Exit from loop if a submatrix of order 1 or 2 has split off.
00254 *
00255          IF( L.GE.I-1 )
00256      $      GO TO 150
00257 *
00258 *        Now the active submatrix is in rows and columns L to I. If
00259 *        eigenvalues only are being computed, only the active submatrix
00260 *        need be transformed.
00261 *
00262          IF( .NOT.WANTT ) THEN
00263             I1 = L
00264             I2 = I
00265          END IF
00266 *
00267          IF( ITS.EQ.10 ) THEN
00268 *
00269 *           Exceptional shift.
00270 *
00271             S = ABS( H( L+1, L ) ) + ABS( H( L+2, L+1 ) )
00272             H11 = DAT1*S + H( L, L )
00273             H12 = DAT2*S
00274             H21 = S
00275             H22 = H11
00276          ELSE IF( ITS.EQ.20 ) THEN
00277 *
00278 *           Exceptional shift.
00279 *
00280             S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) )
00281             H11 = DAT1*S + H( I, I )
00282             H12 = DAT2*S
00283             H21 = S
00284             H22 = H11
00285          ELSE
00286 *
00287 *           Prepare to use Francis' double shift
00288 *           (i.e. 2nd degree generalized Rayleigh quotient)
00289 *
00290             H11 = H( I-1, I-1 )
00291             H21 = H( I, I-1 )
00292             H12 = H( I-1, I )
00293             H22 = H( I, I )
00294          END IF
00295          S = ABS( H11 ) + ABS( H12 ) + ABS( H21 ) + ABS( H22 )
00296          IF( S.EQ.ZERO ) THEN
00297             RT1R = ZERO
00298             RT1I = ZERO
00299             RT2R = ZERO
00300             RT2I = ZERO
00301          ELSE
00302             H11 = H11 / S
00303             H21 = H21 / S
00304             H12 = H12 / S
00305             H22 = H22 / S
00306             TR = ( H11+H22 ) / TWO
00307             DET = ( H11-TR )*( H22-TR ) - H12*H21
00308             RTDISC = SQRT( ABS( DET ) )
00309             IF( DET.GE.ZERO ) THEN
00310 *
00311 *              ==== complex conjugate shifts ====
00312 *
00313                RT1R = TR*S
00314                RT2R = RT1R
00315                RT1I = RTDISC*S
00316                RT2I = -RT1I
00317             ELSE
00318 *
00319 *              ==== real shifts (use only one of them)  ====
00320 *
00321                RT1R = TR + RTDISC
00322                RT2R = TR - RTDISC
00323                IF( ABS( RT1R-H22 ).LE.ABS( RT2R-H22 ) ) THEN
00324                   RT1R = RT1R*S
00325                   RT2R = RT1R
00326                ELSE
00327                   RT2R = RT2R*S
00328                   RT1R = RT2R
00329                END IF
00330                RT1I = ZERO
00331                RT2I = ZERO
00332             END IF
00333          END IF
00334 *
00335 *        Look for two consecutive small subdiagonal elements.
00336 *
00337          DO 50 M = I - 2, L, -1
00338 *           Determine the effect of starting the double-shift QR
00339 *           iteration at row M, and see if this would make H(M,M-1)
00340 *           negligible.  (The following uses scaling to avoid
00341 *           overflows and most underflows.)
00342 *
00343             H21S = H( M+1, M )
00344             S = ABS( H( M, M )-RT2R ) + ABS( RT2I ) + ABS( H21S )
00345             H21S = H( M+1, M ) / S
00346             V( 1 ) = H21S*H( M, M+1 ) + ( H( M, M )-RT1R )*
00347      $               ( ( H( M, M )-RT2R ) / S ) - RT1I*( RT2I / S )
00348             V( 2 ) = H21S*( H( M, M )+H( M+1, M+1 )-RT1R-RT2R )
00349             V( 3 ) = H21S*H( M+2, M+1 )
00350             S = ABS( V( 1 ) ) + ABS( V( 2 ) ) + ABS( V( 3 ) )
00351             V( 1 ) = V( 1 ) / S
00352             V( 2 ) = V( 2 ) / S
00353             V( 3 ) = V( 3 ) / S
00354             IF( M.EQ.L )
00355      $         GO TO 60
00356             IF( ABS( H( M, M-1 ) )*( ABS( V( 2 ) )+ABS( V( 3 ) ) ).LE.
00357      $          ULP*ABS( V( 1 ) )*( ABS( H( M-1, M-1 ) )+ABS( H( M,
00358      $          M ) )+ABS( H( M+1, M+1 ) ) ) )GO TO 60
00359    50    CONTINUE
00360    60    CONTINUE
00361 *
00362 *        Double-shift QR step
00363 *
00364          DO 130 K = M, I - 1
00365 *
00366 *           The first iteration of this loop determines a reflection G
00367 *           from the vector V and applies it from left and right to H,
00368 *           thus creating a nonzero bulge below the subdiagonal.
00369 *
00370 *           Each subsequent iteration determines a reflection G to
00371 *           restore the Hessenberg form in the (K-1)th column, and thus
00372 *           chases the bulge one step toward the bottom of the active
00373 *           submatrix. NR is the order of G.
00374 *
00375             NR = MIN( 3, I-K+1 )
00376             IF( K.GT.M )
00377      $         CALL DCOPY( NR, H( K, K-1 ), 1, V, 1 )
00378             CALL DLARFG( NR, V( 1 ), V( 2 ), 1, T1 )
00379             IF( K.GT.M ) THEN
00380                H( K, K-1 ) = V( 1 )
00381                H( K+1, K-1 ) = ZERO
00382                IF( K.LT.I-1 )
00383      $            H( K+2, K-1 ) = ZERO
00384             ELSE IF( M.GT.L ) THEN
00385 *               ==== Use the following instead of
00386 *               .    H( K, K-1 ) = -H( K, K-1 ) to
00387 *               .    avoid a bug when v(2) and v(3)
00388 *               .    underflow. ====
00389                H( K, K-1 ) = H( K, K-1 )*( ONE-T1 )
00390             END IF
00391             V2 = V( 2 )
00392             T2 = T1*V2
00393             IF( NR.EQ.3 ) THEN
00394                V3 = V( 3 )
00395                T3 = T1*V3
00396 *
00397 *              Apply G from the left to transform the rows of the matrix
00398 *              in columns K to I2.
00399 *
00400                DO 70 J = K, I2
00401                   SUM = H( K, J ) + V2*H( K+1, J ) + V3*H( K+2, J )
00402                   H( K, J ) = H( K, J ) - SUM*T1
00403                   H( K+1, J ) = H( K+1, J ) - SUM*T2
00404                   H( K+2, J ) = H( K+2, J ) - SUM*T3
00405    70          CONTINUE
00406 *
00407 *              Apply G from the right to transform the columns of the
00408 *              matrix in rows I1 to min(K+3,I).
00409 *
00410                DO 80 J = I1, MIN( K+3, I )
00411                   SUM = H( J, K ) + V2*H( J, K+1 ) + V3*H( J, K+2 )
00412                   H( J, K ) = H( J, K ) - SUM*T1
00413                   H( J, K+1 ) = H( J, K+1 ) - SUM*T2
00414                   H( J, K+2 ) = H( J, K+2 ) - SUM*T3
00415    80          CONTINUE
00416 *
00417                IF( WANTZ ) THEN
00418 *
00419 *                 Accumulate transformations in the matrix Z
00420 *
00421                   DO 90 J = ILOZ, IHIZ
00422                      SUM = Z( J, K ) + V2*Z( J, K+1 ) + V3*Z( J, K+2 )
00423                      Z( J, K ) = Z( J, K ) - SUM*T1
00424                      Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2
00425                      Z( J, K+2 ) = Z( J, K+2 ) - SUM*T3
00426    90             CONTINUE
00427                END IF
00428             ELSE IF( NR.EQ.2 ) THEN
00429 *
00430 *              Apply G from the left to transform the rows of the matrix
00431 *              in columns K to I2.
00432 *
00433                DO 100 J = K, I2
00434                   SUM = H( K, J ) + V2*H( K+1, J )
00435                   H( K, J ) = H( K, J ) - SUM*T1
00436                   H( K+1, J ) = H( K+1, J ) - SUM*T2
00437   100          CONTINUE
00438 *
00439 *              Apply G from the right to transform the columns of the
00440 *              matrix in rows I1 to min(K+3,I).
00441 *
00442                DO 110 J = I1, I
00443                   SUM = H( J, K ) + V2*H( J, K+1 )
00444                   H( J, K ) = H( J, K ) - SUM*T1
00445                   H( J, K+1 ) = H( J, K+1 ) - SUM*T2
00446   110          CONTINUE
00447 *
00448                IF( WANTZ ) THEN
00449 *
00450 *                 Accumulate transformations in the matrix Z
00451 *
00452                   DO 120 J = ILOZ, IHIZ
00453                      SUM = Z( J, K ) + V2*Z( J, K+1 )
00454                      Z( J, K ) = Z( J, K ) - SUM*T1
00455                      Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2
00456   120             CONTINUE
00457                END IF
00458             END IF
00459   130    CONTINUE
00460 *
00461   140 CONTINUE
00462 *
00463 *     Failure to converge in remaining number of iterations
00464 *
00465       INFO = I
00466       RETURN
00467 *
00468   150 CONTINUE
00469 *
00470       IF( L.EQ.I ) THEN
00471 *
00472 *        H(I,I-1) is negligible: one eigenvalue has converged.
00473 *
00474          WR( I ) = H( I, I )
00475          WI( I ) = ZERO
00476       ELSE IF( L.EQ.I-1 ) THEN
00477 *
00478 *        H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
00479 *
00480 *        Transform the 2-by-2 submatrix to standard Schur form,
00481 *        and compute and store the eigenvalues.
00482 *
00483          CALL DLANV2( H( I-1, I-1 ), H( I-1, I ), H( I, I-1 ),
00484      $                H( I, I ), WR( I-1 ), WI( I-1 ), WR( I ), WI( I ),
00485      $                CS, SN )
00486 *
00487          IF( WANTT ) THEN
00488 *
00489 *           Apply the transformation to the rest of H.
00490 *
00491             IF( I2.GT.I )
00492      $         CALL DROT( I2-I, H( I-1, I+1 ), LDH, H( I, I+1 ), LDH,
00493      $                    CS, SN )
00494             CALL DROT( I-I1-1, H( I1, I-1 ), 1, H( I1, I ), 1, CS, SN )
00495          END IF
00496          IF( WANTZ ) THEN
00497 *
00498 *           Apply the transformation to Z.
00499 *
00500             CALL DROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, SN )
00501          END IF
00502       END IF
00503 *
00504 *     return to start of the main loop with new value of I.
00505 *
00506       I = L - 1
00507       GO TO 20
00508 *
00509   160 CONTINUE
00510       RETURN
00511 *
00512 *     End of DLAHQR
00513 *
00514       END