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