LAPACK 3.3.0

dhseqr.f

Go to the documentation of this file.
00001       SUBROUTINE DHSEQR( JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z,
00002      $                   LDZ, WORK, LWORK, INFO )
00003 *
00004 *  -- LAPACK computational routine (version 3.2.2) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
00006 *     June 2010
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            IHI, ILO, INFO, LDH, LDZ, LWORK, N
00010       CHARACTER          COMPZ, JOB
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   H( LDH, * ), WI( * ), WORK( * ), WR( * ),
00014      $                   Z( LDZ, * )
00015 *     ..
00016 *     Purpose
00017 *     =======
00018 *
00019 *     DHSEQR computes the eigenvalues of a Hessenberg matrix H
00020 *     and, optionally, the matrices T and Z from the Schur decomposition
00021 *     H = Z T Z**T, where T is an upper quasi-triangular matrix (the
00022 *     Schur form), and Z is the orthogonal matrix of Schur vectors.
00023 *
00024 *     Optionally Z may be postmultiplied into an input orthogonal
00025 *     matrix Q so that this routine can give the Schur factorization
00026 *     of a matrix A which has been reduced to the Hessenberg form H
00027 *     by the orthogonal matrix Q:  A = Q*H*Q**T = (QZ)*T*(QZ)**T.
00028 *
00029 *     Arguments
00030 *     =========
00031 *
00032 *     JOB   (input) CHARACTER*1
00033 *           = 'E':  compute eigenvalues only;
00034 *           = 'S':  compute eigenvalues and the Schur form T.
00035 *
00036 *     COMPZ (input) CHARACTER*1
00037 *           = 'N':  no Schur vectors are computed;
00038 *           = 'I':  Z is initialized to the unit matrix and the matrix Z
00039 *                   of Schur vectors of H is returned;
00040 *           = 'V':  Z must contain an orthogonal matrix Q on entry, and
00041 *                   the product Q*Z is returned.
00042 *
00043 *     N     (input) INTEGER
00044 *           The order of the matrix H.  N .GE. 0.
00045 *
00046 *     ILO   (input) INTEGER
00047 *     IHI   (input) INTEGER
00048 *           It is assumed that H is already upper triangular in rows
00049 *           and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
00050 *           set by a previous call to DGEBAL, and then passed to DGEHRD
00051 *           when the matrix output by DGEBAL is reduced to Hessenberg
00052 *           form. Otherwise ILO and IHI should be set to 1 and N
00053 *           respectively.  If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
00054 *           If N = 0, then ILO = 1 and IHI = 0.
00055 *
00056 *     H     (input/output) DOUBLE PRECISION array, dimension (LDH,N)
00057 *           On entry, the upper Hessenberg matrix H.
00058 *           On exit, if INFO = 0 and JOB = 'S', then H contains the
00059 *           upper quasi-triangular matrix T from the Schur decomposition
00060 *           (the Schur form); 2-by-2 diagonal blocks (corresponding to
00061 *           complex conjugate pairs of eigenvalues) are returned in
00062 *           standard form, with H(i,i) = H(i+1,i+1) and
00063 *           H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and JOB = 'E', the
00064 *           contents of H are unspecified on exit.  (The output value of
00065 *           H when INFO.GT.0 is given under the description of INFO
00066 *           below.)
00067 *
00068 *           Unlike earlier versions of DHSEQR, this subroutine may
00069 *           explicitly H(i,j) = 0 for i.GT.j and j = 1, 2, ... ILO-1
00070 *           or j = IHI+1, IHI+2, ... N.
00071 *
00072 *     LDH   (input) INTEGER
00073 *           The leading dimension of the array H. LDH .GE. max(1,N).
00074 *
00075 *     WR    (output) DOUBLE PRECISION array, dimension (N)
00076 *     WI    (output) DOUBLE PRECISION array, dimension (N)
00077 *           The real and imaginary parts, respectively, of the computed
00078 *           eigenvalues. If two eigenvalues are computed as a complex
00079 *           conjugate pair, they are stored in consecutive elements of
00080 *           WR and WI, say the i-th and (i+1)th, with WI(i) .GT. 0 and
00081 *           WI(i+1) .LT. 0. If JOB = 'S', the eigenvalues are stored in
00082 *           the same order as on the diagonal of the Schur form returned
00083 *           in H, with WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2
00084 *           diagonal block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and
00085 *           WI(i+1) = -WI(i).
00086 *
00087 *     Z     (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
00088 *           If COMPZ = 'N', Z is not referenced.
00089 *           If COMPZ = 'I', on entry Z need not be set and on exit,
00090 *           if INFO = 0, Z contains the orthogonal matrix Z of the Schur
00091 *           vectors of H.  If COMPZ = 'V', on entry Z must contain an
00092 *           N-by-N matrix Q, which is assumed to be equal to the unit
00093 *           matrix except for the submatrix Z(ILO:IHI,ILO:IHI). On exit,
00094 *           if INFO = 0, Z contains Q*Z.
00095 *           Normally Q is the orthogonal matrix generated by DORGHR
00096 *           after the call to DGEHRD which formed the Hessenberg matrix
00097 *           H. (The output value of Z when INFO.GT.0 is given under
00098 *           the description of INFO below.)
00099 *
00100 *     LDZ   (input) INTEGER
00101 *           The leading dimension of the array Z.  if COMPZ = 'I' or
00102 *           COMPZ = 'V', then LDZ.GE.MAX(1,N).  Otherwize, LDZ.GE.1.
00103 *
00104 *     WORK  (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
00105 *           On exit, if INFO = 0, WORK(1) returns an estimate of
00106 *           the optimal value for LWORK.
00107 *
00108 *     LWORK (input) INTEGER
00109 *           The dimension of the array WORK.  LWORK .GE. max(1,N)
00110 *           is sufficient and delivers very good and sometimes
00111 *           optimal performance.  However, LWORK as large as 11*N
00112 *           may be required for optimal performance.  A workspace
00113 *           query is recommended to determine the optimal workspace
00114 *           size.
00115 *
00116 *           If LWORK = -1, then DHSEQR does a workspace query.
00117 *           In this case, DHSEQR checks the input parameters and
00118 *           estimates the optimal workspace size for the given
00119 *           values of N, ILO and IHI.  The estimate is returned
00120 *           in WORK(1).  No error message related to LWORK is
00121 *           issued by XERBLA.  Neither H nor Z are accessed.
00122 *
00123 *
00124 *     INFO  (output) INTEGER
00125 *             =  0:  successful exit
00126 *           .LT. 0:  if INFO = -i, the i-th argument had an illegal
00127 *                    value
00128 *           .GT. 0:  if INFO = i, DHSEQR failed to compute all of
00129 *                the eigenvalues.  Elements 1:ilo-1 and i+1:n of WR
00130 *                and WI contain those eigenvalues which have been
00131 *                successfully computed.  (Failures are rare.)
00132 *
00133 *                If INFO .GT. 0 and JOB = 'E', then on exit, the
00134 *                remaining unconverged eigenvalues are the eigen-
00135 *                values of the upper Hessenberg matrix rows and
00136 *                columns ILO through INFO of the final, output
00137 *                value of H.
00138 *
00139 *                If INFO .GT. 0 and JOB   = 'S', then on exit
00140 *
00141 *           (*)  (initial value of H)*U  = U*(final value of H)
00142 *
00143 *                where U is an orthogonal matrix.  The final
00144 *                value of H is upper Hessenberg and quasi-triangular
00145 *                in rows and columns INFO+1 through IHI.
00146 *
00147 *                If INFO .GT. 0 and COMPZ = 'V', then on exit
00148 *
00149 *                  (final value of Z)  =  (initial value of Z)*U
00150 *
00151 *                where U is the orthogonal matrix in (*) (regard-
00152 *                less of the value of JOB.)
00153 *
00154 *                If INFO .GT. 0 and COMPZ = 'I', then on exit
00155 *                      (final value of Z)  = U
00156 *                where U is the orthogonal matrix in (*) (regard-
00157 *                less of the value of JOB.)
00158 *
00159 *                If INFO .GT. 0 and COMPZ = 'N', then Z is not
00160 *                accessed.
00161 *
00162 *     ================================================================
00163 *             Default values supplied by
00164 *             ILAENV(ISPEC,'DHSEQR',JOB(:1)//COMPZ(:1),N,ILO,IHI,LWORK).
00165 *             It is suggested that these defaults be adjusted in order
00166 *             to attain best performance in each particular
00167 *             computational environment.
00168 *
00169 *            ISPEC=12: The DLAHQR vs DLAQR0 crossover point.
00170 *                      Default: 75. (Must be at least 11.)
00171 *
00172 *            ISPEC=13: Recommended deflation window size.
00173 *                      This depends on ILO, IHI and NS.  NS is the
00174 *                      number of simultaneous shifts returned
00175 *                      by ILAENV(ISPEC=15).  (See ISPEC=15 below.)
00176 *                      The default for (IHI-ILO+1).LE.500 is NS.
00177 *                      The default for (IHI-ILO+1).GT.500 is 3*NS/2.
00178 *
00179 *            ISPEC=14: Nibble crossover point. (See IPARMQ for
00180 *                      details.)  Default: 14% of deflation window
00181 *                      size.
00182 *
00183 *            ISPEC=15: Number of simultaneous shifts in a multishift
00184 *                      QR iteration.
00185 *
00186 *                      If IHI-ILO+1 is ...
00187 *
00188 *                      greater than      ...but less    ... the
00189 *                      or equal to ...      than        default is
00190 *
00191 *                           1               30          NS =   2(+)
00192 *                          30               60          NS =   4(+)
00193 *                          60              150          NS =  10(+)
00194 *                         150              590          NS =  **
00195 *                         590             3000          NS =  64
00196 *                        3000             6000          NS = 128
00197 *                        6000             infinity      NS = 256
00198 *
00199 *                  (+)  By default some or all matrices of this order
00200 *                       are passed to the implicit double shift routine
00201 *                       DLAHQR and this parameter is ignored.  See
00202 *                       ISPEC=12 above and comments in IPARMQ for
00203 *                       details.
00204 *
00205 *                 (**)  The asterisks (**) indicate an ad-hoc
00206 *                       function of N increasing from 10 to 64.
00207 *
00208 *            ISPEC=16: Select structured matrix multiply.
00209 *                      If the number of simultaneous shifts (specified
00210 *                      by ISPEC=15) is less than 14, then the default
00211 *                      for ISPEC=16 is 0.  Otherwise the default for
00212 *                      ISPEC=16 is 2.
00213 *
00214 *     ================================================================
00215 *     Based on contributions by
00216 *        Karen Braman and Ralph Byers, Department of Mathematics,
00217 *        University of Kansas, USA
00218 *
00219 *     ================================================================
00220 *     References:
00221 *       K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
00222 *       Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
00223 *       Performance, SIAM Journal of Matrix Analysis, volume 23, pages
00224 *       929--947, 2002.
00225 *
00226 *       K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
00227 *       Algorithm Part II: Aggressive Early Deflation, SIAM Journal
00228 *       of Matrix Analysis, volume 23, pages 948--973, 2002.
00229 *
00230 *     ================================================================
00231 *     .. Parameters ..
00232 *
00233 *     ==== Matrices of order NTINY or smaller must be processed by
00234 *     .    DLAHQR because of insufficient subdiagonal scratch space.
00235 *     .    (This is a hard limit.) ====
00236       INTEGER            NTINY
00237       PARAMETER          ( NTINY = 11 )
00238 *
00239 *     ==== NL allocates some local workspace to help small matrices
00240 *     .    through a rare DLAHQR failure.  NL .GT. NTINY = 11 is
00241 *     .    required and NL .LE. NMIN = ILAENV(ISPEC=12,...) is recom-
00242 *     .    mended.  (The default value of NMIN is 75.)  Using NL = 49
00243 *     .    allows up to six simultaneous shifts and a 16-by-16
00244 *     .    deflation window.  ====
00245       INTEGER            NL
00246       PARAMETER          ( NL = 49 )
00247       DOUBLE PRECISION   ZERO, ONE
00248       PARAMETER          ( ZERO = 0.0d0, ONE = 1.0d0 )
00249 *     ..
00250 *     .. Local Arrays ..
00251       DOUBLE PRECISION   HL( NL, NL ), WORKL( NL )
00252 *     ..
00253 *     .. Local Scalars ..
00254       INTEGER            I, KBOT, NMIN
00255       LOGICAL            INITZ, LQUERY, WANTT, WANTZ
00256 *     ..
00257 *     .. External Functions ..
00258       INTEGER            ILAENV
00259       LOGICAL            LSAME
00260       EXTERNAL           ILAENV, LSAME
00261 *     ..
00262 *     .. External Subroutines ..
00263       EXTERNAL           DLACPY, DLAHQR, DLAQR0, DLASET, XERBLA
00264 *     ..
00265 *     .. Intrinsic Functions ..
00266       INTRINSIC          DBLE, MAX, MIN
00267 *     ..
00268 *     .. Executable Statements ..
00269 *
00270 *     ==== Decode and check the input parameters. ====
00271 *
00272       WANTT = LSAME( JOB, 'S' )
00273       INITZ = LSAME( COMPZ, 'I' )
00274       WANTZ = INITZ .OR. LSAME( COMPZ, 'V' )
00275       WORK( 1 ) = DBLE( MAX( 1, N ) )
00276       LQUERY = LWORK.EQ.-1
00277 *
00278       INFO = 0
00279       IF( .NOT.LSAME( JOB, 'E' ) .AND. .NOT.WANTT ) THEN
00280          INFO = -1
00281       ELSE IF( .NOT.LSAME( COMPZ, 'N' ) .AND. .NOT.WANTZ ) THEN
00282          INFO = -2
00283       ELSE IF( N.LT.0 ) THEN
00284          INFO = -3
00285       ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN
00286          INFO = -4
00287       ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN
00288          INFO = -5
00289       ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
00290          INFO = -7
00291       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.MAX( 1, N ) ) ) THEN
00292          INFO = -11
00293       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
00294          INFO = -13
00295       END IF
00296 *
00297       IF( INFO.NE.0 ) THEN
00298 *
00299 *        ==== Quick return in case of invalid argument. ====
00300 *
00301          CALL XERBLA( 'DHSEQR', -INFO )
00302          RETURN
00303 *
00304       ELSE IF( N.EQ.0 ) THEN
00305 *
00306 *        ==== Quick return in case N = 0; nothing to do. ====
00307 *
00308          RETURN
00309 *
00310       ELSE IF( LQUERY ) THEN
00311 *
00312 *        ==== Quick return in case of a workspace query ====
00313 *
00314          CALL DLAQR0( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO,
00315      $                IHI, Z, LDZ, WORK, LWORK, INFO )
00316 *        ==== Ensure reported workspace size is backward-compatible with
00317 *        .    previous LAPACK versions. ====
00318          WORK( 1 ) = MAX( DBLE( MAX( 1, N ) ), WORK( 1 ) )
00319          RETURN
00320 *
00321       ELSE
00322 *
00323 *        ==== copy eigenvalues isolated by DGEBAL ====
00324 *
00325          DO 10 I = 1, ILO - 1
00326             WR( I ) = H( I, I )
00327             WI( I ) = ZERO
00328    10    CONTINUE
00329          DO 20 I = IHI + 1, N
00330             WR( I ) = H( I, I )
00331             WI( I ) = ZERO
00332    20    CONTINUE
00333 *
00334 *        ==== Initialize Z, if requested ====
00335 *
00336          IF( INITZ )
00337      $      CALL DLASET( 'A', N, N, ZERO, ONE, Z, LDZ )
00338 *
00339 *        ==== Quick return if possible ====
00340 *
00341          IF( ILO.EQ.IHI ) THEN
00342             WR( ILO ) = H( ILO, ILO )
00343             WI( ILO ) = ZERO
00344             RETURN
00345          END IF
00346 *
00347 *        ==== DLAHQR/DLAQR0 crossover point ====
00348 *
00349          NMIN = ILAENV( 12, 'DHSEQR', JOB( : 1 ) // COMPZ( : 1 ), N,
00350      $          ILO, IHI, LWORK )
00351          NMIN = MAX( NTINY, NMIN )
00352 *
00353 *        ==== DLAQR0 for big matrices; DLAHQR for small ones ====
00354 *
00355          IF( N.GT.NMIN ) THEN
00356             CALL DLAQR0( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO,
00357      $                   IHI, Z, LDZ, WORK, LWORK, INFO )
00358          ELSE
00359 *
00360 *           ==== Small matrix ====
00361 *
00362             CALL DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO,
00363      $                   IHI, Z, LDZ, INFO )
00364 *
00365             IF( INFO.GT.0 ) THEN
00366 *
00367 *              ==== A rare DLAHQR failure!  DLAQR0 sometimes succeeds
00368 *              .    when DLAHQR fails. ====
00369 *
00370                KBOT = INFO
00371 *
00372                IF( N.GE.NL ) THEN
00373 *
00374 *                 ==== Larger matrices have enough subdiagonal scratch
00375 *                 .    space to call DLAQR0 directly. ====
00376 *
00377                   CALL DLAQR0( WANTT, WANTZ, N, ILO, KBOT, H, LDH, WR,
00378      $                         WI, ILO, IHI, Z, LDZ, WORK, LWORK, INFO )
00379 *
00380                ELSE
00381 *
00382 *                 ==== Tiny matrices don't have enough subdiagonal
00383 *                 .    scratch space to benefit from DLAQR0.  Hence,
00384 *                 .    tiny matrices must be copied into a larger
00385 *                 .    array before calling DLAQR0. ====
00386 *
00387                   CALL DLACPY( 'A', N, N, H, LDH, HL, NL )
00388                   HL( N+1, N ) = ZERO
00389                   CALL DLASET( 'A', NL, NL-N, ZERO, ZERO, HL( 1, N+1 ),
00390      $                         NL )
00391                   CALL DLAQR0( WANTT, WANTZ, NL, ILO, KBOT, HL, NL, WR,
00392      $                         WI, ILO, IHI, Z, LDZ, WORKL, NL, INFO )
00393                   IF( WANTT .OR. INFO.NE.0 )
00394      $               CALL DLACPY( 'A', N, N, HL, NL, H, LDH )
00395                END IF
00396             END IF
00397          END IF
00398 *
00399 *        ==== Clear out the trash, if necessary. ====
00400 *
00401          IF( ( WANTT .OR. INFO.NE.0 ) .AND. N.GT.2 )
00402      $      CALL DLASET( 'L', N-2, N-2, ZERO, ZERO, H( 3, 1 ), LDH )
00403 *
00404 *        ==== Ensure reported workspace size is backward-compatible with
00405 *        .    previous LAPACK versions. ====
00406 *
00407          WORK( 1 ) = MAX( DBLE( MAX( 1, N ) ), WORK( 1 ) )
00408       END IF
00409 *
00410 *     ==== End of DHSEQR ====
00411 *
00412       END
 All Files Functions