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 <= NW <= (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 <= 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 <= ILOZ <= IHIZ <= 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,ILOZ:IHIZ) 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 <= 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 193*> The leading dimension of V just as declared in the 194*> calling subroutine. NW <= LDV 195*> \endverbatim 196*> 197*> \param[in] NH 198*> \verbatim 199*> NH is INTEGER 200*> The number of columns of T. NH >= 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 <= 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 >= 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 <= 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*> \ingroup doubleOTHERauxiliary 264* 265*> \par Contributors: 266* ================== 267*> 268*> Karen Braman and Ralph Byers, Department of Mathematics, 269*> University of Kansas, USA 270*> 271* ===================================================================== 272 SUBROUTINE DLAQR3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, 273 $ IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T, 274 $ LDT, NV, WV, LDWV, WORK, LWORK ) 275* 276* -- LAPACK auxiliary routine -- 277* -- LAPACK is a software package provided by Univ. of Tennessee, -- 278* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 279* 280* .. Scalar Arguments .. 281 INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV, 282 $ LDZ, LWORK, N, ND, NH, NS, NV, NW 283 LOGICAL WANTT, WANTZ 284* .. 285* .. Array Arguments .. 286 DOUBLE PRECISION H( LDH, * ), SI( * ), SR( * ), T( LDT, * ), 287 $ V( LDV, * ), WORK( * ), WV( LDWV, * ), 288 $ Z( LDZ, * ) 289* .. 290* 291* ================================================================ 292* .. Parameters .. 293 DOUBLE PRECISION ZERO, ONE 294 PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0 ) 295* .. 296* .. Local Scalars .. 297 DOUBLE PRECISION AA, BB, BETA, CC, CS, DD, EVI, EVK, FOO, S, 298 $ SAFMAX, SAFMIN, SMLNUM, SN, TAU, ULP 299 INTEGER I, IFST, ILST, INFO, INFQR, J, JW, K, KCOL, 300 $ KEND, KLN, KROW, KWTOP, LTOP, LWK1, LWK2, LWK3, 301 $ LWKOPT, NMIN 302 LOGICAL BULGE, SORTED 303* .. 304* .. External Functions .. 305 DOUBLE PRECISION DLAMCH 306 INTEGER ILAENV 307 EXTERNAL DLAMCH, ILAENV 308* .. 309* .. External Subroutines .. 310 EXTERNAL DCOPY, DGEHRD, DGEMM, DLABAD, DLACPY, DLAHQR, 311 $ DLANV2, DLAQR4, DLARF, DLARFG, DLASET, DORMHR, 312 $ DTREXC 313* .. 314* .. Intrinsic Functions .. 315 INTRINSIC ABS, DBLE, INT, MAX, MIN, SQRT 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 DGEHRD ==== 327* 328 CALL DGEHRD( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO ) 329 LWK1 = INT( WORK( 1 ) ) 330* 331* ==== Workspace query call to DORMHR ==== 332* 333 CALL DORMHR( 'R', 'N', JW, JW, 1, JW-1, T, LDT, WORK, V, LDV, 334 $ WORK, -1, INFO ) 335 LWK2 = INT( WORK( 1 ) ) 336* 337* ==== Workspace query call to DLAQR4 ==== 338* 339 CALL DLAQR4( .true., .true., JW, 1, JW, T, LDT, SR, SI, 1, JW, 340 $ V, LDV, WORK, -1, INFQR ) 341 LWK3 = INT( WORK( 1 ) ) 342* 343* ==== Optimal workspace ==== 344* 345 LWKOPT = MAX( JW+MAX( LWK1, LWK2 ), LWK3 ) 346 END IF 347* 348* ==== Quick return in case of workspace query. ==== 349* 350 IF( LWORK.EQ.-1 ) THEN 351 WORK( 1 ) = DBLE( LWKOPT ) 352 RETURN 353 END IF 354* 355* ==== Nothing to do ... 356* ... for an empty active block ... ==== 357 NS = 0 358 ND = 0 359 WORK( 1 ) = ONE 360 IF( KTOP.GT.KBOT ) 361 $ RETURN 362* ... nor for an empty deflation window. ==== 363 IF( NW.LT.1 ) 364 $ RETURN 365* 366* ==== Machine constants ==== 367* 368 SAFMIN = DLAMCH( 'SAFE MINIMUM' ) 369 SAFMAX = ONE / SAFMIN 370 CALL DLABAD( SAFMIN, SAFMAX ) 371 ULP = DLAMCH( 'PRECISION' ) 372 SMLNUM = SAFMIN*( DBLE( N ) / ULP ) 373* 374* ==== Setup deflation window ==== 375* 376 JW = MIN( NW, KBOT-KTOP+1 ) 377 KWTOP = KBOT - JW + 1 378 IF( KWTOP.EQ.KTOP ) THEN 379 S = ZERO 380 ELSE 381 S = H( KWTOP, KWTOP-1 ) 382 END IF 383* 384 IF( KBOT.EQ.KWTOP ) THEN 385* 386* ==== 1-by-1 deflation window: not much to do ==== 387* 388 SR( KWTOP ) = H( KWTOP, KWTOP ) 389 SI( KWTOP ) = ZERO 390 NS = 1 391 ND = 0 392 IF( ABS( S ).LE.MAX( SMLNUM, ULP*ABS( H( KWTOP, KWTOP ) ) ) ) 393 $ THEN 394 NS = 0 395 ND = 1 396 IF( KWTOP.GT.KTOP ) 397 $ H( KWTOP, KWTOP-1 ) = ZERO 398 END IF 399 WORK( 1 ) = ONE 400 RETURN 401 END IF 402* 403* ==== Convert to spike-triangular form. (In case of a 404* . rare QR failure, this routine continues to do 405* . aggressive early deflation using that part of 406* . the deflation window that converged using INFQR 407* . here and there to keep track.) ==== 408* 409 CALL DLACPY( 'U', JW, JW, H( KWTOP, KWTOP ), LDH, T, LDT ) 410 CALL DCOPY( JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ), LDT+1 ) 411* 412 CALL DLASET( 'A', JW, JW, ZERO, ONE, V, LDV ) 413 NMIN = ILAENV( 12, 'DLAQR3', 'SV', JW, 1, JW, LWORK ) 414 IF( JW.GT.NMIN ) THEN 415 CALL DLAQR4( .true., .true., JW, 1, JW, T, LDT, SR( KWTOP ), 416 $ SI( KWTOP ), 1, JW, V, LDV, WORK, LWORK, INFQR ) 417 ELSE 418 CALL DLAHQR( .true., .true., JW, 1, JW, T, LDT, SR( KWTOP ), 419 $ SI( KWTOP ), 1, JW, V, LDV, INFQR ) 420 END IF 421* 422* ==== DTREXC needs a clean margin near the diagonal ==== 423* 424 DO 10 J = 1, JW - 3 425 T( J+2, J ) = ZERO 426 T( J+3, J ) = ZERO 427 10 CONTINUE 428 IF( JW.GT.2 ) 429 $ T( JW, JW-2 ) = ZERO 430* 431* ==== Deflation detection loop ==== 432* 433 NS = JW 434 ILST = INFQR + 1 435 20 CONTINUE 436 IF( ILST.LE.NS ) THEN 437 IF( NS.EQ.1 ) THEN 438 BULGE = .FALSE. 439 ELSE 440 BULGE = T( NS, NS-1 ).NE.ZERO 441 END IF 442* 443* ==== Small spike tip test for deflation ==== 444* 445 IF( .NOT. BULGE ) THEN 446* 447* ==== Real eigenvalue ==== 448* 449 FOO = ABS( T( NS, NS ) ) 450 IF( FOO.EQ.ZERO ) 451 $ FOO = ABS( S ) 452 IF( ABS( S*V( 1, NS ) ).LE.MAX( SMLNUM, ULP*FOO ) ) THEN 453* 454* ==== Deflatable ==== 455* 456 NS = NS - 1 457 ELSE 458* 459* ==== Undeflatable. Move it up out of the way. 460* . (DTREXC can not fail in this case.) ==== 461* 462 IFST = NS 463 CALL DTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK, 464 $ INFO ) 465 ILST = ILST + 1 466 END IF 467 ELSE 468* 469* ==== Complex conjugate pair ==== 470* 471 FOO = ABS( T( NS, NS ) ) + SQRT( ABS( T( NS, NS-1 ) ) )* 472 $ SQRT( ABS( T( NS-1, NS ) ) ) 473 IF( FOO.EQ.ZERO ) 474 $ FOO = ABS( S ) 475 IF( MAX( ABS( S*V( 1, NS ) ), ABS( S*V( 1, NS-1 ) ) ).LE. 476 $ MAX( SMLNUM, ULP*FOO ) ) THEN 477* 478* ==== Deflatable ==== 479* 480 NS = NS - 2 481 ELSE 482* 483* ==== Undeflatable. Move them up out of the way. 484* . Fortunately, DTREXC does the right thing with 485* . ILST in case of a rare exchange failure. ==== 486* 487 IFST = NS 488 CALL DTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK, 489 $ INFO ) 490 ILST = ILST + 2 491 END IF 492 END IF 493* 494* ==== End deflation detection loop ==== 495* 496 GO TO 20 497 END IF 498* 499* ==== Return to Hessenberg form ==== 500* 501 IF( NS.EQ.0 ) 502 $ S = ZERO 503* 504 IF( NS.LT.JW ) THEN 505* 506* ==== sorting diagonal blocks of T improves accuracy for 507* . graded matrices. Bubble sort deals well with 508* . exchange failures. ==== 509* 510 SORTED = .false. 511 I = NS + 1 512 30 CONTINUE 513 IF( SORTED ) 514 $ GO TO 50 515 SORTED = .true. 516* 517 KEND = I - 1 518 I = INFQR + 1 519 IF( I.EQ.NS ) THEN 520 K = I + 1 521 ELSE IF( T( I+1, I ).EQ.ZERO ) THEN 522 K = I + 1 523 ELSE 524 K = I + 2 525 END IF 526 40 CONTINUE 527 IF( K.LE.KEND ) THEN 528 IF( K.EQ.I+1 ) THEN 529 EVI = ABS( T( I, I ) ) 530 ELSE 531 EVI = ABS( T( I, I ) ) + SQRT( ABS( T( I+1, I ) ) )* 532 $ SQRT( ABS( T( I, I+1 ) ) ) 533 END IF 534* 535 IF( K.EQ.KEND ) THEN 536 EVK = ABS( T( K, K ) ) 537 ELSE IF( T( K+1, K ).EQ.ZERO ) THEN 538 EVK = ABS( T( K, K ) ) 539 ELSE 540 EVK = ABS( T( K, K ) ) + SQRT( ABS( T( K+1, K ) ) )* 541 $ SQRT( ABS( T( K, K+1 ) ) ) 542 END IF 543* 544 IF( EVI.GE.EVK ) THEN 545 I = K 546 ELSE 547 SORTED = .false. 548 IFST = I 549 ILST = K 550 CALL DTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK, 551 $ INFO ) 552 IF( INFO.EQ.0 ) THEN 553 I = ILST 554 ELSE 555 I = K 556 END IF 557 END IF 558 IF( I.EQ.KEND ) THEN 559 K = I + 1 560 ELSE IF( T( I+1, I ).EQ.ZERO ) THEN 561 K = I + 1 562 ELSE 563 K = I + 2 564 END IF 565 GO TO 40 566 END IF 567 GO TO 30 568 50 CONTINUE 569 END IF 570* 571* ==== Restore shift/eigenvalue array from T ==== 572* 573 I = JW 574 60 CONTINUE 575 IF( I.GE.INFQR+1 ) THEN 576 IF( I.EQ.INFQR+1 ) THEN 577 SR( KWTOP+I-1 ) = T( I, I ) 578 SI( KWTOP+I-1 ) = ZERO 579 I = I - 1 580 ELSE IF( T( I, I-1 ).EQ.ZERO ) THEN 581 SR( KWTOP+I-1 ) = T( I, I ) 582 SI( KWTOP+I-1 ) = ZERO 583 I = I - 1 584 ELSE 585 AA = T( I-1, I-1 ) 586 CC = T( I, I-1 ) 587 BB = T( I-1, I ) 588 DD = T( I, I ) 589 CALL DLANV2( AA, BB, CC, DD, SR( KWTOP+I-2 ), 590 $ SI( KWTOP+I-2 ), SR( KWTOP+I-1 ), 591 $ SI( KWTOP+I-1 ), CS, SN ) 592 I = I - 2 593 END IF 594 GO TO 60 595 END IF 596* 597 IF( NS.LT.JW .OR. S.EQ.ZERO ) THEN 598 IF( NS.GT.1 .AND. S.NE.ZERO ) THEN 599* 600* ==== Reflect spike back into lower triangle ==== 601* 602 CALL DCOPY( NS, V, LDV, WORK, 1 ) 603 BETA = WORK( 1 ) 604 CALL DLARFG( NS, BETA, WORK( 2 ), 1, TAU ) 605 WORK( 1 ) = ONE 606* 607 CALL DLASET( 'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT ) 608* 609 CALL DLARF( 'L', NS, JW, WORK, 1, TAU, T, LDT, 610 $ WORK( JW+1 ) ) 611 CALL DLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT, 612 $ WORK( JW+1 ) ) 613 CALL DLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV, 614 $ WORK( JW+1 ) ) 615* 616 CALL DGEHRD( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ), 617 $ LWORK-JW, INFO ) 618 END IF 619* 620* ==== Copy updated reduced window into place ==== 621* 622 IF( KWTOP.GT.1 ) 623 $ H( KWTOP, KWTOP-1 ) = S*V( 1, 1 ) 624 CALL DLACPY( 'U', JW, JW, T, LDT, H( KWTOP, KWTOP ), LDH ) 625 CALL DCOPY( JW-1, T( 2, 1 ), LDT+1, H( KWTOP+1, KWTOP ), 626 $ LDH+1 ) 627* 628* ==== Accumulate orthogonal matrix in order update 629* . H and Z, if requested. ==== 630* 631 IF( NS.GT.1 .AND. S.NE.ZERO ) 632 $ CALL DORMHR( 'R', 'N', JW, NS, 1, NS, T, LDT, WORK, V, LDV, 633 $ WORK( JW+1 ), LWORK-JW, INFO ) 634* 635* ==== Update vertical slab in H ==== 636* 637 IF( WANTT ) THEN 638 LTOP = 1 639 ELSE 640 LTOP = KTOP 641 END IF 642 DO 70 KROW = LTOP, KWTOP - 1, NV 643 KLN = MIN( NV, KWTOP-KROW ) 644 CALL DGEMM( 'N', 'N', KLN, JW, JW, ONE, H( KROW, KWTOP ), 645 $ LDH, V, LDV, ZERO, WV, LDWV ) 646 CALL DLACPY( 'A', KLN, JW, WV, LDWV, H( KROW, KWTOP ), LDH ) 647 70 CONTINUE 648* 649* ==== Update horizontal slab in H ==== 650* 651 IF( WANTT ) THEN 652 DO 80 KCOL = KBOT + 1, N, NH 653 KLN = MIN( NH, N-KCOL+1 ) 654 CALL DGEMM( 'C', 'N', JW, KLN, JW, ONE, V, LDV, 655 $ H( KWTOP, KCOL ), LDH, ZERO, T, LDT ) 656 CALL DLACPY( 'A', JW, KLN, T, LDT, H( KWTOP, KCOL ), 657 $ LDH ) 658 80 CONTINUE 659 END IF 660* 661* ==== Update vertical slab in Z ==== 662* 663 IF( WANTZ ) THEN 664 DO 90 KROW = ILOZ, IHIZ, NV 665 KLN = MIN( NV, IHIZ-KROW+1 ) 666 CALL DGEMM( 'N', 'N', KLN, JW, JW, ONE, Z( KROW, KWTOP ), 667 $ LDZ, V, LDV, ZERO, WV, LDWV ) 668 CALL DLACPY( 'A', KLN, JW, WV, LDWV, Z( KROW, KWTOP ), 669 $ LDZ ) 670 90 CONTINUE 671 END IF 672 END IF 673* 674* ==== Return the number of deflations ... ==== 675* 676 ND = JW - NS 677* 678* ==== ... and the number of shifts. (Subtracting 679* . INFQR from the spike length takes care 680* . of the case of a rare QR failure while 681* . calculating eigenvalues of the deflation 682* . window.) ==== 683* 684 NS = NS - INFQR 685* 686* ==== Return optimal workspace. ==== 687* 688 WORK( 1 ) = DBLE( LWKOPT ) 689* 690* ==== End of DLAQR3 ==== 691* 692 END 693