1*> \brief \b DLAQR3 performs the orthogonal 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 DLAQR3 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqr3.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqr3.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqr3.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLAQR3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, 22* IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T, 23* LDT, 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* DOUBLE PRECISION H( LDH, * ), SI( * ), SR( * ), T( LDT, * ), 32* $ V( LDV, * ), WORK( * ), WV( LDWV, * ), 33* $ Z( LDZ, * ) 34* .. 35* 36* 37*> \par Purpose: 38* ============= 39*> 40*> \verbatim 41*> 42*> Aggressive early deflation: 43*> 44*> DLAQR3 accepts as input an upper Hessenberg matrix 45*> H and performs an orthogonal similarity transformation 46*> designed to detect and deflate fully converged eigenvalues from 47*> a trailing principal submatrix. On output H has been over- 48*> written by a new Hessenberg matrix that is a perturbation of 49*> an orthogonal similarity transformation of H. It is to be 50*> hoped that the final version of H has many zero subdiagonal 51*> entries. 52*> \endverbatim 53* 54* Arguments: 55* ========== 56* 57*> \param[in] WANTT 58*> \verbatim 59*> WANTT is LOGICAL 60*> If .TRUE., then the Hessenberg matrix H is fully updated 61*> so that the quasi-triangular Schur factor may be 62*> computed (in cooperation with the calling subroutine). 63*> If .FALSE., then only enough of H is updated to preserve 64*> the eigenvalues. 65*> \endverbatim 66*> 67*> \param[in] WANTZ 68*> \verbatim 69*> WANTZ is LOGICAL 70*> If .TRUE., then the orthogonal matrix Z is updated so 71*> so that the orthogonal Schur factor may be computed 72*> (in cooperation with the calling subroutine). 73*> If .FALSE., then Z is not referenced. 74*> \endverbatim 75*> 76*> \param[in] N 77*> \verbatim 78*> N is INTEGER 79*> The order of the matrix H and (if WANTZ is .TRUE.) the 80*> order of the orthogonal matrix Z. 81*> \endverbatim 82*> 83*> \param[in] KTOP 84*> \verbatim 85*> KTOP is INTEGER 86*> It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0. 87*> KBOT and KTOP together determine an isolated block 88*> along the diagonal of the Hessenberg matrix. 89*> \endverbatim 90*> 91*> \param[in] KBOT 92*> \verbatim 93*> KBOT is INTEGER 94*> It is assumed without a check that either 95*> KBOT = N or H(KBOT+1,KBOT)=0. KBOT and KTOP together 96*> determine an isolated block along the diagonal of the 97*> Hessenberg matrix. 98*> \endverbatim 99*> 100*> \param[in] NW 101*> \verbatim 102*> NW is INTEGER 103*> Deflation window size. 1 .LE. NW .LE. (KBOT-KTOP+1). 104*> \endverbatim 105*> 106*> \param[in,out] H 107*> \verbatim 108*> H is DOUBLE PRECISION array, dimension (LDH,N) 109*> On input the initial N-by-N section of H stores the 110*> Hessenberg matrix undergoing aggressive early deflation. 111*> On output H has been transformed by an orthogonal 112*> similarity transformation, perturbed, and the returned 113*> to Hessenberg form that (it is to be hoped) has some 114*> zero subdiagonal entries. 115*> \endverbatim 116*> 117*> \param[in] LDH 118*> \verbatim 119*> LDH is integer 120*> Leading dimension of H just as declared in the calling 121*> subroutine. N .LE. LDH 122*> \endverbatim 123*> 124*> \param[in] ILOZ 125*> \verbatim 126*> ILOZ is INTEGER 127*> \endverbatim 128*> 129*> \param[in] IHIZ 130*> \verbatim 131*> IHIZ is INTEGER 132*> Specify the rows of Z to which transformations must be 133*> applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N. 134*> \endverbatim 135*> 136*> \param[in,out] Z 137*> \verbatim 138*> Z is DOUBLE PRECISION array, dimension (LDZ,N) 139*> IF WANTZ is .TRUE., then on output, the orthogonal 140*> similarity transformation mentioned above has been 141*> accumulated into Z(ILOZ:IHIZ,ILO:IHI) from the right. 142*> If WANTZ is .FALSE., then Z is unreferenced. 143*> \endverbatim 144*> 145*> \param[in] LDZ 146*> \verbatim 147*> LDZ is integer 148*> The leading dimension of Z just as declared in the 149*> calling subroutine. 1 .LE. LDZ. 150*> \endverbatim 151*> 152*> \param[out] NS 153*> \verbatim 154*> NS is integer 155*> The number of unconverged (ie approximate) eigenvalues 156*> returned in SR and SI that may be used as shifts by the 157*> calling subroutine. 158*> \endverbatim 159*> 160*> \param[out] ND 161*> \verbatim 162*> ND is integer 163*> The number of converged eigenvalues uncovered by this 164*> subroutine. 165*> \endverbatim 166*> 167*> \param[out] SR 168*> \verbatim 169*> SR is DOUBLE PRECISION array, dimension (KBOT) 170*> \endverbatim 171*> 172*> \param[out] SI 173*> \verbatim 174*> SI is DOUBLE PRECISION array, dimension (KBOT) 175*> On output, the real and imaginary parts of approximate 176*> eigenvalues that may be used for shifts are stored in 177*> SR(KBOT-ND-NS+1) through SR(KBOT-ND) and 178*> SI(KBOT-ND-NS+1) through SI(KBOT-ND), respectively. 179*> The real and imaginary parts of converged eigenvalues 180*> are stored in SR(KBOT-ND+1) through SR(KBOT) and 181*> SI(KBOT-ND+1) through SI(KBOT), respectively. 182*> \endverbatim 183*> 184*> \param[out] V 185*> \verbatim 186*> V is DOUBLE PRECISION array, dimension (LDV,NW) 187*> An NW-by-NW work array. 188*> \endverbatim 189*> 190*> \param[in] LDV 191*> \verbatim 192*> LDV is integer scalar 193*> The leading dimension of V just as declared in the 194*> calling subroutine. NW .LE. LDV 195*> \endverbatim 196*> 197*> \param[in] NH 198*> \verbatim 199*> NH is integer scalar 200*> The number of columns of T. NH.GE.NW. 201*> \endverbatim 202*> 203*> \param[out] T 204*> \verbatim 205*> T is DOUBLE PRECISION array, dimension (LDT,NW) 206*> \endverbatim 207*> 208*> \param[in] LDT 209*> \verbatim 210*> LDT is integer 211*> The leading dimension of T just as declared in the 212*> calling subroutine. NW .LE. LDT 213*> \endverbatim 214*> 215*> \param[in] NV 216*> \verbatim 217*> NV is integer 218*> The number of rows of work array WV available for 219*> workspace. NV.GE.NW. 220*> \endverbatim 221*> 222*> \param[out] WV 223*> \verbatim 224*> WV is DOUBLE PRECISION array, dimension (LDWV,NW) 225*> \endverbatim 226*> 227*> \param[in] LDWV 228*> \verbatim 229*> LDWV is integer 230*> The leading dimension of W just as declared in the 231*> calling subroutine. NW .LE. LDV 232*> \endverbatim 233*> 234*> \param[out] WORK 235*> \verbatim 236*> WORK is DOUBLE PRECISION array, dimension (LWORK) 237*> On exit, WORK(1) is set to an estimate of the optimal value 238*> of LWORK for the given values of N, NW, KTOP and KBOT. 239*> \endverbatim 240*> 241*> \param[in] LWORK 242*> \verbatim 243*> LWORK is integer 244*> The dimension of the work array WORK. LWORK = 2*NW 245*> suffices, but greater efficiency may result from larger 246*> values of LWORK. 247*> 248*> If LWORK = -1, then a workspace query is assumed; DLAQR3 249*> only estimates the optimal workspace size for the given 250*> values of N, NW, KTOP and KBOT. The estimate is returned 251*> in WORK(1). No error message related to LWORK is issued 252*> by XERBLA. Neither H nor Z are accessed. 253*> \endverbatim 254* 255* Authors: 256* ======== 257* 258*> \author Univ. of Tennessee 259*> \author Univ. of California Berkeley 260*> \author Univ. of Colorado Denver 261*> \author NAG Ltd. 262* 263*> \date September 2012 264* 265*> \ingroup doubleOTHERauxiliary 266* 267*> \par Contributors: 268* ================== 269*> 270*> Karen Braman and Ralph Byers, Department of Mathematics, 271*> University of Kansas, USA 272*> 273* ===================================================================== 274 SUBROUTINE DLAQR3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, 275 $ IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T, 276 $ LDT, NV, WV, LDWV, WORK, LWORK ) 277* 278* -- LAPACK auxiliary routine (version 3.4.2) -- 279* -- LAPACK is a software package provided by Univ. of Tennessee, -- 280* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 281* September 2012 282* 283* .. Scalar Arguments .. 284 INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV, 285 $ LDZ, LWORK, N, ND, NH, NS, NV, NW 286 LOGICAL WANTT, WANTZ 287* .. 288* .. Array Arguments .. 289 DOUBLE PRECISION H( LDH, * ), SI( * ), SR( * ), T( LDT, * ), 290 $ V( LDV, * ), WORK( * ), WV( LDWV, * ), 291 $ Z( LDZ, * ) 292* .. 293* 294* ================================================================ 295* .. Parameters .. 296 DOUBLE PRECISION ZERO, ONE 297 PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0 ) 298* .. 299* .. Local Scalars .. 300 DOUBLE PRECISION AA, BB, BETA, CC, CS, DD, EVI, EVK, FOO, S, 301 $ SAFMAX, SAFMIN, SMLNUM, SN, TAU, ULP 302 INTEGER I, IFST, ILST, INFO, INFQR, J, JW, K, KCOL, 303 $ KEND, KLN, KROW, KWTOP, LTOP, LWK1, LWK2, LWK3, 304 $ LWKOPT, NMIN 305 LOGICAL BULGE, SORTED 306* .. 307* .. External Functions .. 308 DOUBLE PRECISION DLAMCH 309 INTEGER ILAENV 310 EXTERNAL DLAMCH, ILAENV 311* .. 312* .. External Subroutines .. 313 EXTERNAL DCOPY, DGEHRD, DGEMM, DLABAD, DLACPY, DLAHQR, 314 $ DLANV2, DLAQR4, DLARF, DLARFG, DLASET, DORMHR, 315 $ DTREXC 316* .. 317* .. Intrinsic Functions .. 318 INTRINSIC ABS, DBLE, INT, MAX, MIN, SQRT 319* .. 320* .. Executable Statements .. 321* 322* ==== Estimate optimal workspace. ==== 323* 324 JW = MIN( NW, KBOT-KTOP+1 ) 325 IF( JW.LE.2 ) THEN 326 LWKOPT = 1 327 ELSE 328* 329* ==== Workspace query call to DGEHRD ==== 330* 331 CALL DGEHRD( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO ) 332 LWK1 = INT( WORK( 1 ) ) 333* 334* ==== Workspace query call to DORMHR ==== 335* 336 CALL DORMHR( 'R', 'N', JW, JW, 1, JW-1, T, LDT, WORK, V, LDV, 337 $ WORK, -1, INFO ) 338 LWK2 = INT( WORK( 1 ) ) 339* 340* ==== Workspace query call to DLAQR4 ==== 341* 342 CALL DLAQR4( .true., .true., JW, 1, JW, T, LDT, SR, SI, 1, JW, 343 $ V, LDV, WORK, -1, INFQR ) 344 LWK3 = INT( WORK( 1 ) ) 345* 346* ==== Optimal workspace ==== 347* 348 LWKOPT = MAX( JW+MAX( LWK1, LWK2 ), LWK3 ) 349 END IF 350* 351* ==== Quick return in case of workspace query. ==== 352* 353 IF( LWORK.EQ.-1 ) THEN 354 WORK( 1 ) = DBLE( LWKOPT ) 355 RETURN 356 END IF 357* 358* ==== Nothing to do ... 359* ... for an empty active block ... ==== 360 NS = 0 361 ND = 0 362 WORK( 1 ) = ONE 363 IF( KTOP.GT.KBOT ) 364 $ RETURN 365* ... nor for an empty deflation window. ==== 366 IF( NW.LT.1 ) 367 $ RETURN 368* 369* ==== Machine constants ==== 370* 371 SAFMIN = DLAMCH( 'SAFE MINIMUM' ) 372 SAFMAX = ONE / SAFMIN 373 CALL DLABAD( SAFMIN, SAFMAX ) 374 ULP = DLAMCH( 'PRECISION' ) 375 SMLNUM = SAFMIN*( DBLE( N ) / ULP ) 376* 377* ==== Setup deflation window ==== 378* 379 JW = MIN( NW, KBOT-KTOP+1 ) 380 KWTOP = KBOT - JW + 1 381 IF( KWTOP.EQ.KTOP ) THEN 382 S = ZERO 383 ELSE 384 S = H( KWTOP, KWTOP-1 ) 385 END IF 386* 387 IF( KBOT.EQ.KWTOP ) THEN 388* 389* ==== 1-by-1 deflation window: not much to do ==== 390* 391 SR( KWTOP ) = H( KWTOP, KWTOP ) 392 SI( KWTOP ) = ZERO 393 NS = 1 394 ND = 0 395 IF( ABS( S ).LE.MAX( SMLNUM, ULP*ABS( H( KWTOP, KWTOP ) ) ) ) 396 $ THEN 397 NS = 0 398 ND = 1 399 IF( KWTOP.GT.KTOP ) 400 $ H( KWTOP, KWTOP-1 ) = ZERO 401 END IF 402 WORK( 1 ) = ONE 403 RETURN 404 END IF 405* 406* ==== Convert to spike-triangular form. (In case of a 407* . rare QR failure, this routine continues to do 408* . aggressive early deflation using that part of 409* . the deflation window that converged using INFQR 410* . here and there to keep track.) ==== 411* 412 CALL DLACPY( 'U', JW, JW, H( KWTOP, KWTOP ), LDH, T, LDT ) 413 CALL DCOPY( JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ), LDT+1 ) 414* 415 CALL DLASET( 'A', JW, JW, ZERO, ONE, V, LDV ) 416 NMIN = ILAENV( 12, 'DLAQR3', 'SV', JW, 1, JW, LWORK ) 417 IF( JW.GT.NMIN ) THEN 418 CALL DLAQR4( .true., .true., JW, 1, JW, T, LDT, SR( KWTOP ), 419 $ SI( KWTOP ), 1, JW, V, LDV, WORK, LWORK, INFQR ) 420 ELSE 421 CALL DLAHQR( .true., .true., JW, 1, JW, T, LDT, SR( KWTOP ), 422 $ SI( KWTOP ), 1, JW, V, LDV, INFQR ) 423 END IF 424* 425* ==== DTREXC needs a clean margin near the diagonal ==== 426* 427 DO 10 J = 1, JW - 3 428 T( J+2, J ) = ZERO 429 T( J+3, J ) = ZERO 430 10 CONTINUE 431 IF( JW.GT.2 ) 432 $ T( JW, JW-2 ) = ZERO 433* 434* ==== Deflation detection loop ==== 435* 436 NS = JW 437 ILST = INFQR + 1 438 20 CONTINUE 439 IF( ILST.LE.NS ) THEN 440 IF( NS.EQ.1 ) THEN 441 BULGE = .FALSE. 442 ELSE 443 BULGE = T( NS, NS-1 ).NE.ZERO 444 END IF 445* 446* ==== Small spike tip test for deflation ==== 447* 448 IF( .NOT. BULGE ) THEN 449* 450* ==== Real eigenvalue ==== 451* 452 FOO = ABS( T( NS, NS ) ) 453 IF( FOO.EQ.ZERO ) 454 $ FOO = ABS( S ) 455 IF( ABS( S*V( 1, NS ) ).LE.MAX( SMLNUM, ULP*FOO ) ) THEN 456* 457* ==== Deflatable ==== 458* 459 NS = NS - 1 460 ELSE 461* 462* ==== Undeflatable. Move it up out of the way. 463* . (DTREXC can not fail in this case.) ==== 464* 465 IFST = NS 466 CALL DTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK, 467 $ INFO ) 468 ILST = ILST + 1 469 END IF 470 ELSE 471* 472* ==== Complex conjugate pair ==== 473* 474 FOO = ABS( T( NS, NS ) ) + SQRT( ABS( T( NS, NS-1 ) ) )* 475 $ SQRT( ABS( T( NS-1, NS ) ) ) 476 IF( FOO.EQ.ZERO ) 477 $ FOO = ABS( S ) 478 IF( MAX( ABS( S*V( 1, NS ) ), ABS( S*V( 1, NS-1 ) ) ).LE. 479 $ MAX( SMLNUM, ULP*FOO ) ) THEN 480* 481* ==== Deflatable ==== 482* 483 NS = NS - 2 484 ELSE 485* 486* ==== Undeflatable. Move them up out of the way. 487* . Fortunately, DTREXC does the right thing with 488* . ILST in case of a rare exchange failure. ==== 489* 490 IFST = NS 491 CALL DTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK, 492 $ INFO ) 493 ILST = ILST + 2 494 END IF 495 END IF 496* 497* ==== End deflation detection loop ==== 498* 499 GO TO 20 500 END IF 501* 502* ==== Return to Hessenberg form ==== 503* 504 IF( NS.EQ.0 ) 505 $ S = ZERO 506* 507 IF( NS.LT.JW ) THEN 508* 509* ==== sorting diagonal blocks of T improves accuracy for 510* . graded matrices. Bubble sort deals well with 511* . exchange failures. ==== 512* 513 SORTED = .false. 514 I = NS + 1 515 30 CONTINUE 516 IF( SORTED ) 517 $ GO TO 50 518 SORTED = .true. 519* 520 KEND = I - 1 521 I = INFQR + 1 522 IF( I.EQ.NS ) THEN 523 K = I + 1 524 ELSE IF( T( I+1, I ).EQ.ZERO ) THEN 525 K = I + 1 526 ELSE 527 K = I + 2 528 END IF 529 40 CONTINUE 530 IF( K.LE.KEND ) THEN 531 IF( K.EQ.I+1 ) THEN 532 EVI = ABS( T( I, I ) ) 533 ELSE 534 EVI = ABS( T( I, I ) ) + SQRT( ABS( T( I+1, I ) ) )* 535 $ SQRT( ABS( T( I, I+1 ) ) ) 536 END IF 537* 538 IF( K.EQ.KEND ) THEN 539 EVK = ABS( T( K, K ) ) 540 ELSE IF( T( K+1, K ).EQ.ZERO ) THEN 541 EVK = ABS( T( K, K ) ) 542 ELSE 543 EVK = ABS( T( K, K ) ) + SQRT( ABS( T( K+1, K ) ) )* 544 $ SQRT( ABS( T( K, K+1 ) ) ) 545 END IF 546* 547 IF( EVI.GE.EVK ) THEN 548 I = K 549 ELSE 550 SORTED = .false. 551 IFST = I 552 ILST = K 553 CALL DTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK, 554 $ INFO ) 555 IF( INFO.EQ.0 ) THEN 556 I = ILST 557 ELSE 558 I = K 559 END IF 560 END IF 561 IF( I.EQ.KEND ) THEN 562 K = I + 1 563 ELSE IF( T( I+1, I ).EQ.ZERO ) THEN 564 K = I + 1 565 ELSE 566 K = I + 2 567 END IF 568 GO TO 40 569 END IF 570 GO TO 30 571 50 CONTINUE 572 END IF 573* 574* ==== Restore shift/eigenvalue array from T ==== 575* 576 I = JW 577 60 CONTINUE 578 IF( I.GE.INFQR+1 ) THEN 579 IF( I.EQ.INFQR+1 ) THEN 580 SR( KWTOP+I-1 ) = T( I, I ) 581 SI( KWTOP+I-1 ) = ZERO 582 I = I - 1 583 ELSE IF( T( I, I-1 ).EQ.ZERO ) THEN 584 SR( KWTOP+I-1 ) = T( I, I ) 585 SI( KWTOP+I-1 ) = ZERO 586 I = I - 1 587 ELSE 588 AA = T( I-1, I-1 ) 589 CC = T( I, I-1 ) 590 BB = T( I-1, I ) 591 DD = T( I, I ) 592 CALL DLANV2( AA, BB, CC, DD, SR( KWTOP+I-2 ), 593 $ SI( KWTOP+I-2 ), SR( KWTOP+I-1 ), 594 $ SI( KWTOP+I-1 ), CS, SN ) 595 I = I - 2 596 END IF 597 GO TO 60 598 END IF 599* 600 IF( NS.LT.JW .OR. S.EQ.ZERO ) THEN 601 IF( NS.GT.1 .AND. S.NE.ZERO ) THEN 602* 603* ==== Reflect spike back into lower triangle ==== 604* 605 CALL DCOPY( NS, V, LDV, WORK, 1 ) 606 BETA = WORK( 1 ) 607 CALL DLARFG( NS, BETA, WORK( 2 ), 1, TAU ) 608 WORK( 1 ) = ONE 609* 610 CALL DLASET( 'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT ) 611* 612 CALL DLARF( 'L', NS, JW, WORK, 1, TAU, T, LDT, 613 $ WORK( JW+1 ) ) 614 CALL DLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT, 615 $ WORK( JW+1 ) ) 616 CALL DLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV, 617 $ WORK( JW+1 ) ) 618* 619 CALL DGEHRD( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ), 620 $ LWORK-JW, INFO ) 621 END IF 622* 623* ==== Copy updated reduced window into place ==== 624* 625 IF( KWTOP.GT.1 ) 626 $ H( KWTOP, KWTOP-1 ) = S*V( 1, 1 ) 627 CALL DLACPY( 'U', JW, JW, T, LDT, H( KWTOP, KWTOP ), LDH ) 628 CALL DCOPY( JW-1, T( 2, 1 ), LDT+1, H( KWTOP+1, KWTOP ), 629 $ LDH+1 ) 630* 631* ==== Accumulate orthogonal matrix in order update 632* . H and Z, if requested. ==== 633* 634 IF( NS.GT.1 .AND. S.NE.ZERO ) 635 $ CALL DORMHR( 'R', 'N', JW, NS, 1, NS, T, LDT, WORK, V, LDV, 636 $ WORK( JW+1 ), LWORK-JW, INFO ) 637* 638* ==== Update vertical slab in H ==== 639* 640 IF( WANTT ) THEN 641 LTOP = 1 642 ELSE 643 LTOP = KTOP 644 END IF 645 DO 70 KROW = LTOP, KWTOP - 1, NV 646 KLN = MIN( NV, KWTOP-KROW ) 647 CALL DGEMM( 'N', 'N', KLN, JW, JW, ONE, H( KROW, KWTOP ), 648 $ LDH, V, LDV, ZERO, WV, LDWV ) 649 CALL DLACPY( 'A', KLN, JW, WV, LDWV, H( KROW, KWTOP ), LDH ) 650 70 CONTINUE 651* 652* ==== Update horizontal slab in H ==== 653* 654 IF( WANTT ) THEN 655 DO 80 KCOL = KBOT + 1, N, NH 656 KLN = MIN( NH, N-KCOL+1 ) 657 CALL DGEMM( 'C', 'N', JW, KLN, JW, ONE, V, LDV, 658 $ H( KWTOP, KCOL ), LDH, ZERO, T, LDT ) 659 CALL DLACPY( 'A', JW, KLN, T, LDT, H( KWTOP, KCOL ), 660 $ LDH ) 661 80 CONTINUE 662 END IF 663* 664* ==== Update vertical slab in Z ==== 665* 666 IF( WANTZ ) THEN 667 DO 90 KROW = ILOZ, IHIZ, NV 668 KLN = MIN( NV, IHIZ-KROW+1 ) 669 CALL DGEMM( 'N', 'N', KLN, JW, JW, ONE, Z( KROW, KWTOP ), 670 $ LDZ, V, LDV, ZERO, WV, LDWV ) 671 CALL DLACPY( 'A', KLN, JW, WV, LDWV, Z( KROW, KWTOP ), 672 $ LDZ ) 673 90 CONTINUE 674 END IF 675 END IF 676* 677* ==== Return the number of deflations ... ==== 678* 679 ND = JW - NS 680* 681* ==== ... and the number of shifts. (Subtracting 682* . INFQR from the spike length takes care 683* . of the case of a rare QR failure while 684* . calculating eigenvalues of the deflation 685* . window.) ==== 686* 687 NS = NS - INFQR 688* 689* ==== Return optimal workspace. ==== 690* 691 WORK( 1 ) = DBLE( LWKOPT ) 692* 693* ==== End of DLAQR3 ==== 694* 695 END 696c $Id$ 697