1*> \brief \b STREVC 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download STREVC + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/strevc.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/strevc.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/strevc.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE STREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, 22* LDVR, MM, M, WORK, INFO ) 23* 24* .. Scalar Arguments .. 25* CHARACTER HOWMNY, SIDE 26* INTEGER INFO, LDT, LDVL, LDVR, M, MM, N 27* .. 28* .. Array Arguments .. 29* LOGICAL SELECT( * ) 30* REAL T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), 31* $ WORK( * ) 32* .. 33* 34* 35*> \par Purpose: 36* ============= 37*> 38*> \verbatim 39*> 40*> STREVC computes some or all of the right and/or left eigenvectors of 41*> a real upper quasi-triangular matrix T. 42*> Matrices of this type are produced by the Schur factorization of 43*> a real general matrix: A = Q*T*Q**T, as computed by SHSEQR. 44*> 45*> The right eigenvector x and the left eigenvector y of T corresponding 46*> to an eigenvalue w are defined by: 47*> 48*> T*x = w*x, (y**H)*T = w*(y**H) 49*> 50*> where y**H denotes the conjugate transpose of y. 51*> The eigenvalues are not input to this routine, but are read directly 52*> from the diagonal blocks of T. 53*> 54*> This routine returns the matrices X and/or Y of right and left 55*> eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an 56*> input matrix. If Q is the orthogonal factor that reduces a matrix 57*> A to Schur form T, then Q*X and Q*Y are the matrices of right and 58*> left eigenvectors of A. 59*> \endverbatim 60* 61* Arguments: 62* ========== 63* 64*> \param[in] SIDE 65*> \verbatim 66*> SIDE is CHARACTER*1 67*> = 'R': compute right eigenvectors only; 68*> = 'L': compute left eigenvectors only; 69*> = 'B': compute both right and left eigenvectors. 70*> \endverbatim 71*> 72*> \param[in] HOWMNY 73*> \verbatim 74*> HOWMNY is CHARACTER*1 75*> = 'A': compute all right and/or left eigenvectors; 76*> = 'B': compute all right and/or left eigenvectors, 77*> backtransformed by the matrices in VR and/or VL; 78*> = 'S': compute selected right and/or left eigenvectors, 79*> as indicated by the logical array SELECT. 80*> \endverbatim 81*> 82*> \param[in,out] SELECT 83*> \verbatim 84*> SELECT is LOGICAL array, dimension (N) 85*> If HOWMNY = 'S', SELECT specifies the eigenvectors to be 86*> computed. 87*> If w(j) is a real eigenvalue, the corresponding real 88*> eigenvector is computed if SELECT(j) is .TRUE.. 89*> If w(j) and w(j+1) are the real and imaginary parts of a 90*> complex eigenvalue, the corresponding complex eigenvector is 91*> computed if either SELECT(j) or SELECT(j+1) is .TRUE., and 92*> on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to 93*> .FALSE.. 94*> Not referenced if HOWMNY = 'A' or 'B'. 95*> \endverbatim 96*> 97*> \param[in] N 98*> \verbatim 99*> N is INTEGER 100*> The order of the matrix T. N >= 0. 101*> \endverbatim 102*> 103*> \param[in] T 104*> \verbatim 105*> T is REAL array, dimension (LDT,N) 106*> The upper quasi-triangular matrix T in Schur canonical form. 107*> \endverbatim 108*> 109*> \param[in] LDT 110*> \verbatim 111*> LDT is INTEGER 112*> The leading dimension of the array T. LDT >= max(1,N). 113*> \endverbatim 114*> 115*> \param[in,out] VL 116*> \verbatim 117*> VL is REAL array, dimension (LDVL,MM) 118*> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must 119*> contain an N-by-N matrix Q (usually the orthogonal matrix Q 120*> of Schur vectors returned by SHSEQR). 121*> On exit, if SIDE = 'L' or 'B', VL contains: 122*> if HOWMNY = 'A', the matrix Y of left eigenvectors of T; 123*> if HOWMNY = 'B', the matrix Q*Y; 124*> if HOWMNY = 'S', the left eigenvectors of T specified by 125*> SELECT, stored consecutively in the columns 126*> of VL, in the same order as their 127*> eigenvalues. 128*> A complex eigenvector corresponding to a complex eigenvalue 129*> is stored in two consecutive columns, the first holding the 130*> real part, and the second the imaginary part. 131*> Not referenced if SIDE = 'R'. 132*> \endverbatim 133*> 134*> \param[in] LDVL 135*> \verbatim 136*> LDVL is INTEGER 137*> The leading dimension of the array VL. LDVL >= 1, and if 138*> SIDE = 'L' or 'B', LDVL >= N. 139*> \endverbatim 140*> 141*> \param[in,out] VR 142*> \verbatim 143*> VR is REAL array, dimension (LDVR,MM) 144*> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must 145*> contain an N-by-N matrix Q (usually the orthogonal matrix Q 146*> of Schur vectors returned by SHSEQR). 147*> On exit, if SIDE = 'R' or 'B', VR contains: 148*> if HOWMNY = 'A', the matrix X of right eigenvectors of T; 149*> if HOWMNY = 'B', the matrix Q*X; 150*> if HOWMNY = 'S', the right eigenvectors of T specified by 151*> SELECT, stored consecutively in the columns 152*> of VR, in the same order as their 153*> eigenvalues. 154*> A complex eigenvector corresponding to a complex eigenvalue 155*> is stored in two consecutive columns, the first holding the 156*> real part and the second the imaginary part. 157*> Not referenced if SIDE = 'L'. 158*> \endverbatim 159*> 160*> \param[in] LDVR 161*> \verbatim 162*> LDVR is INTEGER 163*> The leading dimension of the array VR. LDVR >= 1, and if 164*> SIDE = 'R' or 'B', LDVR >= N. 165*> \endverbatim 166*> 167*> \param[in] MM 168*> \verbatim 169*> MM is INTEGER 170*> The number of columns in the arrays VL and/or VR. MM >= M. 171*> \endverbatim 172*> 173*> \param[out] M 174*> \verbatim 175*> M is INTEGER 176*> The number of columns in the arrays VL and/or VR actually 177*> used to store the eigenvectors. 178*> If HOWMNY = 'A' or 'B', M is set to N. 179*> Each selected real eigenvector occupies one column and each 180*> selected complex eigenvector occupies two columns. 181*> \endverbatim 182*> 183*> \param[out] WORK 184*> \verbatim 185*> WORK is REAL array, dimension (3*N) 186*> \endverbatim 187*> 188*> \param[out] INFO 189*> \verbatim 190*> INFO is INTEGER 191*> = 0: successful exit 192*> < 0: if INFO = -i, the i-th argument had an illegal value 193*> \endverbatim 194* 195* Authors: 196* ======== 197* 198*> \author Univ. of Tennessee 199*> \author Univ. of California Berkeley 200*> \author Univ. of Colorado Denver 201*> \author NAG Ltd. 202* 203*> \ingroup realOTHERcomputational 204* 205*> \par Further Details: 206* ===================== 207*> 208*> \verbatim 209*> 210*> The algorithm used in this program is basically backward (forward) 211*> substitution, with scaling to make the the code robust against 212*> possible overflow. 213*> 214*> Each eigenvector is normalized so that the element of largest 215*> magnitude has magnitude 1; here the magnitude of a complex number 216*> (x,y) is taken to be |x| + |y|. 217*> \endverbatim 218*> 219* ===================================================================== 220 SUBROUTINE STREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, 221 $ LDVR, MM, M, WORK, INFO ) 222* 223* -- LAPACK computational routine -- 224* -- LAPACK is a software package provided by Univ. of Tennessee, -- 225* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 226* 227* .. Scalar Arguments .. 228 CHARACTER HOWMNY, SIDE 229 INTEGER INFO, LDT, LDVL, LDVR, M, MM, N 230* .. 231* .. Array Arguments .. 232 LOGICAL SELECT( * ) 233 REAL T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), 234 $ WORK( * ) 235* .. 236* 237* ===================================================================== 238* 239* .. Parameters .. 240 REAL ZERO, ONE 241 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 242* .. 243* .. Local Scalars .. 244 LOGICAL ALLV, BOTHV, LEFTV, OVER, PAIR, RIGHTV, SOMEV 245 INTEGER I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI, N2 246 REAL BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE, 247 $ SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR, 248 $ XNORM 249* .. 250* .. External Functions .. 251 LOGICAL LSAME 252 INTEGER ISAMAX 253 REAL SDOT, SLAMCH 254 EXTERNAL LSAME, ISAMAX, SDOT, SLAMCH 255* .. 256* .. External Subroutines .. 257 EXTERNAL SAXPY, SCOPY, SGEMV, SLABAD, SLALN2, SSCAL, 258 $ XERBLA 259* .. 260* .. Intrinsic Functions .. 261 INTRINSIC ABS, MAX, SQRT 262* .. 263* .. Local Arrays .. 264 REAL X( 2, 2 ) 265* .. 266* .. Executable Statements .. 267* 268* Decode and test the input parameters 269* 270 BOTHV = LSAME( SIDE, 'B' ) 271 RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV 272 LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV 273* 274 ALLV = LSAME( HOWMNY, 'A' ) 275 OVER = LSAME( HOWMNY, 'B' ) 276 SOMEV = LSAME( HOWMNY, 'S' ) 277* 278 INFO = 0 279 IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN 280 INFO = -1 281 ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN 282 INFO = -2 283 ELSE IF( N.LT.0 ) THEN 284 INFO = -4 285 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN 286 INFO = -6 287 ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN 288 INFO = -8 289 ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN 290 INFO = -10 291 ELSE 292* 293* Set M to the number of columns required to store the selected 294* eigenvectors, standardize the array SELECT if necessary, and 295* test MM. 296* 297 IF( SOMEV ) THEN 298 M = 0 299 PAIR = .FALSE. 300 DO 10 J = 1, N 301 IF( PAIR ) THEN 302 PAIR = .FALSE. 303 SELECT( J ) = .FALSE. 304 ELSE 305 IF( J.LT.N ) THEN 306 IF( T( J+1, J ).EQ.ZERO ) THEN 307 IF( SELECT( J ) ) 308 $ M = M + 1 309 ELSE 310 PAIR = .TRUE. 311 IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN 312 SELECT( J ) = .TRUE. 313 M = M + 2 314 END IF 315 END IF 316 ELSE 317 IF( SELECT( N ) ) 318 $ M = M + 1 319 END IF 320 END IF 321 10 CONTINUE 322 ELSE 323 M = N 324 END IF 325* 326 IF( MM.LT.M ) THEN 327 INFO = -11 328 END IF 329 END IF 330 IF( INFO.NE.0 ) THEN 331 CALL XERBLA( 'STREVC', -INFO ) 332 RETURN 333 END IF 334* 335* Quick return if possible. 336* 337 IF( N.EQ.0 ) 338 $ RETURN 339* 340* Set the constants to control overflow. 341* 342 UNFL = SLAMCH( 'Safe minimum' ) 343 OVFL = ONE / UNFL 344 CALL SLABAD( UNFL, OVFL ) 345 ULP = SLAMCH( 'Precision' ) 346 SMLNUM = UNFL*( N / ULP ) 347 BIGNUM = ( ONE-ULP ) / SMLNUM 348* 349* Compute 1-norm of each column of strictly upper triangular 350* part of T to control overflow in triangular solver. 351* 352 WORK( 1 ) = ZERO 353 DO 30 J = 2, N 354 WORK( J ) = ZERO 355 DO 20 I = 1, J - 1 356 WORK( J ) = WORK( J ) + ABS( T( I, J ) ) 357 20 CONTINUE 358 30 CONTINUE 359* 360* Index IP is used to specify the real or complex eigenvalue: 361* IP = 0, real eigenvalue, 362* 1, first of conjugate complex pair: (wr,wi) 363* -1, second of conjugate complex pair: (wr,wi) 364* 365 N2 = 2*N 366* 367 IF( RIGHTV ) THEN 368* 369* Compute right eigenvectors. 370* 371 IP = 0 372 IS = M 373 DO 140 KI = N, 1, -1 374* 375 IF( IP.EQ.1 ) 376 $ GO TO 130 377 IF( KI.EQ.1 ) 378 $ GO TO 40 379 IF( T( KI, KI-1 ).EQ.ZERO ) 380 $ GO TO 40 381 IP = -1 382* 383 40 CONTINUE 384 IF( SOMEV ) THEN 385 IF( IP.EQ.0 ) THEN 386 IF( .NOT.SELECT( KI ) ) 387 $ GO TO 130 388 ELSE 389 IF( .NOT.SELECT( KI-1 ) ) 390 $ GO TO 130 391 END IF 392 END IF 393* 394* Compute the KI-th eigenvalue (WR,WI). 395* 396 WR = T( KI, KI ) 397 WI = ZERO 398 IF( IP.NE.0 ) 399 $ WI = SQRT( ABS( T( KI, KI-1 ) ) )* 400 $ SQRT( ABS( T( KI-1, KI ) ) ) 401 SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM ) 402* 403 IF( IP.EQ.0 ) THEN 404* 405* Real right eigenvector 406* 407 WORK( KI+N ) = ONE 408* 409* Form right-hand side 410* 411 DO 50 K = 1, KI - 1 412 WORK( K+N ) = -T( K, KI ) 413 50 CONTINUE 414* 415* Solve the upper quasi-triangular system: 416* (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK. 417* 418 JNXT = KI - 1 419 DO 60 J = KI - 1, 1, -1 420 IF( J.GT.JNXT ) 421 $ GO TO 60 422 J1 = J 423 J2 = J 424 JNXT = J - 1 425 IF( J.GT.1 ) THEN 426 IF( T( J, J-1 ).NE.ZERO ) THEN 427 J1 = J - 1 428 JNXT = J - 2 429 END IF 430 END IF 431* 432 IF( J1.EQ.J2 ) THEN 433* 434* 1-by-1 diagonal block 435* 436 CALL SLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ), 437 $ LDT, ONE, ONE, WORK( J+N ), N, WR, 438 $ ZERO, X, 2, SCALE, XNORM, IERR ) 439* 440* Scale X(1,1) to avoid overflow when updating 441* the right-hand side. 442* 443 IF( XNORM.GT.ONE ) THEN 444 IF( WORK( J ).GT.BIGNUM / XNORM ) THEN 445 X( 1, 1 ) = X( 1, 1 ) / XNORM 446 SCALE = SCALE / XNORM 447 END IF 448 END IF 449* 450* Scale if necessary 451* 452 IF( SCALE.NE.ONE ) 453 $ CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 ) 454 WORK( J+N ) = X( 1, 1 ) 455* 456* Update right-hand side 457* 458 CALL SAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1, 459 $ WORK( 1+N ), 1 ) 460* 461 ELSE 462* 463* 2-by-2 diagonal block 464* 465 CALL SLALN2( .FALSE., 2, 1, SMIN, ONE, 466 $ T( J-1, J-1 ), LDT, ONE, ONE, 467 $ WORK( J-1+N ), N, WR, ZERO, X, 2, 468 $ SCALE, XNORM, IERR ) 469* 470* Scale X(1,1) and X(2,1) to avoid overflow when 471* updating the right-hand side. 472* 473 IF( XNORM.GT.ONE ) THEN 474 BETA = MAX( WORK( J-1 ), WORK( J ) ) 475 IF( BETA.GT.BIGNUM / XNORM ) THEN 476 X( 1, 1 ) = X( 1, 1 ) / XNORM 477 X( 2, 1 ) = X( 2, 1 ) / XNORM 478 SCALE = SCALE / XNORM 479 END IF 480 END IF 481* 482* Scale if necessary 483* 484 IF( SCALE.NE.ONE ) 485 $ CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 ) 486 WORK( J-1+N ) = X( 1, 1 ) 487 WORK( J+N ) = X( 2, 1 ) 488* 489* Update right-hand side 490* 491 CALL SAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1, 492 $ WORK( 1+N ), 1 ) 493 CALL SAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1, 494 $ WORK( 1+N ), 1 ) 495 END IF 496 60 CONTINUE 497* 498* Copy the vector x or Q*x to VR and normalize. 499* 500 IF( .NOT.OVER ) THEN 501 CALL SCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 ) 502* 503 II = ISAMAX( KI, VR( 1, IS ), 1 ) 504 REMAX = ONE / ABS( VR( II, IS ) ) 505 CALL SSCAL( KI, REMAX, VR( 1, IS ), 1 ) 506* 507 DO 70 K = KI + 1, N 508 VR( K, IS ) = ZERO 509 70 CONTINUE 510 ELSE 511 IF( KI.GT.1 ) 512 $ CALL SGEMV( 'N', N, KI-1, ONE, VR, LDVR, 513 $ WORK( 1+N ), 1, WORK( KI+N ), 514 $ VR( 1, KI ), 1 ) 515* 516 II = ISAMAX( N, VR( 1, KI ), 1 ) 517 REMAX = ONE / ABS( VR( II, KI ) ) 518 CALL SSCAL( N, REMAX, VR( 1, KI ), 1 ) 519 END IF 520* 521 ELSE 522* 523* Complex right eigenvector. 524* 525* Initial solve 526* [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0. 527* [ (T(KI,KI-1) T(KI,KI) ) ] 528* 529 IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN 530 WORK( KI-1+N ) = ONE 531 WORK( KI+N2 ) = WI / T( KI-1, KI ) 532 ELSE 533 WORK( KI-1+N ) = -WI / T( KI, KI-1 ) 534 WORK( KI+N2 ) = ONE 535 END IF 536 WORK( KI+N ) = ZERO 537 WORK( KI-1+N2 ) = ZERO 538* 539* Form right-hand side 540* 541 DO 80 K = 1, KI - 2 542 WORK( K+N ) = -WORK( KI-1+N )*T( K, KI-1 ) 543 WORK( K+N2 ) = -WORK( KI+N2 )*T( K, KI ) 544 80 CONTINUE 545* 546* Solve upper quasi-triangular system: 547* (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2) 548* 549 JNXT = KI - 2 550 DO 90 J = KI - 2, 1, -1 551 IF( J.GT.JNXT ) 552 $ GO TO 90 553 J1 = J 554 J2 = J 555 JNXT = J - 1 556 IF( J.GT.1 ) THEN 557 IF( T( J, J-1 ).NE.ZERO ) THEN 558 J1 = J - 1 559 JNXT = J - 2 560 END IF 561 END IF 562* 563 IF( J1.EQ.J2 ) THEN 564* 565* 1-by-1 diagonal block 566* 567 CALL SLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ), 568 $ LDT, ONE, ONE, WORK( J+N ), N, WR, WI, 569 $ X, 2, SCALE, XNORM, IERR ) 570* 571* Scale X(1,1) and X(1,2) to avoid overflow when 572* updating the right-hand side. 573* 574 IF( XNORM.GT.ONE ) THEN 575 IF( WORK( J ).GT.BIGNUM / XNORM ) THEN 576 X( 1, 1 ) = X( 1, 1 ) / XNORM 577 X( 1, 2 ) = X( 1, 2 ) / XNORM 578 SCALE = SCALE / XNORM 579 END IF 580 END IF 581* 582* Scale if necessary 583* 584 IF( SCALE.NE.ONE ) THEN 585 CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 ) 586 CALL SSCAL( KI, SCALE, WORK( 1+N2 ), 1 ) 587 END IF 588 WORK( J+N ) = X( 1, 1 ) 589 WORK( J+N2 ) = X( 1, 2 ) 590* 591* Update the right-hand side 592* 593 CALL SAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1, 594 $ WORK( 1+N ), 1 ) 595 CALL SAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1, 596 $ WORK( 1+N2 ), 1 ) 597* 598 ELSE 599* 600* 2-by-2 diagonal block 601* 602 CALL SLALN2( .FALSE., 2, 2, SMIN, ONE, 603 $ T( J-1, J-1 ), LDT, ONE, ONE, 604 $ WORK( J-1+N ), N, WR, WI, X, 2, SCALE, 605 $ XNORM, IERR ) 606* 607* Scale X to avoid overflow when updating 608* the right-hand side. 609* 610 IF( XNORM.GT.ONE ) THEN 611 BETA = MAX( WORK( J-1 ), WORK( J ) ) 612 IF( BETA.GT.BIGNUM / XNORM ) THEN 613 REC = ONE / XNORM 614 X( 1, 1 ) = X( 1, 1 )*REC 615 X( 1, 2 ) = X( 1, 2 )*REC 616 X( 2, 1 ) = X( 2, 1 )*REC 617 X( 2, 2 ) = X( 2, 2 )*REC 618 SCALE = SCALE*REC 619 END IF 620 END IF 621* 622* Scale if necessary 623* 624 IF( SCALE.NE.ONE ) THEN 625 CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 ) 626 CALL SSCAL( KI, SCALE, WORK( 1+N2 ), 1 ) 627 END IF 628 WORK( J-1+N ) = X( 1, 1 ) 629 WORK( J+N ) = X( 2, 1 ) 630 WORK( J-1+N2 ) = X( 1, 2 ) 631 WORK( J+N2 ) = X( 2, 2 ) 632* 633* Update the right-hand side 634* 635 CALL SAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1, 636 $ WORK( 1+N ), 1 ) 637 CALL SAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1, 638 $ WORK( 1+N ), 1 ) 639 CALL SAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1, 640 $ WORK( 1+N2 ), 1 ) 641 CALL SAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1, 642 $ WORK( 1+N2 ), 1 ) 643 END IF 644 90 CONTINUE 645* 646* Copy the vector x or Q*x to VR and normalize. 647* 648 IF( .NOT.OVER ) THEN 649 CALL SCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 ) 650 CALL SCOPY( KI, WORK( 1+N2 ), 1, VR( 1, IS ), 1 ) 651* 652 EMAX = ZERO 653 DO 100 K = 1, KI 654 EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+ 655 $ ABS( VR( K, IS ) ) ) 656 100 CONTINUE 657* 658 REMAX = ONE / EMAX 659 CALL SSCAL( KI, REMAX, VR( 1, IS-1 ), 1 ) 660 CALL SSCAL( KI, REMAX, VR( 1, IS ), 1 ) 661* 662 DO 110 K = KI + 1, N 663 VR( K, IS-1 ) = ZERO 664 VR( K, IS ) = ZERO 665 110 CONTINUE 666* 667 ELSE 668* 669 IF( KI.GT.2 ) THEN 670 CALL SGEMV( 'N', N, KI-2, ONE, VR, LDVR, 671 $ WORK( 1+N ), 1, WORK( KI-1+N ), 672 $ VR( 1, KI-1 ), 1 ) 673 CALL SGEMV( 'N', N, KI-2, ONE, VR, LDVR, 674 $ WORK( 1+N2 ), 1, WORK( KI+N2 ), 675 $ VR( 1, KI ), 1 ) 676 ELSE 677 CALL SSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 ) 678 CALL SSCAL( N, WORK( KI+N2 ), VR( 1, KI ), 1 ) 679 END IF 680* 681 EMAX = ZERO 682 DO 120 K = 1, N 683 EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+ 684 $ ABS( VR( K, KI ) ) ) 685 120 CONTINUE 686 REMAX = ONE / EMAX 687 CALL SSCAL( N, REMAX, VR( 1, KI-1 ), 1 ) 688 CALL SSCAL( N, REMAX, VR( 1, KI ), 1 ) 689 END IF 690 END IF 691* 692 IS = IS - 1 693 IF( IP.NE.0 ) 694 $ IS = IS - 1 695 130 CONTINUE 696 IF( IP.EQ.1 ) 697 $ IP = 0 698 IF( IP.EQ.-1 ) 699 $ IP = 1 700 140 CONTINUE 701 END IF 702* 703 IF( LEFTV ) THEN 704* 705* Compute left eigenvectors. 706* 707 IP = 0 708 IS = 1 709 DO 260 KI = 1, N 710* 711 IF( IP.EQ.-1 ) 712 $ GO TO 250 713 IF( KI.EQ.N ) 714 $ GO TO 150 715 IF( T( KI+1, KI ).EQ.ZERO ) 716 $ GO TO 150 717 IP = 1 718* 719 150 CONTINUE 720 IF( SOMEV ) THEN 721 IF( .NOT.SELECT( KI ) ) 722 $ GO TO 250 723 END IF 724* 725* Compute the KI-th eigenvalue (WR,WI). 726* 727 WR = T( KI, KI ) 728 WI = ZERO 729 IF( IP.NE.0 ) 730 $ WI = SQRT( ABS( T( KI, KI+1 ) ) )* 731 $ SQRT( ABS( T( KI+1, KI ) ) ) 732 SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM ) 733* 734 IF( IP.EQ.0 ) THEN 735* 736* Real left eigenvector. 737* 738 WORK( KI+N ) = ONE 739* 740* Form right-hand side 741* 742 DO 160 K = KI + 1, N 743 WORK( K+N ) = -T( KI, K ) 744 160 CONTINUE 745* 746* Solve the quasi-triangular system: 747* (T(KI+1:N,KI+1:N) - WR)**T*X = SCALE*WORK 748* 749 VMAX = ONE 750 VCRIT = BIGNUM 751* 752 JNXT = KI + 1 753 DO 170 J = KI + 1, N 754 IF( J.LT.JNXT ) 755 $ GO TO 170 756 J1 = J 757 J2 = J 758 JNXT = J + 1 759 IF( J.LT.N ) THEN 760 IF( T( J+1, J ).NE.ZERO ) THEN 761 J2 = J + 1 762 JNXT = J + 2 763 END IF 764 END IF 765* 766 IF( J1.EQ.J2 ) THEN 767* 768* 1-by-1 diagonal block 769* 770* Scale if necessary to avoid overflow when forming 771* the right-hand side. 772* 773 IF( WORK( J ).GT.VCRIT ) THEN 774 REC = ONE / VMAX 775 CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) 776 VMAX = ONE 777 VCRIT = BIGNUM 778 END IF 779* 780 WORK( J+N ) = WORK( J+N ) - 781 $ SDOT( J-KI-1, T( KI+1, J ), 1, 782 $ WORK( KI+1+N ), 1 ) 783* 784* Solve (T(J,J)-WR)**T*X = WORK 785* 786 CALL SLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ), 787 $ LDT, ONE, ONE, WORK( J+N ), N, WR, 788 $ ZERO, X, 2, SCALE, XNORM, IERR ) 789* 790* Scale if necessary 791* 792 IF( SCALE.NE.ONE ) 793 $ CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) 794 WORK( J+N ) = X( 1, 1 ) 795 VMAX = MAX( ABS( WORK( J+N ) ), VMAX ) 796 VCRIT = BIGNUM / VMAX 797* 798 ELSE 799* 800* 2-by-2 diagonal block 801* 802* Scale if necessary to avoid overflow when forming 803* the right-hand side. 804* 805 BETA = MAX( WORK( J ), WORK( J+1 ) ) 806 IF( BETA.GT.VCRIT ) THEN 807 REC = ONE / VMAX 808 CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) 809 VMAX = ONE 810 VCRIT = BIGNUM 811 END IF 812* 813 WORK( J+N ) = WORK( J+N ) - 814 $ SDOT( J-KI-1, T( KI+1, J ), 1, 815 $ WORK( KI+1+N ), 1 ) 816* 817 WORK( J+1+N ) = WORK( J+1+N ) - 818 $ SDOT( J-KI-1, T( KI+1, J+1 ), 1, 819 $ WORK( KI+1+N ), 1 ) 820* 821* Solve 822* [T(J,J)-WR T(J,J+1) ]**T* X = SCALE*( WORK1 ) 823* [T(J+1,J) T(J+1,J+1)-WR] ( WORK2 ) 824* 825 CALL SLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ), 826 $ LDT, ONE, ONE, WORK( J+N ), N, WR, 827 $ ZERO, X, 2, SCALE, XNORM, IERR ) 828* 829* Scale if necessary 830* 831 IF( SCALE.NE.ONE ) 832 $ CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) 833 WORK( J+N ) = X( 1, 1 ) 834 WORK( J+1+N ) = X( 2, 1 ) 835* 836 VMAX = MAX( ABS( WORK( J+N ) ), 837 $ ABS( WORK( J+1+N ) ), VMAX ) 838 VCRIT = BIGNUM / VMAX 839* 840 END IF 841 170 CONTINUE 842* 843* Copy the vector x or Q*x to VL and normalize. 844* 845 IF( .NOT.OVER ) THEN 846 CALL SCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 ) 847* 848 II = ISAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1 849 REMAX = ONE / ABS( VL( II, IS ) ) 850 CALL SSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 ) 851* 852 DO 180 K = 1, KI - 1 853 VL( K, IS ) = ZERO 854 180 CONTINUE 855* 856 ELSE 857* 858 IF( KI.LT.N ) 859 $ CALL SGEMV( 'N', N, N-KI, ONE, VL( 1, KI+1 ), LDVL, 860 $ WORK( KI+1+N ), 1, WORK( KI+N ), 861 $ VL( 1, KI ), 1 ) 862* 863 II = ISAMAX( N, VL( 1, KI ), 1 ) 864 REMAX = ONE / ABS( VL( II, KI ) ) 865 CALL SSCAL( N, REMAX, VL( 1, KI ), 1 ) 866* 867 END IF 868* 869 ELSE 870* 871* Complex left eigenvector. 872* 873* Initial solve: 874* ((T(KI,KI) T(KI,KI+1) )**T - (WR - I* WI))*X = 0. 875* ((T(KI+1,KI) T(KI+1,KI+1)) ) 876* 877 IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN 878 WORK( KI+N ) = WI / T( KI, KI+1 ) 879 WORK( KI+1+N2 ) = ONE 880 ELSE 881 WORK( KI+N ) = ONE 882 WORK( KI+1+N2 ) = -WI / T( KI+1, KI ) 883 END IF 884 WORK( KI+1+N ) = ZERO 885 WORK( KI+N2 ) = ZERO 886* 887* Form right-hand side 888* 889 DO 190 K = KI + 2, N 890 WORK( K+N ) = -WORK( KI+N )*T( KI, K ) 891 WORK( K+N2 ) = -WORK( KI+1+N2 )*T( KI+1, K ) 892 190 CONTINUE 893* 894* Solve complex quasi-triangular system: 895* ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2 896* 897 VMAX = ONE 898 VCRIT = BIGNUM 899* 900 JNXT = KI + 2 901 DO 200 J = KI + 2, N 902 IF( J.LT.JNXT ) 903 $ GO TO 200 904 J1 = J 905 J2 = J 906 JNXT = J + 1 907 IF( J.LT.N ) THEN 908 IF( T( J+1, J ).NE.ZERO ) THEN 909 J2 = J + 1 910 JNXT = J + 2 911 END IF 912 END IF 913* 914 IF( J1.EQ.J2 ) THEN 915* 916* 1-by-1 diagonal block 917* 918* Scale if necessary to avoid overflow when 919* forming the right-hand side elements. 920* 921 IF( WORK( J ).GT.VCRIT ) THEN 922 REC = ONE / VMAX 923 CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) 924 CALL SSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 ) 925 VMAX = ONE 926 VCRIT = BIGNUM 927 END IF 928* 929 WORK( J+N ) = WORK( J+N ) - 930 $ SDOT( J-KI-2, T( KI+2, J ), 1, 931 $ WORK( KI+2+N ), 1 ) 932 WORK( J+N2 ) = WORK( J+N2 ) - 933 $ SDOT( J-KI-2, T( KI+2, J ), 1, 934 $ WORK( KI+2+N2 ), 1 ) 935* 936* Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2 937* 938 CALL SLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ), 939 $ LDT, ONE, ONE, WORK( J+N ), N, WR, 940 $ -WI, X, 2, SCALE, XNORM, IERR ) 941* 942* Scale if necessary 943* 944 IF( SCALE.NE.ONE ) THEN 945 CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) 946 CALL SSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 ) 947 END IF 948 WORK( J+N ) = X( 1, 1 ) 949 WORK( J+N2 ) = X( 1, 2 ) 950 VMAX = MAX( ABS( WORK( J+N ) ), 951 $ ABS( WORK( J+N2 ) ), VMAX ) 952 VCRIT = BIGNUM / VMAX 953* 954 ELSE 955* 956* 2-by-2 diagonal block 957* 958* Scale if necessary to avoid overflow when forming 959* the right-hand side elements. 960* 961 BETA = MAX( WORK( J ), WORK( J+1 ) ) 962 IF( BETA.GT.VCRIT ) THEN 963 REC = ONE / VMAX 964 CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) 965 CALL SSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 ) 966 VMAX = ONE 967 VCRIT = BIGNUM 968 END IF 969* 970 WORK( J+N ) = WORK( J+N ) - 971 $ SDOT( J-KI-2, T( KI+2, J ), 1, 972 $ WORK( KI+2+N ), 1 ) 973* 974 WORK( J+N2 ) = WORK( J+N2 ) - 975 $ SDOT( J-KI-2, T( KI+2, J ), 1, 976 $ WORK( KI+2+N2 ), 1 ) 977* 978 WORK( J+1+N ) = WORK( J+1+N ) - 979 $ SDOT( J-KI-2, T( KI+2, J+1 ), 1, 980 $ WORK( KI+2+N ), 1 ) 981* 982 WORK( J+1+N2 ) = WORK( J+1+N2 ) - 983 $ SDOT( J-KI-2, T( KI+2, J+1 ), 1, 984 $ WORK( KI+2+N2 ), 1 ) 985* 986* Solve 2-by-2 complex linear equation 987* ([T(j,j) T(j,j+1) ]**T-(wr-i*wi)*I)*X = SCALE*B 988* ([T(j+1,j) T(j+1,j+1)] ) 989* 990 CALL SLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ), 991 $ LDT, ONE, ONE, WORK( J+N ), N, WR, 992 $ -WI, X, 2, SCALE, XNORM, IERR ) 993* 994* Scale if necessary 995* 996 IF( SCALE.NE.ONE ) THEN 997 CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) 998 CALL SSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 ) 999 END IF 1000 WORK( J+N ) = X( 1, 1 ) 1001 WORK( J+N2 ) = X( 1, 2 ) 1002 WORK( J+1+N ) = X( 2, 1 ) 1003 WORK( J+1+N2 ) = X( 2, 2 ) 1004 VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ), 1005 $ ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ), VMAX ) 1006 VCRIT = BIGNUM / VMAX 1007* 1008 END IF 1009 200 CONTINUE 1010* 1011* Copy the vector x or Q*x to VL and normalize. 1012* 1013 IF( .NOT.OVER ) THEN 1014 CALL SCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 ) 1015 CALL SCOPY( N-KI+1, WORK( KI+N2 ), 1, VL( KI, IS+1 ), 1016 $ 1 ) 1017* 1018 EMAX = ZERO 1019 DO 220 K = KI, N 1020 EMAX = MAX( EMAX, ABS( VL( K, IS ) )+ 1021 $ ABS( VL( K, IS+1 ) ) ) 1022 220 CONTINUE 1023 REMAX = ONE / EMAX 1024 CALL SSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 ) 1025 CALL SSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 ) 1026* 1027 DO 230 K = 1, KI - 1 1028 VL( K, IS ) = ZERO 1029 VL( K, IS+1 ) = ZERO 1030 230 CONTINUE 1031 ELSE 1032 IF( KI.LT.N-1 ) THEN 1033 CALL SGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ), 1034 $ LDVL, WORK( KI+2+N ), 1, WORK( KI+N ), 1035 $ VL( 1, KI ), 1 ) 1036 CALL SGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ), 1037 $ LDVL, WORK( KI+2+N2 ), 1, 1038 $ WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 ) 1039 ELSE 1040 CALL SSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 ) 1041 CALL SSCAL( N, WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 ) 1042 END IF 1043* 1044 EMAX = ZERO 1045 DO 240 K = 1, N 1046 EMAX = MAX( EMAX, ABS( VL( K, KI ) )+ 1047 $ ABS( VL( K, KI+1 ) ) ) 1048 240 CONTINUE 1049 REMAX = ONE / EMAX 1050 CALL SSCAL( N, REMAX, VL( 1, KI ), 1 ) 1051 CALL SSCAL( N, REMAX, VL( 1, KI+1 ), 1 ) 1052* 1053 END IF 1054* 1055 END IF 1056* 1057 IS = IS + 1 1058 IF( IP.NE.0 ) 1059 $ IS = IS + 1 1060 250 CONTINUE 1061 IF( IP.EQ.-1 ) 1062 $ IP = 0 1063 IF( IP.EQ.1 ) 1064 $ IP = -1 1065* 1066 260 CONTINUE 1067* 1068 END IF 1069* 1070 RETURN 1071* 1072* End of STREVC 1073* 1074 END 1075