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