1*> \brief \b SHSEQR 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download SHSEQR + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/shseqr.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/shseqr.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/shseqr.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SHSEQR( 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* REAL H( LDH, * ), WI( * ), WORK( * ), WR( * ), 30* $ Z( LDZ, * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> SHSEQR 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 >= 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 SGEBAL, and then passed to ZGEHRD 88*> when the matrix output by SGEBAL is reduced to Hessenberg 89*> form. Otherwise ILO and IHI should be set to 1 and N 90*> respectively. If N > 0, then 1 <= ILO <= IHI <= N. 91*> If N = 0, then ILO = 1 and IHI = 0. 92*> \endverbatim 93*> 94*> \param[in,out] H 95*> \verbatim 96*> H is REAL 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) < 0. If INFO = 0 and JOB = 'E', the 104*> contents of H are unspecified on exit. (The output value of 105*> H when INFO > 0 is given under the description of INFO 106*> below.) 107*> 108*> Unlike earlier versions of SHSEQR, this subroutine may 109*> explicitly H(i,j) = 0 for i > 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 >= max(1,N). 117*> \endverbatim 118*> 119*> \param[out] WR 120*> \verbatim 121*> WR is REAL array, dimension (N) 122*> \endverbatim 123*> 124*> \param[out] WI 125*> \verbatim 126*> WI is REAL 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) > 0 and 132*> WI(i+1) < 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 REAL 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 SORGHR 150*> after the call to SGEHRD which formed the Hessenberg matrix 151*> H. (The output value of Z when INFO > 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 >= MAX(1,N). Otherwise, LDZ >= 1. 160*> \endverbatim 161*> 162*> \param[out] WORK 163*> \verbatim 164*> WORK is REAL 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 >= 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 SHSEQR does a workspace query. 180*> In this case, SHSEQR 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*> < 0: if INFO = -i, the i-th argument had an illegal 192*> value 193*> > 0: if INFO = i, SHSEQR 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 > 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 > 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 > 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 > 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 > 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*> \ingroup realOTHERcomputational 237* 238*> \par Contributors: 239* ================== 240*> 241*> Karen Braman and Ralph Byers, Department of Mathematics, 242*> University of Kansas, USA 243* 244*> \par Further Details: 245* ===================== 246*> 247*> \verbatim 248*> 249*> Default values supplied by 250*> ILAENV(ISPEC,'SHSEQR',JOB(:1)//COMPZ(:1),N,ILO,IHI,LWORK). 251*> It is suggested that these defaults be adjusted in order 252*> to attain best performance in each particular 253*> computational environment. 254*> 255*> ISPEC=12: The SLAHQR vs SLAQR0 crossover point. 256*> Default: 75. (Must be at least 11.) 257*> 258*> ISPEC=13: Recommended deflation window size. 259*> This depends on ILO, IHI and NS. NS is the 260*> number of simultaneous shifts returned 261*> by ILAENV(ISPEC=15). (See ISPEC=15 below.) 262*> The default for (IHI-ILO+1) <= 500 is NS. 263*> The default for (IHI-ILO+1) > 500 is 3*NS/2. 264*> 265*> ISPEC=14: Nibble crossover point. (See IPARMQ for 266*> details.) Default: 14% of deflation window 267*> size. 268*> 269*> ISPEC=15: Number of simultaneous shifts in a multishift 270*> QR iteration. 271*> 272*> If IHI-ILO+1 is ... 273*> 274*> greater than ...but less ... the 275*> or equal to ... than default is 276*> 277*> 1 30 NS = 2(+) 278*> 30 60 NS = 4(+) 279*> 60 150 NS = 10(+) 280*> 150 590 NS = ** 281*> 590 3000 NS = 64 282*> 3000 6000 NS = 128 283*> 6000 infinity NS = 256 284*> 285*> (+) By default some or all matrices of this order 286*> are passed to the implicit double shift routine 287*> SLAHQR and this parameter is ignored. See 288*> ISPEC=12 above and comments in IPARMQ for 289*> details. 290*> 291*> (**) The asterisks (**) indicate an ad-hoc 292*> function of N increasing from 10 to 64. 293*> 294*> ISPEC=16: Select structured matrix multiply. 295*> If the number of simultaneous shifts (specified 296*> by ISPEC=15) is less than 14, then the default 297*> for ISPEC=16 is 0. Otherwise the default for 298*> ISPEC=16 is 2. 299*> \endverbatim 300* 301*> \par References: 302* ================ 303*> 304*> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR 305*> Algorithm Part I: Maintaining Well Focused Shifts, and Level 3 306*> Performance, SIAM Journal of Matrix Analysis, volume 23, pages 307*> 929--947, 2002. 308*> \n 309*> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR 310*> Algorithm Part II: Aggressive Early Deflation, SIAM Journal 311*> of Matrix Analysis, volume 23, pages 948--973, 2002. 312* 313* ===================================================================== 314 SUBROUTINE SHSEQR( JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, 315 $ LDZ, WORK, LWORK, INFO ) 316* 317* -- LAPACK computational routine -- 318* -- LAPACK is a software package provided by Univ. of Tennessee, -- 319* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 320* 321* .. Scalar Arguments .. 322 INTEGER IHI, ILO, INFO, LDH, LDZ, LWORK, N 323 CHARACTER COMPZ, JOB 324* .. 325* .. Array Arguments .. 326 REAL H( LDH, * ), WI( * ), WORK( * ), WR( * ), 327 $ Z( LDZ, * ) 328* .. 329* 330* ===================================================================== 331* 332* .. Parameters .. 333* 334* ==== Matrices of order NTINY or smaller must be processed by 335* . SLAHQR because of insufficient subdiagonal scratch space. 336* . (This is a hard limit.) ==== 337 INTEGER NTINY 338 PARAMETER ( NTINY = 15 ) 339* 340* ==== NL allocates some local workspace to help small matrices 341* . through a rare SLAHQR failure. NL > NTINY = 15 is 342* . required and NL <= NMIN = ILAENV(ISPEC=12,...) is recom- 343* . mended. (The default value of NMIN is 75.) Using NL = 49 344* . allows up to six simultaneous shifts and a 16-by-16 345* . deflation window. ==== 346 INTEGER NL 347 PARAMETER ( NL = 49 ) 348 REAL ZERO, ONE 349 PARAMETER ( ZERO = 0.0e0, ONE = 1.0e0 ) 350* .. 351* .. Local Arrays .. 352 REAL HL( NL, NL ), WORKL( NL ) 353* .. 354* .. Local Scalars .. 355 INTEGER I, KBOT, NMIN 356 LOGICAL INITZ, LQUERY, WANTT, WANTZ 357* .. 358* .. External Functions .. 359 INTEGER ILAENV 360 LOGICAL LSAME 361 EXTERNAL ILAENV, LSAME 362* .. 363* .. External Subroutines .. 364 EXTERNAL SLACPY, SLAHQR, SLAQR0, SLASET, XERBLA 365* .. 366* .. Intrinsic Functions .. 367 INTRINSIC MAX, MIN, REAL 368* .. 369* .. Executable Statements .. 370* 371* ==== Decode and check the input parameters. ==== 372* 373 WANTT = LSAME( JOB, 'S' ) 374 INITZ = LSAME( COMPZ, 'I' ) 375 WANTZ = INITZ .OR. LSAME( COMPZ, 'V' ) 376 WORK( 1 ) = REAL( MAX( 1, N ) ) 377 LQUERY = LWORK.EQ.-1 378* 379 INFO = 0 380 IF( .NOT.LSAME( JOB, 'E' ) .AND. .NOT.WANTT ) THEN 381 INFO = -1 382 ELSE IF( .NOT.LSAME( COMPZ, 'N' ) .AND. .NOT.WANTZ ) THEN 383 INFO = -2 384 ELSE IF( N.LT.0 ) THEN 385 INFO = -3 386 ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN 387 INFO = -4 388 ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN 389 INFO = -5 390 ELSE IF( LDH.LT.MAX( 1, N ) ) THEN 391 INFO = -7 392 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.MAX( 1, N ) ) ) THEN 393 INFO = -11 394 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN 395 INFO = -13 396 END IF 397* 398 IF( INFO.NE.0 ) THEN 399* 400* ==== Quick return in case of invalid argument. ==== 401* 402 CALL XERBLA( 'SHSEQR', -INFO ) 403 RETURN 404* 405 ELSE IF( N.EQ.0 ) THEN 406* 407* ==== Quick return in case N = 0; nothing to do. ==== 408* 409 RETURN 410* 411 ELSE IF( LQUERY ) THEN 412* 413* ==== Quick return in case of a workspace query ==== 414* 415 CALL SLAQR0( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO, 416 $ IHI, Z, LDZ, WORK, LWORK, INFO ) 417* ==== Ensure reported workspace size is backward-compatible with 418* . previous LAPACK versions. ==== 419 WORK( 1 ) = MAX( REAL( MAX( 1, N ) ), WORK( 1 ) ) 420 RETURN 421* 422 ELSE 423* 424* ==== copy eigenvalues isolated by SGEBAL ==== 425* 426 DO 10 I = 1, ILO - 1 427 WR( I ) = H( I, I ) 428 WI( I ) = ZERO 429 10 CONTINUE 430 DO 20 I = IHI + 1, N 431 WR( I ) = H( I, I ) 432 WI( I ) = ZERO 433 20 CONTINUE 434* 435* ==== Initialize Z, if requested ==== 436* 437 IF( INITZ ) 438 $ CALL SLASET( 'A', N, N, ZERO, ONE, Z, LDZ ) 439* 440* ==== Quick return if possible ==== 441* 442 IF( ILO.EQ.IHI ) THEN 443 WR( ILO ) = H( ILO, ILO ) 444 WI( ILO ) = ZERO 445 RETURN 446 END IF 447* 448* ==== SLAHQR/SLAQR0 crossover point ==== 449* 450 NMIN = ILAENV( 12, 'SHSEQR', JOB( : 1 ) // COMPZ( : 1 ), N, 451 $ ILO, IHI, LWORK ) 452 NMIN = MAX( NTINY, NMIN ) 453* 454* ==== SLAQR0 for big matrices; SLAHQR for small ones ==== 455* 456 IF( N.GT.NMIN ) THEN 457 CALL SLAQR0( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO, 458 $ IHI, Z, LDZ, WORK, LWORK, INFO ) 459 ELSE 460* 461* ==== Small matrix ==== 462* 463 CALL SLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO, 464 $ IHI, Z, LDZ, INFO ) 465* 466 IF( INFO.GT.0 ) THEN 467* 468* ==== A rare SLAHQR failure! SLAQR0 sometimes succeeds 469* . when SLAHQR fails. ==== 470* 471 KBOT = INFO 472* 473 IF( N.GE.NL ) THEN 474* 475* ==== Larger matrices have enough subdiagonal scratch 476* . space to call SLAQR0 directly. ==== 477* 478 CALL SLAQR0( WANTT, WANTZ, N, ILO, KBOT, H, LDH, WR, 479 $ WI, ILO, IHI, Z, LDZ, WORK, LWORK, INFO ) 480* 481 ELSE 482* 483* ==== Tiny matrices don't have enough subdiagonal 484* . scratch space to benefit from SLAQR0. Hence, 485* . tiny matrices must be copied into a larger 486* . array before calling SLAQR0. ==== 487* 488 CALL SLACPY( 'A', N, N, H, LDH, HL, NL ) 489 HL( N+1, N ) = ZERO 490 CALL SLASET( 'A', NL, NL-N, ZERO, ZERO, HL( 1, N+1 ), 491 $ NL ) 492 CALL SLAQR0( WANTT, WANTZ, NL, ILO, KBOT, HL, NL, WR, 493 $ WI, ILO, IHI, Z, LDZ, WORKL, NL, INFO ) 494 IF( WANTT .OR. INFO.NE.0 ) 495 $ CALL SLACPY( 'A', N, N, HL, NL, H, LDH ) 496 END IF 497 END IF 498 END IF 499* 500* ==== Clear out the trash, if necessary. ==== 501* 502 IF( ( WANTT .OR. INFO.NE.0 ) .AND. N.GT.2 ) 503 $ CALL SLASET( 'L', N-2, N-2, ZERO, ZERO, H( 3, 1 ), LDH ) 504* 505* ==== Ensure reported workspace size is backward-compatible with 506* . previous LAPACK versions. ==== 507* 508 WORK( 1 ) = MAX( REAL( MAX( 1, N ) ), WORK( 1 ) ) 509 END IF 510* 511* ==== End of SHSEQR ==== 512* 513 END 514