1*> \brief \b ZLAQR2 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation). 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZLAQR2 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaqr2.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaqr2.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaqr2.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZLAQR2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, 22* IHIZ, Z, LDZ, NS, ND, SH, V, LDV, NH, T, LDT, 23* NV, WV, LDWV, WORK, LWORK ) 24* 25* .. Scalar Arguments .. 26* INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV, 27* $ LDZ, LWORK, N, ND, NH, NS, NV, NW 28* LOGICAL WANTT, WANTZ 29* .. 30* .. Array Arguments .. 31* COMPLEX*16 H( LDH, * ), SH( * ), T( LDT, * ), V( LDV, * ), 32* $ WORK( * ), WV( LDWV, * ), Z( LDZ, * ) 33* .. 34* 35* 36*> \par Purpose: 37* ============= 38*> 39*> \verbatim 40*> 41*> ZLAQR2 is identical to ZLAQR3 except that it avoids 42*> recursion by calling ZLAHQR instead of ZLAQR4. 43*> 44*> Aggressive early deflation: 45*> 46*> ZLAQR2 accepts as input an upper Hessenberg matrix 47*> H and performs an unitary similarity transformation 48*> designed to detect and deflate fully converged eigenvalues from 49*> a trailing principal submatrix. On output H has been over- 50*> written by a new Hessenberg matrix that is a perturbation of 51*> an unitary similarity transformation of H. It is to be 52*> hoped that the final version of H has many zero subdiagonal 53*> entries. 54*> 55*> \endverbatim 56* 57* Arguments: 58* ========== 59* 60*> \param[in] WANTT 61*> \verbatim 62*> WANTT is LOGICAL 63*> If .TRUE., then the Hessenberg matrix H is fully updated 64*> so that the triangular Schur factor may be 65*> computed (in cooperation with the calling subroutine). 66*> If .FALSE., then only enough of H is updated to preserve 67*> the eigenvalues. 68*> \endverbatim 69*> 70*> \param[in] WANTZ 71*> \verbatim 72*> WANTZ is LOGICAL 73*> If .TRUE., then the unitary matrix Z is updated so 74*> so that the unitary Schur factor may be computed 75*> (in cooperation with the calling subroutine). 76*> If .FALSE., then Z is not referenced. 77*> \endverbatim 78*> 79*> \param[in] N 80*> \verbatim 81*> N is INTEGER 82*> The order of the matrix H and (if WANTZ is .TRUE.) the 83*> order of the unitary matrix Z. 84*> \endverbatim 85*> 86*> \param[in] KTOP 87*> \verbatim 88*> KTOP is INTEGER 89*> It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0. 90*> KBOT and KTOP together determine an isolated block 91*> along the diagonal of the Hessenberg matrix. 92*> \endverbatim 93*> 94*> \param[in] KBOT 95*> \verbatim 96*> KBOT is INTEGER 97*> It is assumed without a check that either 98*> KBOT = N or H(KBOT+1,KBOT)=0. KBOT and KTOP together 99*> determine an isolated block along the diagonal of the 100*> Hessenberg matrix. 101*> \endverbatim 102*> 103*> \param[in] NW 104*> \verbatim 105*> NW is INTEGER 106*> Deflation window size. 1 <= NW <= (KBOT-KTOP+1). 107*> \endverbatim 108*> 109*> \param[in,out] H 110*> \verbatim 111*> H is COMPLEX*16 array, dimension (LDH,N) 112*> On input the initial N-by-N section of H stores the 113*> Hessenberg matrix undergoing aggressive early deflation. 114*> On output H has been transformed by a unitary 115*> similarity transformation, perturbed, and the returned 116*> to Hessenberg form that (it is to be hoped) has some 117*> zero subdiagonal entries. 118*> \endverbatim 119*> 120*> \param[in] LDH 121*> \verbatim 122*> LDH is INTEGER 123*> Leading dimension of H just as declared in the calling 124*> subroutine. N <= LDH 125*> \endverbatim 126*> 127*> \param[in] ILOZ 128*> \verbatim 129*> ILOZ is INTEGER 130*> \endverbatim 131*> 132*> \param[in] IHIZ 133*> \verbatim 134*> IHIZ is INTEGER 135*> Specify the rows of Z to which transformations must be 136*> applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N. 137*> \endverbatim 138*> 139*> \param[in,out] Z 140*> \verbatim 141*> Z is COMPLEX*16 array, dimension (LDZ,N) 142*> IF WANTZ is .TRUE., then on output, the unitary 143*> similarity transformation mentioned above has been 144*> accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right. 145*> If WANTZ is .FALSE., then Z is unreferenced. 146*> \endverbatim 147*> 148*> \param[in] LDZ 149*> \verbatim 150*> LDZ is INTEGER 151*> The leading dimension of Z just as declared in the 152*> calling subroutine. 1 <= LDZ. 153*> \endverbatim 154*> 155*> \param[out] NS 156*> \verbatim 157*> NS is INTEGER 158*> The number of unconverged (ie approximate) eigenvalues 159*> returned in SR and SI that may be used as shifts by the 160*> calling subroutine. 161*> \endverbatim 162*> 163*> \param[out] ND 164*> \verbatim 165*> ND is INTEGER 166*> The number of converged eigenvalues uncovered by this 167*> subroutine. 168*> \endverbatim 169*> 170*> \param[out] SH 171*> \verbatim 172*> SH is COMPLEX*16 array, dimension (KBOT) 173*> On output, approximate eigenvalues that may 174*> be used for shifts are stored in SH(KBOT-ND-NS+1) 175*> through SR(KBOT-ND). Converged eigenvalues are 176*> stored in SH(KBOT-ND+1) through SH(KBOT). 177*> \endverbatim 178*> 179*> \param[out] V 180*> \verbatim 181*> V is COMPLEX*16 array, dimension (LDV,NW) 182*> An NW-by-NW work array. 183*> \endverbatim 184*> 185*> \param[in] LDV 186*> \verbatim 187*> LDV is INTEGER 188*> The leading dimension of V just as declared in the 189*> calling subroutine. NW <= LDV 190*> \endverbatim 191*> 192*> \param[in] NH 193*> \verbatim 194*> NH is INTEGER 195*> The number of columns of T. NH >= NW. 196*> \endverbatim 197*> 198*> \param[out] T 199*> \verbatim 200*> T is COMPLEX*16 array, dimension (LDT,NW) 201*> \endverbatim 202*> 203*> \param[in] LDT 204*> \verbatim 205*> LDT is INTEGER 206*> The leading dimension of T just as declared in the 207*> calling subroutine. NW <= LDT 208*> \endverbatim 209*> 210*> \param[in] NV 211*> \verbatim 212*> NV is INTEGER 213*> The number of rows of work array WV available for 214*> workspace. NV >= NW. 215*> \endverbatim 216*> 217*> \param[out] WV 218*> \verbatim 219*> WV is COMPLEX*16 array, dimension (LDWV,NW) 220*> \endverbatim 221*> 222*> \param[in] LDWV 223*> \verbatim 224*> LDWV is INTEGER 225*> The leading dimension of W just as declared in the 226*> calling subroutine. NW <= LDV 227*> \endverbatim 228*> 229*> \param[out] WORK 230*> \verbatim 231*> WORK is COMPLEX*16 array, dimension (LWORK) 232*> On exit, WORK(1) is set to an estimate of the optimal value 233*> of LWORK for the given values of N, NW, KTOP and KBOT. 234*> \endverbatim 235*> 236*> \param[in] LWORK 237*> \verbatim 238*> LWORK is INTEGER 239*> The dimension of the work array WORK. LWORK = 2*NW 240*> suffices, but greater efficiency may result from larger 241*> values of LWORK. 242*> 243*> If LWORK = -1, then a workspace query is assumed; ZLAQR2 244*> only estimates the optimal workspace size for the given 245*> values of N, NW, KTOP and KBOT. The estimate is returned 246*> in WORK(1). No error message related to LWORK is issued 247*> by XERBLA. Neither H nor Z are accessed. 248*> \endverbatim 249* 250* Authors: 251* ======== 252* 253*> \author Univ. of Tennessee 254*> \author Univ. of California Berkeley 255*> \author Univ. of Colorado Denver 256*> \author NAG Ltd. 257* 258*> \ingroup complex16OTHERauxiliary 259* 260*> \par Contributors: 261* ================== 262*> 263*> Karen Braman and Ralph Byers, Department of Mathematics, 264*> University of Kansas, USA 265*> 266* ===================================================================== 267 SUBROUTINE ZLAQR2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, 268 $ IHIZ, Z, LDZ, NS, ND, SH, V, LDV, NH, T, LDT, 269 $ NV, WV, LDWV, WORK, LWORK ) 270* 271* -- LAPACK auxiliary routine -- 272* -- LAPACK is a software package provided by Univ. of Tennessee, -- 273* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 274* 275* .. Scalar Arguments .. 276 INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV, 277 $ LDZ, LWORK, N, ND, NH, NS, NV, NW 278 LOGICAL WANTT, WANTZ 279* .. 280* .. Array Arguments .. 281 COMPLEX*16 H( LDH, * ), SH( * ), T( LDT, * ), V( LDV, * ), 282 $ WORK( * ), WV( LDWV, * ), Z( LDZ, * ) 283* .. 284* 285* ================================================================ 286* 287* .. Parameters .. 288 COMPLEX*16 ZERO, ONE 289 PARAMETER ( ZERO = ( 0.0d0, 0.0d0 ), 290 $ ONE = ( 1.0d0, 0.0d0 ) ) 291 DOUBLE PRECISION RZERO, RONE 292 PARAMETER ( RZERO = 0.0d0, RONE = 1.0d0 ) 293* .. 294* .. Local Scalars .. 295 COMPLEX*16 BETA, CDUM, S, TAU 296 DOUBLE PRECISION FOO, SAFMAX, SAFMIN, SMLNUM, ULP 297 INTEGER I, IFST, ILST, INFO, INFQR, J, JW, KCOL, KLN, 298 $ KNT, KROW, KWTOP, LTOP, LWK1, LWK2, LWKOPT 299* .. 300* .. External Functions .. 301 DOUBLE PRECISION DLAMCH 302 EXTERNAL DLAMCH 303* .. 304* .. External Subroutines .. 305 EXTERNAL DLABAD, ZCOPY, ZGEHRD, ZGEMM, ZLACPY, ZLAHQR, 306 $ ZLARF, ZLARFG, ZLASET, ZTREXC, ZUNMHR 307* .. 308* .. Intrinsic Functions .. 309 INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, MAX, MIN 310* .. 311* .. Statement Functions .. 312 DOUBLE PRECISION CABS1 313* .. 314* .. Statement Function definitions .. 315 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) ) 316* .. 317* .. Executable Statements .. 318* 319* ==== Estimate optimal workspace. ==== 320* 321 JW = MIN( NW, KBOT-KTOP+1 ) 322 IF( JW.LE.2 ) THEN 323 LWKOPT = 1 324 ELSE 325* 326* ==== Workspace query call to ZGEHRD ==== 327* 328 CALL ZGEHRD( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO ) 329 LWK1 = INT( WORK( 1 ) ) 330* 331* ==== Workspace query call to ZUNMHR ==== 332* 333 CALL ZUNMHR( 'R', 'N', JW, JW, 1, JW-1, T, LDT, WORK, V, LDV, 334 $ WORK, -1, INFO ) 335 LWK2 = INT( WORK( 1 ) ) 336* 337* ==== Optimal workspace ==== 338* 339 LWKOPT = JW + MAX( LWK1, LWK2 ) 340 END IF 341* 342* ==== Quick return in case of workspace query. ==== 343* 344 IF( LWORK.EQ.-1 ) THEN 345 WORK( 1 ) = DCMPLX( LWKOPT, 0 ) 346 RETURN 347 END IF 348* 349* ==== Nothing to do ... 350* ... for an empty active block ... ==== 351 NS = 0 352 ND = 0 353 WORK( 1 ) = ONE 354 IF( KTOP.GT.KBOT ) 355 $ RETURN 356* ... nor for an empty deflation window. ==== 357 IF( NW.LT.1 ) 358 $ RETURN 359* 360* ==== Machine constants ==== 361* 362 SAFMIN = DLAMCH( 'SAFE MINIMUM' ) 363 SAFMAX = RONE / SAFMIN 364 CALL DLABAD( SAFMIN, SAFMAX ) 365 ULP = DLAMCH( 'PRECISION' ) 366 SMLNUM = SAFMIN*( DBLE( N ) / ULP ) 367* 368* ==== Setup deflation window ==== 369* 370 JW = MIN( NW, KBOT-KTOP+1 ) 371 KWTOP = KBOT - JW + 1 372 IF( KWTOP.EQ.KTOP ) THEN 373 S = ZERO 374 ELSE 375 S = H( KWTOP, KWTOP-1 ) 376 END IF 377* 378 IF( KBOT.EQ.KWTOP ) THEN 379* 380* ==== 1-by-1 deflation window: not much to do ==== 381* 382 SH( KWTOP ) = H( KWTOP, KWTOP ) 383 NS = 1 384 ND = 0 385 IF( CABS1( S ).LE.MAX( SMLNUM, ULP*CABS1( H( KWTOP, 386 $ KWTOP ) ) ) ) THEN 387 NS = 0 388 ND = 1 389 IF( KWTOP.GT.KTOP ) 390 $ H( KWTOP, KWTOP-1 ) = ZERO 391 END IF 392 WORK( 1 ) = ONE 393 RETURN 394 END IF 395* 396* ==== Convert to spike-triangular form. (In case of a 397* . rare QR failure, this routine continues to do 398* . aggressive early deflation using that part of 399* . the deflation window that converged using INFQR 400* . here and there to keep track.) ==== 401* 402 CALL ZLACPY( 'U', JW, JW, H( KWTOP, KWTOP ), LDH, T, LDT ) 403 CALL ZCOPY( JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ), LDT+1 ) 404* 405 CALL ZLASET( 'A', JW, JW, ZERO, ONE, V, LDV ) 406 CALL ZLAHQR( .true., .true., JW, 1, JW, T, LDT, SH( KWTOP ), 1, 407 $ JW, V, LDV, INFQR ) 408* 409* ==== Deflation detection loop ==== 410* 411 NS = JW 412 ILST = INFQR + 1 413 DO 10 KNT = INFQR + 1, JW 414* 415* ==== Small spike tip deflation test ==== 416* 417 FOO = CABS1( T( NS, NS ) ) 418 IF( FOO.EQ.RZERO ) 419 $ FOO = CABS1( S ) 420 IF( CABS1( S )*CABS1( V( 1, NS ) ).LE.MAX( SMLNUM, ULP*FOO ) ) 421 $ THEN 422* 423* ==== One more converged eigenvalue ==== 424* 425 NS = NS - 1 426 ELSE 427* 428* ==== One undeflatable eigenvalue. Move it up out of the 429* . way. (ZTREXC can not fail in this case.) ==== 430* 431 IFST = NS 432 CALL ZTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, INFO ) 433 ILST = ILST + 1 434 END IF 435 10 CONTINUE 436* 437* ==== Return to Hessenberg form ==== 438* 439 IF( NS.EQ.0 ) 440 $ S = ZERO 441* 442 IF( NS.LT.JW ) THEN 443* 444* ==== sorting the diagonal of T improves accuracy for 445* . graded matrices. ==== 446* 447 DO 30 I = INFQR + 1, NS 448 IFST = I 449 DO 20 J = I + 1, NS 450 IF( CABS1( T( J, J ) ).GT.CABS1( T( IFST, IFST ) ) ) 451 $ IFST = J 452 20 CONTINUE 453 ILST = I 454 IF( IFST.NE.ILST ) 455 $ CALL ZTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, INFO ) 456 30 CONTINUE 457 END IF 458* 459* ==== Restore shift/eigenvalue array from T ==== 460* 461 DO 40 I = INFQR + 1, JW 462 SH( KWTOP+I-1 ) = T( I, I ) 463 40 CONTINUE 464* 465* 466 IF( NS.LT.JW .OR. S.EQ.ZERO ) THEN 467 IF( NS.GT.1 .AND. S.NE.ZERO ) THEN 468* 469* ==== Reflect spike back into lower triangle ==== 470* 471 CALL ZCOPY( NS, V, LDV, WORK, 1 ) 472 DO 50 I = 1, NS 473 WORK( I ) = DCONJG( WORK( I ) ) 474 50 CONTINUE 475 BETA = WORK( 1 ) 476 CALL ZLARFG( NS, BETA, WORK( 2 ), 1, TAU ) 477 WORK( 1 ) = ONE 478* 479 CALL ZLASET( 'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT ) 480* 481 CALL ZLARF( 'L', NS, JW, WORK, 1, DCONJG( TAU ), T, LDT, 482 $ WORK( JW+1 ) ) 483 CALL ZLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT, 484 $ WORK( JW+1 ) ) 485 CALL ZLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV, 486 $ WORK( JW+1 ) ) 487* 488 CALL ZGEHRD( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ), 489 $ LWORK-JW, INFO ) 490 END IF 491* 492* ==== Copy updated reduced window into place ==== 493* 494 IF( KWTOP.GT.1 ) 495 $ H( KWTOP, KWTOP-1 ) = S*DCONJG( V( 1, 1 ) ) 496 CALL ZLACPY( 'U', JW, JW, T, LDT, H( KWTOP, KWTOP ), LDH ) 497 CALL ZCOPY( JW-1, T( 2, 1 ), LDT+1, H( KWTOP+1, KWTOP ), 498 $ LDH+1 ) 499* 500* ==== Accumulate orthogonal matrix in order update 501* . H and Z, if requested. ==== 502* 503 IF( NS.GT.1 .AND. S.NE.ZERO ) 504 $ CALL ZUNMHR( 'R', 'N', JW, NS, 1, NS, T, LDT, WORK, V, LDV, 505 $ WORK( JW+1 ), LWORK-JW, INFO ) 506* 507* ==== Update vertical slab in H ==== 508* 509 IF( WANTT ) THEN 510 LTOP = 1 511 ELSE 512 LTOP = KTOP 513 END IF 514 DO 60 KROW = LTOP, KWTOP - 1, NV 515 KLN = MIN( NV, KWTOP-KROW ) 516 CALL ZGEMM( 'N', 'N', KLN, JW, JW, ONE, H( KROW, KWTOP ), 517 $ LDH, V, LDV, ZERO, WV, LDWV ) 518 CALL ZLACPY( 'A', KLN, JW, WV, LDWV, H( KROW, KWTOP ), LDH ) 519 60 CONTINUE 520* 521* ==== Update horizontal slab in H ==== 522* 523 IF( WANTT ) THEN 524 DO 70 KCOL = KBOT + 1, N, NH 525 KLN = MIN( NH, N-KCOL+1 ) 526 CALL ZGEMM( 'C', 'N', JW, KLN, JW, ONE, V, LDV, 527 $ H( KWTOP, KCOL ), LDH, ZERO, T, LDT ) 528 CALL ZLACPY( 'A', JW, KLN, T, LDT, H( KWTOP, KCOL ), 529 $ LDH ) 530 70 CONTINUE 531 END IF 532* 533* ==== Update vertical slab in Z ==== 534* 535 IF( WANTZ ) THEN 536 DO 80 KROW = ILOZ, IHIZ, NV 537 KLN = MIN( NV, IHIZ-KROW+1 ) 538 CALL ZGEMM( 'N', 'N', KLN, JW, JW, ONE, Z( KROW, KWTOP ), 539 $ LDZ, V, LDV, ZERO, WV, LDWV ) 540 CALL ZLACPY( 'A', KLN, JW, WV, LDWV, Z( KROW, KWTOP ), 541 $ LDZ ) 542 80 CONTINUE 543 END IF 544 END IF 545* 546* ==== Return the number of deflations ... ==== 547* 548 ND = JW - NS 549* 550* ==== ... and the number of shifts. (Subtracting 551* . INFQR from the spike length takes care 552* . of the case of a rare QR failure while 553* . calculating eigenvalues of the deflation 554* . window.) ==== 555* 556 NS = NS - INFQR 557* 558* ==== Return optimal workspace. ==== 559* 560 WORK( 1 ) = DCMPLX( LWKOPT, 0 ) 561* 562* ==== End of ZLAQR2 ==== 563* 564 END 565