1*> \brief \b SLAQR2 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 SLAQR2 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaqr2.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaqr2.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaqr2.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SLAQR2( 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* REAL 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*> SLAQR2 is identical to SLAQR3 except that it avoids 43*> recursion by calling SLAHQR instead of SLAQR4. 44*> 45*> Aggressive early deflation: 46*> 47*> This subroutine accepts as input an upper Hessenberg matrix 48*> H and performs an orthogonal similarity transformation 49*> designed to detect and deflate fully converged eigenvalues from 50*> a trailing principal submatrix. On output H has been over- 51*> written by a new Hessenberg matrix that is a perturbation of 52*> an orthogonal similarity transformation of H. It is to be 53*> hoped that the final version of H has many zero subdiagonal 54*> entries. 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 quasi-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 orthogonal matrix Z is updated so 74*> so that the orthogonal 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 orthogonal 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 .LE. NW .LE. (KBOT-KTOP+1). 107*> \endverbatim 108*> 109*> \param[in,out] H 110*> \verbatim 111*> H is REAL 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 an orthogonal 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 .LE. 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 .LE. ILOZ .LE. IHIZ .LE. N. 137*> \endverbatim 138*> 139*> \param[in,out] Z 140*> \verbatim 141*> Z is REAL array, dimension (LDZ,N) 142*> IF WANTZ is .TRUE., then on output, the orthogonal 143*> similarity transformation mentioned above has been 144*> accumulated into Z(ILOZ:IHIZ,ILO:IHI) 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 .LE. 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] SR 171*> \verbatim 172*> SR is REAL array, dimension KBOT 173*> \endverbatim 174*> 175*> \param[out] SI 176*> \verbatim 177*> SI is REAL array, dimension KBOT 178*> On output, the real and imaginary parts of approximate 179*> eigenvalues that may be used for shifts are stored in 180*> SR(KBOT-ND-NS+1) through SR(KBOT-ND) and 181*> SI(KBOT-ND-NS+1) through SI(KBOT-ND), respectively. 182*> The real and imaginary parts of converged eigenvalues 183*> are stored in SR(KBOT-ND+1) through SR(KBOT) and 184*> SI(KBOT-ND+1) through SI(KBOT), respectively. 185*> \endverbatim 186*> 187*> \param[out] V 188*> \verbatim 189*> V is REAL array, dimension (LDV,NW) 190*> An NW-by-NW work array. 191*> \endverbatim 192*> 193*> \param[in] LDV 194*> \verbatim 195*> LDV is integer scalar 196*> The leading dimension of V just as declared in the 197*> calling subroutine. NW .LE. LDV 198*> \endverbatim 199*> 200*> \param[in] NH 201*> \verbatim 202*> NH is integer scalar 203*> The number of columns of T. NH.GE.NW. 204*> \endverbatim 205*> 206*> \param[out] T 207*> \verbatim 208*> T is REAL array, dimension (LDT,NW) 209*> \endverbatim 210*> 211*> \param[in] LDT 212*> \verbatim 213*> LDT is integer 214*> The leading dimension of T just as declared in the 215*> calling subroutine. NW .LE. LDT 216*> \endverbatim 217*> 218*> \param[in] NV 219*> \verbatim 220*> NV is integer 221*> The number of rows of work array WV available for 222*> workspace. NV.GE.NW. 223*> \endverbatim 224*> 225*> \param[out] WV 226*> \verbatim 227*> WV is REAL array, dimension (LDWV,NW) 228*> \endverbatim 229*> 230*> \param[in] LDWV 231*> \verbatim 232*> LDWV is integer 233*> The leading dimension of W just as declared in the 234*> calling subroutine. NW .LE. LDV 235*> \endverbatim 236*> 237*> \param[out] WORK 238*> \verbatim 239*> WORK is REAL array, dimension LWORK. 240*> On exit, WORK(1) is set to an estimate of the optimal value 241*> of LWORK for the given values of N, NW, KTOP and KBOT. 242*> \endverbatim 243*> 244*> \param[in] LWORK 245*> \verbatim 246*> LWORK is integer 247*> The dimension of the work array WORK. LWORK = 2*NW 248*> suffices, but greater efficiency may result from larger 249*> values of LWORK. 250*> 251*> If LWORK = -1, then a workspace query is assumed; SLAQR2 252*> only estimates the optimal workspace size for the given 253*> values of N, NW, KTOP and KBOT. The estimate is returned 254*> in WORK(1). No error message related to LWORK is issued 255*> by XERBLA. Neither H nor Z are accessed. 256*> \endverbatim 257* 258* Authors: 259* ======== 260* 261*> \author Univ. of Tennessee 262*> \author Univ. of California Berkeley 263*> \author Univ. of Colorado Denver 264*> \author NAG Ltd. 265* 266*> \date September 2012 267* 268*> \ingroup realOTHERauxiliary 269* 270*> \par Contributors: 271* ================== 272*> 273*> Karen Braman and Ralph Byers, Department of Mathematics, 274*> University of Kansas, USA 275*> 276* ===================================================================== 277 SUBROUTINE SLAQR2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, 278 $ IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T, 279 $ LDT, NV, WV, LDWV, WORK, LWORK ) 280* 281* -- LAPACK auxiliary routine (version 3.4.2) -- 282* -- LAPACK is a software package provided by Univ. of Tennessee, -- 283* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 284* September 2012 285* 286* .. Scalar Arguments .. 287 INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV, 288 $ LDZ, LWORK, N, ND, NH, NS, NV, NW 289 LOGICAL WANTT, WANTZ 290* .. 291* .. Array Arguments .. 292 REAL H( LDH, * ), SI( * ), SR( * ), T( LDT, * ), 293 $ V( LDV, * ), WORK( * ), WV( LDWV, * ), 294 $ Z( LDZ, * ) 295* .. 296* 297* ================================================================ 298* .. Parameters .. 299 REAL ZERO, ONE 300 PARAMETER ( ZERO = 0.0e0, ONE = 1.0e0 ) 301* .. 302* .. Local Scalars .. 303 REAL AA, BB, BETA, CC, CS, DD, EVI, EVK, FOO, S, 304 $ SAFMAX, SAFMIN, SMLNUM, SN, TAU, ULP 305 INTEGER I, IFST, ILST, INFO, INFQR, J, JW, K, KCOL, 306 $ KEND, KLN, KROW, KWTOP, LTOP, LWK1, LWK2, 307 $ LWKOPT 308 LOGICAL BULGE, SORTED 309* .. 310* .. External Functions .. 311 REAL SLAMCH 312 EXTERNAL SLAMCH 313* .. 314* .. External Subroutines .. 315 EXTERNAL SCOPY, SGEHRD, SGEMM, SLABAD, SLACPY, SLAHQR, 316 $ SLANV2, SLARF, SLARFG, SLASET, SORMHR, STREXC 317* .. 318* .. Intrinsic Functions .. 319 INTRINSIC ABS, INT, MAX, MIN, REAL, SQRT 320* .. 321* .. Executable Statements .. 322* 323* ==== Estimate optimal workspace. ==== 324* 325 JW = MIN( NW, KBOT-KTOP+1 ) 326 IF( JW.LE.2 ) THEN 327 LWKOPT = 1 328 ELSE 329* 330* ==== Workspace query call to SGEHRD ==== 331* 332 CALL SGEHRD( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO ) 333 LWK1 = INT( WORK( 1 ) ) 334* 335* ==== Workspace query call to SORMHR ==== 336* 337 CALL SORMHR( 'R', 'N', JW, JW, 1, JW-1, T, LDT, WORK, V, LDV, 338 $ WORK, -1, INFO ) 339 LWK2 = INT( WORK( 1 ) ) 340* 341* ==== Optimal workspace ==== 342* 343 LWKOPT = JW + MAX( LWK1, LWK2 ) 344 END IF 345* 346* ==== Quick return in case of workspace query. ==== 347* 348 IF( LWORK.EQ.-1 ) THEN 349 WORK( 1 ) = REAL( LWKOPT ) 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 = ONE / 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 SR( KWTOP ) = H( KWTOP, KWTOP ) 387 SI( KWTOP ) = ZERO 388 NS = 1 389 ND = 0 390 IF( ABS( S ).LE.MAX( SMLNUM, ULP*ABS( H( KWTOP, KWTOP ) ) ) ) 391 $ THEN 392 NS = 0 393 ND = 1 394 IF( KWTOP.GT.KTOP ) 395 $ H( KWTOP, KWTOP-1 ) = ZERO 396 END IF 397 WORK( 1 ) = ONE 398 RETURN 399 END IF 400* 401* ==== Convert to spike-triangular form. (In case of a 402* . rare QR failure, this routine continues to do 403* . aggressive early deflation using that part of 404* . the deflation window that converged using INFQR 405* . here and there to keep track.) ==== 406* 407 CALL SLACPY( 'U', JW, JW, H( KWTOP, KWTOP ), LDH, T, LDT ) 408 CALL SCOPY( JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ), LDT+1 ) 409* 410 CALL SLASET( 'A', JW, JW, ZERO, ONE, V, LDV ) 411 CALL SLAHQR( .true., .true., JW, 1, JW, T, LDT, SR( KWTOP ), 412 $ SI( KWTOP ), 1, JW, V, LDV, INFQR ) 413* 414* ==== STREXC needs a clean margin near the diagonal ==== 415* 416 DO 10 J = 1, JW - 3 417 T( J+2, J ) = ZERO 418 T( J+3, J ) = ZERO 419 10 CONTINUE 420 IF( JW.GT.2 ) 421 $ T( JW, JW-2 ) = ZERO 422* 423* ==== Deflation detection loop ==== 424* 425 NS = JW 426 ILST = INFQR + 1 427 20 CONTINUE 428 IF( ILST.LE.NS ) THEN 429 IF( NS.EQ.1 ) THEN 430 BULGE = .FALSE. 431 ELSE 432 BULGE = T( NS, NS-1 ).NE.ZERO 433 END IF 434* 435* ==== Small spike tip test for deflation ==== 436* 437 IF( .NOT.BULGE ) THEN 438* 439* ==== Real eigenvalue ==== 440* 441 FOO = ABS( T( NS, NS ) ) 442 IF( FOO.EQ.ZERO ) 443 $ FOO = ABS( S ) 444 IF( ABS( S*V( 1, NS ) ).LE.MAX( SMLNUM, ULP*FOO ) ) THEN 445* 446* ==== Deflatable ==== 447* 448 NS = NS - 1 449 ELSE 450* 451* ==== Undeflatable. Move it up out of the way. 452* . (STREXC can not fail in this case.) ==== 453* 454 IFST = NS 455 CALL STREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK, 456 $ INFO ) 457 ILST = ILST + 1 458 END IF 459 ELSE 460* 461* ==== Complex conjugate pair ==== 462* 463 FOO = ABS( T( NS, NS ) ) + SQRT( ABS( T( NS, NS-1 ) ) )* 464 $ SQRT( ABS( T( NS-1, NS ) ) ) 465 IF( FOO.EQ.ZERO ) 466 $ FOO = ABS( S ) 467 IF( MAX( ABS( S*V( 1, NS ) ), ABS( S*V( 1, NS-1 ) ) ).LE. 468 $ MAX( SMLNUM, ULP*FOO ) ) THEN 469* 470* ==== Deflatable ==== 471* 472 NS = NS - 2 473 ELSE 474* 475* ==== Undeflatable. Move them up out of the way. 476* . Fortunately, STREXC does the right thing with 477* . ILST in case of a rare exchange failure. ==== 478* 479 IFST = NS 480 CALL STREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK, 481 $ INFO ) 482 ILST = ILST + 2 483 END IF 484 END IF 485* 486* ==== End deflation detection loop ==== 487* 488 GO TO 20 489 END IF 490* 491* ==== Return to Hessenberg form ==== 492* 493 IF( NS.EQ.0 ) 494 $ S = ZERO 495* 496 IF( NS.LT.JW ) THEN 497* 498* ==== sorting diagonal blocks of T improves accuracy for 499* . graded matrices. Bubble sort deals well with 500* . exchange failures. ==== 501* 502 SORTED = .false. 503 I = NS + 1 504 30 CONTINUE 505 IF( SORTED ) 506 $ GO TO 50 507 SORTED = .true. 508* 509 KEND = I - 1 510 I = INFQR + 1 511 IF( I.EQ.NS ) THEN 512 K = I + 1 513 ELSE IF( T( I+1, I ).EQ.ZERO ) THEN 514 K = I + 1 515 ELSE 516 K = I + 2 517 END IF 518 40 CONTINUE 519 IF( K.LE.KEND ) THEN 520 IF( K.EQ.I+1 ) THEN 521 EVI = ABS( T( I, I ) ) 522 ELSE 523 EVI = ABS( T( I, I ) ) + SQRT( ABS( T( I+1, I ) ) )* 524 $ SQRT( ABS( T( I, I+1 ) ) ) 525 END IF 526* 527 IF( K.EQ.KEND ) THEN 528 EVK = ABS( T( K, K ) ) 529 ELSE IF( T( K+1, K ).EQ.ZERO ) THEN 530 EVK = ABS( T( K, K ) ) 531 ELSE 532 EVK = ABS( T( K, K ) ) + SQRT( ABS( T( K+1, K ) ) )* 533 $ SQRT( ABS( T( K, K+1 ) ) ) 534 END IF 535* 536 IF( EVI.GE.EVK ) THEN 537 I = K 538 ELSE 539 SORTED = .false. 540 IFST = I 541 ILST = K 542 CALL STREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK, 543 $ INFO ) 544 IF( INFO.EQ.0 ) THEN 545 I = ILST 546 ELSE 547 I = K 548 END IF 549 END IF 550 IF( I.EQ.KEND ) THEN 551 K = I + 1 552 ELSE IF( T( I+1, I ).EQ.ZERO ) THEN 553 K = I + 1 554 ELSE 555 K = I + 2 556 END IF 557 GO TO 40 558 END IF 559 GO TO 30 560 50 CONTINUE 561 END IF 562* 563* ==== Restore shift/eigenvalue array from T ==== 564* 565 I = JW 566 60 CONTINUE 567 IF( I.GE.INFQR+1 ) THEN 568 IF( I.EQ.INFQR+1 ) THEN 569 SR( KWTOP+I-1 ) = T( I, I ) 570 SI( KWTOP+I-1 ) = ZERO 571 I = I - 1 572 ELSE IF( T( I, I-1 ).EQ.ZERO ) THEN 573 SR( KWTOP+I-1 ) = T( I, I ) 574 SI( KWTOP+I-1 ) = ZERO 575 I = I - 1 576 ELSE 577 AA = T( I-1, I-1 ) 578 CC = T( I, I-1 ) 579 BB = T( I-1, I ) 580 DD = T( I, I ) 581 CALL SLANV2( AA, BB, CC, DD, SR( KWTOP+I-2 ), 582 $ SI( KWTOP+I-2 ), SR( KWTOP+I-1 ), 583 $ SI( KWTOP+I-1 ), CS, SN ) 584 I = I - 2 585 END IF 586 GO TO 60 587 END IF 588* 589 IF( NS.LT.JW .OR. S.EQ.ZERO ) THEN 590 IF( NS.GT.1 .AND. S.NE.ZERO ) THEN 591* 592* ==== Reflect spike back into lower triangle ==== 593* 594 CALL SCOPY( NS, V, LDV, WORK, 1 ) 595 BETA = WORK( 1 ) 596 CALL SLARFG( NS, BETA, WORK( 2 ), 1, TAU ) 597 WORK( 1 ) = ONE 598* 599 CALL SLASET( 'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT ) 600* 601 CALL SLARF( 'L', NS, JW, WORK, 1, TAU, T, LDT, 602 $ WORK( JW+1 ) ) 603 CALL SLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT, 604 $ WORK( JW+1 ) ) 605 CALL SLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV, 606 $ WORK( JW+1 ) ) 607* 608 CALL SGEHRD( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ), 609 $ LWORK-JW, INFO ) 610 END IF 611* 612* ==== Copy updated reduced window into place ==== 613* 614 IF( KWTOP.GT.1 ) 615 $ H( KWTOP, KWTOP-1 ) = S*V( 1, 1 ) 616 CALL SLACPY( 'U', JW, JW, T, LDT, H( KWTOP, KWTOP ), LDH ) 617 CALL SCOPY( JW-1, T( 2, 1 ), LDT+1, H( KWTOP+1, KWTOP ), 618 $ LDH+1 ) 619* 620* ==== Accumulate orthogonal matrix in order update 621* . H and Z, if requested. ==== 622* 623 IF( NS.GT.1 .AND. S.NE.ZERO ) 624 $ CALL SORMHR( 'R', 'N', JW, NS, 1, NS, T, LDT, WORK, V, LDV, 625 $ WORK( JW+1 ), LWORK-JW, INFO ) 626* 627* ==== Update vertical slab in H ==== 628* 629 IF( WANTT ) THEN 630 LTOP = 1 631 ELSE 632 LTOP = KTOP 633 END IF 634 DO 70 KROW = LTOP, KWTOP - 1, NV 635 KLN = MIN( NV, KWTOP-KROW ) 636 CALL SGEMM( 'N', 'N', KLN, JW, JW, ONE, H( KROW, KWTOP ), 637 $ LDH, V, LDV, ZERO, WV, LDWV ) 638 CALL SLACPY( 'A', KLN, JW, WV, LDWV, H( KROW, KWTOP ), LDH ) 639 70 CONTINUE 640* 641* ==== Update horizontal slab in H ==== 642* 643 IF( WANTT ) THEN 644 DO 80 KCOL = KBOT + 1, N, NH 645 KLN = MIN( NH, N-KCOL+1 ) 646 CALL SGEMM( 'C', 'N', JW, KLN, JW, ONE, V, LDV, 647 $ H( KWTOP, KCOL ), LDH, ZERO, T, LDT ) 648 CALL SLACPY( 'A', JW, KLN, T, LDT, H( KWTOP, KCOL ), 649 $ LDH ) 650 80 CONTINUE 651 END IF 652* 653* ==== Update vertical slab in Z ==== 654* 655 IF( WANTZ ) THEN 656 DO 90 KROW = ILOZ, IHIZ, NV 657 KLN = MIN( NV, IHIZ-KROW+1 ) 658 CALL SGEMM( 'N', 'N', KLN, JW, JW, ONE, Z( KROW, KWTOP ), 659 $ LDZ, V, LDV, ZERO, WV, LDWV ) 660 CALL SLACPY( 'A', KLN, JW, WV, LDWV, Z( KROW, KWTOP ), 661 $ LDZ ) 662 90 CONTINUE 663 END IF 664 END IF 665* 666* ==== Return the number of deflations ... ==== 667* 668 ND = JW - NS 669* 670* ==== ... and the number of shifts. (Subtracting 671* . INFQR from the spike length takes care 672* . of the case of a rare QR failure while 673* . calculating eigenvalues of the deflation 674* . window.) ==== 675* 676 NS = NS - INFQR 677* 678* ==== Return optimal workspace. ==== 679* 680 WORK( 1 ) = REAL( LWKOPT ) 681* 682* ==== End of SLAQR2 ==== 683* 684 END 685c $Id$ 686