SRC/dhseqr.f

   1 *> \brief \b DHSEQR
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DHSEQR + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dhseqr.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dhseqr.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dhseqr.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DHSEQR( JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z,
  22 *                          LDZ, WORK, LWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            IHI, ILO, INFO, LDH, LDZ, LWORK, N
  26 *       CHARACTER          COMPZ, JOB
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   H( LDH, * ), WI( * ), WORK( * ), WR( * ),
  30 *      $                   Z( LDZ, * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *>    DHSEQR computes the eigenvalues of a Hessenberg matrix H
  40 *>    and, optionally, the matrices T and Z from the Schur decomposition
  41 *>    H = Z T Z**T, where T is an upper quasi-triangular matrix (the
  42 *>    Schur form), and Z is the orthogonal matrix of Schur vectors.
  43 *>
  44 *>    Optionally Z may be postmultiplied into an input orthogonal
  45 *>    matrix Q so that this routine can give the Schur factorization
  46 *>    of a matrix A which has been reduced to the Hessenberg form H
  47 *>    by the orthogonal matrix Q:  A = Q*H*Q**T = (QZ)*T*(QZ)**T.
  48 *> \endverbatim
  49 *
  50 *  Arguments:
  51 *  ==========
  52 *
  53 *> \param[in] JOB
  54 *> \verbatim
  55 *>          JOB is CHARACTER*1
  56 *>           = 'E':  compute eigenvalues only;
  57 *>           = 'S':  compute eigenvalues and the Schur form T.
  58 *> \endverbatim
  59 *>
  60 *> \param[in] COMPZ
  61 *> \verbatim
  62 *>          COMPZ is CHARACTER*1
  63 *>           = 'N':  no Schur vectors are computed;
  64 *>           = 'I':  Z is initialized to the unit matrix and the matrix Z
  65 *>                   of Schur vectors of H is returned;
  66 *>           = 'V':  Z must contain an orthogonal matrix Q on entry, and
  67 *>                   the product Q*Z is returned.
  68 *> \endverbatim
  69 *>
  70 *> \param[in] N
  71 *> \verbatim
  72 *>          N is INTEGER
  73 *>           The order of the matrix H.  N .GE. 0.
  74 *> \endverbatim
  75 *>
  76 *> \param[in] ILO
  77 *> \verbatim
  78 *>          ILO is INTEGER
  79 *> \endverbatim
  80 *>
  81 *> \param[in] IHI
  82 *> \verbatim
  83 *>          IHI is INTEGER
  84 *>
  85 *>           It is assumed that H is already upper triangular in rows
  86 *>           and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
  87 *>           set by a previous call to DGEBAL, and then passed to ZGEHRD
  88 *>           when the matrix output by DGEBAL is reduced to Hessenberg
  89 *>           form. Otherwise ILO and IHI should be set to 1 and N
  90 *>           respectively.  If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
  91 *>           If N = 0, then ILO = 1 and IHI = 0.
  92 *> \endverbatim
  93 *>
  94 *> \param[in,out] H
  95 *> \verbatim
  96 *>          H is DOUBLE PRECISION array, dimension (LDH,N)
  97 *>           On entry, the upper Hessenberg matrix H.
  98 *>           On exit, if INFO = 0 and JOB = 'S', then H contains the
  99 *>           upper quasi-triangular matrix T from the Schur decomposition
 100 *>           (the Schur form); 2-by-2 diagonal blocks (corresponding to
 101 *>           complex conjugate pairs of eigenvalues) are returned in
 102 *>           standard form, with H(i,i) = H(i+1,i+1) and
 103 *>           H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and JOB = 'E', the
 104 *>           contents of H are unspecified on exit.  (The output value of
 105 *>           H when INFO.GT.0 is given under the description of INFO
 106 *>           below.)
 107 *>
 108 *>           Unlike earlier versions of DHSEQR, this subroutine may
 109 *>           explicitly H(i,j) = 0 for i.GT.j and j = 1, 2, ... ILO-1
 110 *>           or j = IHI+1, IHI+2, ... N.
 111 *> \endverbatim
 112 *>
 113 *> \param[in] LDH
 114 *> \verbatim
 115 *>          LDH is INTEGER
 116 *>           The leading dimension of the array H. LDH .GE. max(1,N).
 117 *> \endverbatim
 118 *>
 119 *> \param[out] WR
 120 *> \verbatim
 121 *>          WR is DOUBLE PRECISION array, dimension (N)
 122 *> \endverbatim
 123 *>
 124 *> \param[out] WI
 125 *> \verbatim
 126 *>          WI is DOUBLE PRECISION array, dimension (N)
 127 *>
 128 *>           The real and imaginary parts, respectively, of the computed
 129 *>           eigenvalues. If two eigenvalues are computed as a complex
 130 *>           conjugate pair, they are stored in consecutive elements of
 131 *>           WR and WI, say the i-th and (i+1)th, with WI(i) .GT. 0 and
 132 *>           WI(i+1) .LT. 0. If JOB = 'S', the eigenvalues are stored in
 133 *>           the same order as on the diagonal of the Schur form returned
 134 *>           in H, with WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2
 135 *>           diagonal block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and
 136 *>           WI(i+1) = -WI(i).
 137 *> \endverbatim
 138 *>
 139 *> \param[in,out] Z
 140 *> \verbatim
 141 *>          Z is DOUBLE PRECISION array, dimension (LDZ,N)
 142 *>           If COMPZ = 'N', Z is not referenced.
 143 *>           If COMPZ = 'I', on entry Z need not be set and on exit,
 144 *>           if INFO = 0, Z contains the orthogonal matrix Z of the Schur
 145 *>           vectors of H.  If COMPZ = 'V', on entry Z must contain an
 146 *>           N-by-N matrix Q, which is assumed to be equal to the unit
 147 *>           matrix except for the submatrix Z(ILO:IHI,ILO:IHI). On exit,
 148 *>           if INFO = 0, Z contains Q*Z.
 149 *>           Normally Q is the orthogonal matrix generated by DORGHR
 150 *>           after the call to DGEHRD which formed the Hessenberg matrix
 151 *>           H. (The output value of Z when INFO.GT.0 is given under
 152 *>           the description of INFO below.)
 153 *> \endverbatim
 154 *>
 155 *> \param[in] LDZ
 156 *> \verbatim
 157 *>          LDZ is INTEGER
 158 *>           The leading dimension of the array Z.  if COMPZ = 'I' or
 159 *>           COMPZ = 'V', then LDZ.GE.MAX(1,N).  Otherwize, LDZ.GE.1.
 160 *> \endverbatim
 161 *>
 162 *> \param[out] WORK
 163 *> \verbatim
 164 *>          WORK is DOUBLE PRECISION array, dimension (LWORK)
 165 *>           On exit, if INFO = 0, WORK(1) returns an estimate of
 166 *>           the optimal value for LWORK.
 167 *> \endverbatim
 168 *>
 169 *> \param[in] LWORK
 170 *> \verbatim
 171 *>          LWORK is INTEGER
 172 *>           The dimension of the array WORK.  LWORK .GE. max(1,N)
 173 *>           is sufficient and delivers very good and sometimes
 174 *>           optimal performance.  However, LWORK as large as 11*N
 175 *>           may be required for optimal performance.  A workspace
 176 *>           query is recommended to determine the optimal workspace
 177 *>           size.
 178 *>
 179 *>           If LWORK = -1, then DHSEQR does a workspace query.
 180 *>           In this case, DHSEQR checks the input parameters and
 181 *>           estimates the optimal workspace size for the given
 182 *>           values of N, ILO and IHI.  The estimate is returned
 183 *>           in WORK(1).  No error message related to LWORK is
 184 *>           issued by XERBLA.  Neither H nor Z are accessed.
 185 *> \endverbatim
 186 *>
 187 *> \param[out] INFO
 188 *> \verbatim
 189 *>          INFO is INTEGER
 190 *>             =  0:  successful exit
 191 *>           .LT. 0:  if INFO = -i, the i-th argument had an illegal
 192 *>                    value
 193 *>           .GT. 0:  if INFO = i, DHSEQR failed to compute all of
 194 *>                the eigenvalues.  Elements 1:ilo-1 and i+1:n of WR
 195 *>                and WI contain those eigenvalues which have been
 196 *>                successfully computed.  (Failures are rare.)
 197 *>
 198 *>                If INFO .GT. 0 and JOB = 'E', then on exit, the
 199 *>                remaining unconverged eigenvalues are the eigen-
 200 *>                values of the upper Hessenberg matrix rows and
 201 *>                columns ILO through INFO of the final, output
 202 *>                value of H.
 203 *>
 204 *>                If INFO .GT. 0 and JOB   = 'S', then on exit
 205 *>
 206 *>           (*)  (initial value of H)*U  = U*(final value of H)
 207 *>
 208 *>                where U is an orthogonal matrix.  The final
 209 *>                value of H is upper Hessenberg and quasi-triangular
 210 *>                in rows and columns INFO+1 through IHI.
 211 *>
 212 *>                If INFO .GT. 0 and COMPZ = 'V', then on exit
 213 *>
 214 *>                  (final value of Z)  =  (initial value of Z)*U
 215 *>
 216 *>                where U is the orthogonal matrix in (*) (regard-
 217 *>                less of the value of JOB.)
 218 *>
 219 *>                If INFO .GT. 0 and COMPZ = 'I', then on exit
 220 *>                      (final value of Z)  = U
 221 *>                where U is the orthogonal matrix in (*) (regard-
 222 *>                less of the value of JOB.)
 223 *>
 224 *>                If INFO .GT. 0 and COMPZ = 'N', then Z is not
 225 *>                accessed.
 226 *> \endverbatim
 227 *
 228 *  Authors:
 229 *  ========
 230 *
 231 *> \author Univ. of Tennessee
 232 *> \author Univ. of California Berkeley
 233 *> \author Univ. of Colorado Denver
 234 *> \author NAG Ltd.
 235 *
 236 *> \date November 2011
 237 *
 238 *> \ingroup doubleOTHERcomputational
 239 *
 240 *> \par Contributors:
 241 *  ==================
 242 *>
 243 *>       Karen Braman and Ralph Byers, Department of Mathematics,
 244 *>       University of Kansas, USA
 245 *
 246 *> \par Further Details:
 247 *  =====================
 248 *>
 249 *> \verbatim
 250 *>
 251 *>             Default values supplied by
 252 *>             ILAENV(ISPEC,'DHSEQR',JOB(:1)//COMPZ(:1),N,ILO,IHI,LWORK).
 253 *>             It is suggested that these defaults be adjusted in order
 254 *>             to attain best performance in each particular
 255 *>             computational environment.
 256 *>
 257 *>            ISPEC=12: The DLAHQR vs DLAQR0 crossover point.
 258 *>                      Default: 75. (Must be at least 11.)
 259 *>
 260 *>            ISPEC=13: Recommended deflation window size.
 261 *>                      This depends on ILO, IHI and NS.  NS is the
 262 *>                      number of simultaneous shifts returned
 263 *>                      by ILAENV(ISPEC=15).  (See ISPEC=15 below.)
 264 *>                      The default for (IHI-ILO+1).LE.500 is NS.
 265 *>                      The default for (IHI-ILO+1).GT.500 is 3*NS/2.
 266 *>
 267 *>            ISPEC=14: Nibble crossover point. (See IPARMQ for
 268 *>                      details.)  Default: 14% of deflation window
 269 *>                      size.
 270 *>
 271 *>            ISPEC=15: Number of simultaneous shifts in a multishift
 272 *>                      QR iteration.
 273 *>
 274 *>                      If IHI-ILO+1 is ...
 275 *>
 276 *>                      greater than      ...but less    ... the
 277 *>                      or equal to ...      than        default is
 278 *>
 279 *>                           1               30          NS =   2(+)
 280 *>                          30               60          NS =   4(+)
 281 *>                          60              150          NS =  10(+)
 282 *>                         150              590          NS =  **
 283 *>                         590             3000          NS =  64
 284 *>                        3000             6000          NS = 128
 285 *>                        6000             infinity      NS = 256
 286 *>
 287 *>                  (+)  By default some or all matrices of this order
 288 *>                       are passed to the implicit double shift routine
 289 *>                       DLAHQR and this parameter is ignored.  See
 290 *>                       ISPEC=12 above and comments in IPARMQ for
 291 *>                       details.
 292 *>
 293 *>                 (**)  The asterisks (**) indicate an ad-hoc
 294 *>                       function of N increasing from 10 to 64.
 295 *>
 296 *>            ISPEC=16: Select structured matrix multiply.
 297 *>                      If the number of simultaneous shifts (specified
 298 *>                      by ISPEC=15) is less than 14, then the default
 299 *>                      for ISPEC=16 is 0.  Otherwise the default for
 300 *>                      ISPEC=16 is 2.
 301 *> \endverbatim
 302 *
 303 *> \par References:
 304 *  ================
 305 *>
 306 *>       K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
 307 *>       Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
 308 *>       Performance, SIAM Journal of Matrix Analysis, volume 23, pages
 309 *>       929--947, 2002.
 310 *> \n
 311 *>       K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
 312 *>       Algorithm Part II: Aggressive Early Deflation, SIAM Journal
 313 *>       of Matrix Analysis, volume 23, pages 948--973, 2002.
 314 *
 315 *  =====================================================================
 316       SUBROUTINE DHSEQR( JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z,
 317      $                   LDZ, WORK, LWORK, INFO )
 318 *
 319 *  -- LAPACK computational routine (version 3.4.0) --
 320 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 321 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 322 *     November 2011
 323 *
 324 *     .. Scalar Arguments ..
 325       INTEGER            IHI, ILO, INFO, LDH, LDZ, LWORK, N
 326       CHARACTER          COMPZ, JOB
 327 *     ..
 328 *     .. Array Arguments ..
 329       DOUBLE PRECISION   H( LDH, * ), WI( * ), WORK( * ), WR( * ),
 330      $                   Z( LDZ, * )
 331 *     ..
 332 *
 333 *  =====================================================================
 334 *
 335 *     .. Parameters ..
 336 *
 337 *     ==== Matrices of order NTINY or smaller must be processed by
 338 *     .    DLAHQR because of insufficient subdiagonal scratch space.
 339 *     .    (This is a hard limit.) ====
 340       INTEGER            NTINY
 341       PARAMETER          ( NTINY = 11 )
 342 *
 343 *     ==== NL allocates some local workspace to help small matrices
 344 *     .    through a rare DLAHQR failure.  NL .GT. NTINY = 11 is
 345 *     .    required and NL .LE. NMIN = ILAENV(ISPEC=12,...) is recom-
 346 *     .    mended.  (The default value of NMIN is 75.)  Using NL = 49
 347 *     .    allows up to six simultaneous shifts and a 16-by-16
 348 *     .    deflation window.  ====
 349       INTEGER            NL
 350       PARAMETER          ( NL = 49 )
 351       DOUBLE PRECISION   ZERO, ONE
 352       PARAMETER          ( ZERO = 0.0d0, ONE = 1.0d0 )
 353 *     ..
 354 *     .. Local Arrays ..
 355       DOUBLE PRECISION   HL( NL, NL ), WORKL( NL )
 356 *     ..
 357 *     .. Local Scalars ..
 358       INTEGER            I, KBOT, NMIN
 359       LOGICAL            INITZ, LQUERY, WANTT, WANTZ
 360 *     ..
 361 *     .. External Functions ..
 362       INTEGER            ILAENV
 363       LOGICAL            LSAME
 364       EXTERNAL           ILAENV, LSAME
 365 *     ..
 366 *     .. External Subroutines ..
 367       EXTERNAL           DLACPY, DLAHQR, DLAQR0, DLASET, XERBLA
 368 *     ..
 369 *     .. Intrinsic Functions ..
 370       INTRINSIC          DBLE, MAX, MIN
 371 *     ..
 372 *     .. Executable Statements ..
 373 *
 374 *     ==== Decode and check the input parameters. ====
 375 *
 376       WANTT = LSAME( JOB, 'S' )
 377       INITZ = LSAME( COMPZ, 'I' )
 378       WANTZ = INITZ .OR. LSAME( COMPZ, 'V' )
 379       WORK( 1 ) = DBLE( MAX( 1, N ) )
 380       LQUERY = LWORK.EQ.-1
 381 *
 382       INFO = 0
 383       IF( .NOT.LSAME( JOB, 'E' ) .AND. .NOT.WANTT ) THEN
 384          INFO = -1
 385       ELSE IF( .NOT.LSAME( COMPZ, 'N' ) .AND. .NOT.WANTZ ) THEN
 386          INFO = -2
 387       ELSE IF( N.LT.0 ) THEN
 388          INFO = -3
 389       ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN
 390          INFO = -4
 391       ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN
 392          INFO = -5
 393       ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
 394          INFO = -7
 395       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.MAX( 1, N ) ) ) THEN
 396          INFO = -11
 397       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
 398          INFO = -13
 399       END IF
 400 *
 401       IF( INFO.NE.0 ) THEN
 402 *
 403 *        ==== Quick return in case of invalid argument. ====
 404 *
 405          CALL XERBLA( 'DHSEQR', -INFO )
 406          RETURN
 407 *
 408       ELSE IF( N.EQ.0 ) THEN
 409 *
 410 *        ==== Quick return in case N = 0; nothing to do. ====
 411 *
 412          RETURN
 413 *
 414       ELSE IF( LQUERY ) THEN
 415 *
 416 *        ==== Quick return in case of a workspace query ====
 417 *
 418          CALL DLAQR0( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO,
 419      $                IHI, Z, LDZ, WORK, LWORK, INFO )
 420 *        ==== Ensure reported workspace size is backward-compatible with
 421 *        .    previous LAPACK versions. ====
 422          WORK( 1 ) = MAX( DBLE( MAX( 1, N ) ), WORK( 1 ) )
 423          RETURN
 424 *
 425       ELSE
 426 *
 427 *        ==== copy eigenvalues isolated by DGEBAL ====
 428 *
 429          DO 10 I = 1, ILO - 1
 430             WR( I ) = H( I, I )
 431             WI( I ) = ZERO
 432    10    CONTINUE
 433          DO 20 I = IHI + 1, N
 434             WR( I ) = H( I, I )
 435             WI( I ) = ZERO
 436    20    CONTINUE
 437 *
 438 *        ==== Initialize Z, if requested ====
 439 *
 440          IF( INITZ )
 441      $      CALL DLASET( 'A', N, N, ZERO, ONE, Z, LDZ )
 442 *
 443 *        ==== Quick return if possible ====
 444 *
 445          IF( ILO.EQ.IHI ) THEN
 446             WR( ILO ) = H( ILO, ILO )
 447             WI( ILO ) = ZERO
 448             RETURN
 449          END IF
 450 *
 451 *        ==== DLAHQR/DLAQR0 crossover point ====
 452 *
 453          NMIN = ILAENV( 12, 'DHSEQR', JOB( : 1 ) // COMPZ( : 1 ), N,
 454      $          ILO, IHI, LWORK )
 455          NMIN = MAX( NTINY, NMIN )
 456 *
 457 *        ==== DLAQR0 for big matrices; DLAHQR for small ones ====
 458 *
 459          IF( N.GT.NMIN ) THEN
 460             CALL DLAQR0( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO,
 461      $                   IHI, Z, LDZ, WORK, LWORK, INFO )
 462          ELSE
 463 *
 464 *           ==== Small matrix ====
 465 *
 466             CALL DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO,
 467      $                   IHI, Z, LDZ, INFO )
 468 *
 469             IF( INFO.GT.0 ) THEN
 470 *
 471 *              ==== A rare DLAHQR failure!  DLAQR0 sometimes succeeds
 472 *              .    when DLAHQR fails. ====
 473 *
 474                KBOT = INFO
 475 *
 476                IF( N.GE.NL ) THEN
 477 *
 478 *                 ==== Larger matrices have enough subdiagonal scratch
 479 *                 .    space to call DLAQR0 directly. ====
 480 *
 481                   CALL DLAQR0( WANTT, WANTZ, N, ILO, KBOT, H, LDH, WR,
 482      $                         WI, ILO, IHI, Z, LDZ, WORK, LWORK, INFO )
 483 *
 484                ELSE
 485 *
 486 *                 ==== Tiny matrices don't have enough subdiagonal
 487 *                 .    scratch space to benefit from DLAQR0.  Hence,
 488 *                 .    tiny matrices must be copied into a larger
 489 *                 .    array before calling DLAQR0. ====
 490 *
 491                   CALL DLACPY( 'A', N, N, H, LDH, HL, NL )
 492                   HL( N+1, N ) = ZERO
 493                   CALL DLASET( 'A', NL, NL-N, ZERO, ZERO, HL( 1, N+1 ),
 494      $                         NL )
 495                   CALL DLAQR0( WANTT, WANTZ, NL, ILO, KBOT, HL, NL, WR,
 496      $                         WI, ILO, IHI, Z, LDZ, WORKL, NL, INFO )
 497                   IF( WANTT .OR. INFO.NE.0 )
 498      $               CALL DLACPY( 'A', N, N, HL, NL, H, LDH )
 499                END IF
 500             END IF
 501          END IF
 502 *
 503 *        ==== Clear out the trash, if necessary. ====
 504 *
 505          IF( ( WANTT .OR. INFO.NE.0 ) .AND. N.GT.2 )
 506      $      CALL DLASET( 'L', N-2, N-2, ZERO, ZERO, H( 3, 1 ), LDH )
 507 *
 508 *        ==== Ensure reported workspace size is backward-compatible with
 509 *        .    previous LAPACK versions. ====
 510 *
 511          WORK( 1 ) = MAX( DBLE( MAX( 1, N ) ), WORK( 1 ) )
 512       END IF
 513 *
 514 *     ==== End of DHSEQR ====
 515 *
 516       END