SRC/strevc.f

   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**T)*T = w*(y**T)
  49 *>
  50 *> where y**T denotes the 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 *> \date November 2011
 204 *
 205 *> \ingroup realOTHERcomputational
 206 *
 207 *> \par Further Details:
 208 *  =====================
 209 *>
 210 *> \verbatim
 211 *>
 212 *>  The algorithm used in this program is basically backward (forward)
 213 *>  substitution, with scaling to make the the code robust against
 214 *>  possible overflow.
 215 *>
 216 *>  Each eigenvector is normalized so that the element of largest
 217 *>  magnitude has magnitude 1; here the magnitude of a complex number
 218 *>  (x,y) is taken to be |x| + |y|.
 219 *> \endverbatim
 220 *>
 221 *  =====================================================================
 222       SUBROUTINE STREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
 223      $                   LDVR, MM, M, WORK, INFO )
 224 *
 225 *  -- LAPACK computational routine (version 3.4.0) --
 226 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 227 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 228 *     November 2011
 229 *
 230 *     .. Scalar Arguments ..
 231       CHARACTER          HOWMNY, SIDE
 232       INTEGER            INFO, LDT, LDVL, LDVR, M, MM, N
 233 *     ..
 234 *     .. Array Arguments ..
 235       LOGICAL            SELECT( * )
 236       REAL               T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
 237      $                   WORK( * )
 238 *     ..
 239 *
 240 *  =====================================================================
 241 *
 242 *     .. Parameters ..
 243       REAL               ZERO, ONE
 244       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 245 *     ..
 246 *     .. Local Scalars ..
 247       LOGICAL            ALLV, BOTHV, LEFTV, OVER, PAIR, RIGHTV, SOMEV
 248       INTEGER            I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI, N2
 249       REAL               BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
 250      $                   SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR,
 251      $                   XNORM
 252 *     ..
 253 *     .. External Functions ..
 254       LOGICAL            LSAME
 255       INTEGER            ISAMAX
 256       REAL               SDOT, SLAMCH
 257       EXTERNAL           LSAME, ISAMAX, SDOT, SLAMCH
 258 *     ..
 259 *     .. External Subroutines ..
 260       EXTERNAL           SAXPY, SCOPY, SGEMV, SLABAD, SLALN2, SSCAL,
 261      $                   XERBLA
 262 *     ..
 263 *     .. Intrinsic Functions ..
 264       INTRINSIC          ABS, MAX, SQRT
 265 *     ..
 266 *     .. Local Arrays ..
 267       REAL               X( 2, 2 )
 268 *     ..
 269 *     .. Executable Statements ..
 270 *
 271 *     Decode and test the input parameters
 272 *
 273       BOTHV = LSAME( SIDE, 'B' )
 274       RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
 275       LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
 276 *
 277       ALLV = LSAME( HOWMNY, 'A' )
 278       OVER = LSAME( HOWMNY, 'B' )
 279       SOMEV = LSAME( HOWMNY, 'S' )
 280 *
 281       INFO = 0
 282       IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
 283          INFO = -1
 284       ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
 285          INFO = -2
 286       ELSE IF( N.LT.0 ) THEN
 287          INFO = -4
 288       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
 289          INFO = -6
 290       ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
 291          INFO = -8
 292       ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
 293          INFO = -10
 294       ELSE
 295 *
 296 *        Set M to the number of columns required to store the selected
 297 *        eigenvectors, standardize the array SELECT if necessary, and
 298 *        test MM.
 299 *
 300          IF( SOMEV ) THEN
 301             M = 0
 302             PAIR = .FALSE.
 303             DO 10 J = 1, N
 304                IF( PAIR ) THEN
 305                   PAIR = .FALSE.
 306                   SELECT( J ) = .FALSE.
 307                ELSE
 308                   IF( J.LT.N ) THEN
 309                      IF( T( J+1, J ).EQ.ZERO ) THEN
 310                         IF( SELECT( J ) )
 311      $                     M = M + 1
 312                      ELSE
 313                         PAIR = .TRUE.
 314                         IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN
 315                            SELECT( J ) = .TRUE.
 316                            M = M + 2
 317                         END IF
 318                      END IF
 319                   ELSE
 320                      IF( SELECT( N ) )
 321      $                  M = M + 1
 322                   END IF
 323                END IF
 324    10       CONTINUE
 325          ELSE
 326             M = N
 327          END IF
 328 *
 329          IF( MM.LT.M ) THEN
 330             INFO = -11
 331          END IF
 332       END IF
 333       IF( INFO.NE.0 ) THEN
 334          CALL XERBLA( 'STREVC', -INFO )
 335          RETURN
 336       END IF
 337 *
 338 *     Quick return if possible.
 339 *
 340       IF( N.EQ.0 )
 341      $   RETURN
 342 *
 343 *     Set the constants to control overflow.
 344 *
 345       UNFL = SLAMCH( 'Safe minimum' )
 346       OVFL = ONE / UNFL
 347       CALL SLABAD( UNFL, OVFL )
 348       ULP = SLAMCH( 'Precision' )
 349       SMLNUM = UNFL*( N / ULP )
 350       BIGNUM = ( ONE-ULP ) / SMLNUM
 351 *
 352 *     Compute 1-norm of each column of strictly upper triangular
 353 *     part of T to control overflow in triangular solver.
 354 *
 355       WORK( 1 ) = ZERO
 356       DO 30 J = 2, N
 357          WORK( J ) = ZERO
 358          DO 20 I = 1, J - 1
 359             WORK( J ) = WORK( J ) + ABS( T( I, J ) )
 360    20    CONTINUE
 361    30 CONTINUE
 362 *
 363 *     Index IP is used to specify the real or complex eigenvalue:
 364 *       IP = 0, real eigenvalue,
 365 *            1, first of conjugate complex pair: (wr,wi)
 366 *           -1, second of conjugate complex pair: (wr,wi)
 367 *
 368       N2 = 2*N
 369 *
 370       IF( RIGHTV ) THEN
 371 *
 372 *        Compute right eigenvectors.
 373 *
 374          IP = 0
 375          IS = M
 376          DO 140 KI = N, 1, -1
 377 *
 378             IF( IP.EQ.1 )
 379      $         GO TO 130
 380             IF( KI.EQ.1 )
 381      $         GO TO 40
 382             IF( T( KI, KI-1 ).EQ.ZERO )
 383      $         GO TO 40
 384             IP = -1
 385 *
 386    40       CONTINUE
 387             IF( SOMEV ) THEN
 388                IF( IP.EQ.0 ) THEN
 389                   IF( .NOT.SELECT( KI ) )
 390      $               GO TO 130
 391                ELSE
 392                   IF( .NOT.SELECT( KI-1 ) )
 393      $               GO TO 130
 394                END IF
 395             END IF
 396 *
 397 *           Compute the KI-th eigenvalue (WR,WI).
 398 *
 399             WR = T( KI, KI )
 400             WI = ZERO
 401             IF( IP.NE.0 )
 402      $         WI = SQRT( ABS( T( KI, KI-1 ) ) )*
 403      $              SQRT( ABS( T( KI-1, KI ) ) )
 404             SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
 405 *
 406             IF( IP.EQ.0 ) THEN
 407 *
 408 *              Real right eigenvector
 409 *
 410                WORK( KI+N ) = ONE
 411 *
 412 *              Form right-hand side
 413 *
 414                DO 50 K = 1, KI - 1
 415                   WORK( K+N ) = -T( K, KI )
 416    50          CONTINUE
 417 *
 418 *              Solve the upper quasi-triangular system:
 419 *                 (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK.
 420 *
 421                JNXT = KI - 1
 422                DO 60 J = KI - 1, 1, -1
 423                   IF( J.GT.JNXT )
 424      $               GO TO 60
 425                   J1 = J
 426                   J2 = J
 427                   JNXT = J - 1
 428                   IF( J.GT.1 ) THEN
 429                      IF( T( J, J-1 ).NE.ZERO ) THEN
 430                         J1 = J - 1
 431                         JNXT = J - 2
 432                      END IF
 433                   END IF
 434 *
 435                   IF( J1.EQ.J2 ) THEN
 436 *
 437 *                    1-by-1 diagonal block
 438 *
 439                      CALL SLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
 440      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
 441      $                            ZERO, X, 2, SCALE, XNORM, IERR )
 442 *
 443 *                    Scale X(1,1) to avoid overflow when updating
 444 *                    the right-hand side.
 445 *
 446                      IF( XNORM.GT.ONE ) THEN
 447                         IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
 448                            X( 1, 1 ) = X( 1, 1 ) / XNORM
 449                            SCALE = SCALE / XNORM
 450                         END IF
 451                      END IF
 452 *
 453 *                    Scale if necessary
 454 *
 455                      IF( SCALE.NE.ONE )
 456      $                  CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 )
 457                      WORK( J+N ) = X( 1, 1 )
 458 *
 459 *                    Update right-hand side
 460 *
 461                      CALL SAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
 462      $                           WORK( 1+N ), 1 )
 463 *
 464                   ELSE
 465 *
 466 *                    2-by-2 diagonal block
 467 *
 468                      CALL SLALN2( .FALSE., 2, 1, SMIN, ONE,
 469      $                            T( J-1, J-1 ), LDT, ONE, ONE,
 470      $                            WORK( J-1+N ), N, WR, ZERO, X, 2,
 471      $                            SCALE, XNORM, IERR )
 472 *
 473 *                    Scale X(1,1) and X(2,1) to avoid overflow when
 474 *                    updating the right-hand side.
 475 *
 476                      IF( XNORM.GT.ONE ) THEN
 477                         BETA = MAX( WORK( J-1 ), WORK( J ) )
 478                         IF( BETA.GT.BIGNUM / XNORM ) THEN
 479                            X( 1, 1 ) = X( 1, 1 ) / XNORM
 480                            X( 2, 1 ) = X( 2, 1 ) / XNORM
 481                            SCALE = SCALE / XNORM
 482                         END IF
 483                      END IF
 484 *
 485 *                    Scale if necessary
 486 *
 487                      IF( SCALE.NE.ONE )
 488      $                  CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 )
 489                      WORK( J-1+N ) = X( 1, 1 )
 490                      WORK( J+N ) = X( 2, 1 )
 491 *
 492 *                    Update right-hand side
 493 *
 494                      CALL SAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
 495      $                           WORK( 1+N ), 1 )
 496                      CALL SAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
 497      $                           WORK( 1+N ), 1 )
 498                   END IF
 499    60          CONTINUE
 500 *
 501 *              Copy the vector x or Q*x to VR and normalize.
 502 *
 503                IF( .NOT.OVER ) THEN
 504                   CALL SCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 )
 505 *
 506                   II = ISAMAX( KI, VR( 1, IS ), 1 )
 507                   REMAX = ONE / ABS( VR( II, IS ) )
 508                   CALL SSCAL( KI, REMAX, VR( 1, IS ), 1 )
 509 *
 510                   DO 70 K = KI + 1, N
 511                      VR( K, IS ) = ZERO
 512    70             CONTINUE
 513                ELSE
 514                   IF( KI.GT.1 )
 515      $               CALL SGEMV( 'N', N, KI-1, ONE, VR, LDVR,
 516      $                           WORK( 1+N ), 1, WORK( KI+N ),
 517      $                           VR( 1, KI ), 1 )
 518 *
 519                   II = ISAMAX( N, VR( 1, KI ), 1 )
 520                   REMAX = ONE / ABS( VR( II, KI ) )
 521                   CALL SSCAL( N, REMAX, VR( 1, KI ), 1 )
 522                END IF
 523 *
 524             ELSE
 525 *
 526 *              Complex right eigenvector.
 527 *
 528 *              Initial solve
 529 *                [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.
 530 *                [ (T(KI,KI-1)   T(KI,KI)   )               ]
 531 *
 532                IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN
 533                   WORK( KI-1+N ) = ONE
 534                   WORK( KI+N2 ) = WI / T( KI-1, KI )
 535                ELSE
 536                   WORK( KI-1+N ) = -WI / T( KI, KI-1 )
 537                   WORK( KI+N2 ) = ONE
 538                END IF
 539                WORK( KI+N ) = ZERO
 540                WORK( KI-1+N2 ) = ZERO
 541 *
 542 *              Form right-hand side
 543 *
 544                DO 80 K = 1, KI - 2
 545                   WORK( K+N ) = -WORK( KI-1+N )*T( K, KI-1 )
 546                   WORK( K+N2 ) = -WORK( KI+N2 )*T( K, KI )
 547    80          CONTINUE
 548 *
 549 *              Solve upper quasi-triangular system:
 550 *              (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2)
 551 *
 552                JNXT = KI - 2
 553                DO 90 J = KI - 2, 1, -1
 554                   IF( J.GT.JNXT )
 555      $               GO TO 90
 556                   J1 = J
 557                   J2 = J
 558                   JNXT = J - 1
 559                   IF( J.GT.1 ) THEN
 560                      IF( T( J, J-1 ).NE.ZERO ) THEN
 561                         J1 = J - 1
 562                         JNXT = J - 2
 563                      END IF
 564                   END IF
 565 *
 566                   IF( J1.EQ.J2 ) THEN
 567 *
 568 *                    1-by-1 diagonal block
 569 *
 570                      CALL SLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
 571      $                            LDT, ONE, ONE, WORK( J+N ), N, WR, WI,
 572      $                            X, 2, SCALE, XNORM, IERR )
 573 *
 574 *                    Scale X(1,1) and X(1,2) to avoid overflow when
 575 *                    updating the right-hand side.
 576 *
 577                      IF( XNORM.GT.ONE ) THEN
 578                         IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
 579                            X( 1, 1 ) = X( 1, 1 ) / XNORM
 580                            X( 1, 2 ) = X( 1, 2 ) / XNORM
 581                            SCALE = SCALE / XNORM
 582                         END IF
 583                      END IF
 584 *
 585 *                    Scale if necessary
 586 *
 587                      IF( SCALE.NE.ONE ) THEN
 588                         CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 )
 589                         CALL SSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
 590                      END IF
 591                      WORK( J+N ) = X( 1, 1 )
 592                      WORK( J+N2 ) = X( 1, 2 )
 593 *
 594 *                    Update the right-hand side
 595 *
 596                      CALL SAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
 597      $                           WORK( 1+N ), 1 )
 598                      CALL SAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1,
 599      $                           WORK( 1+N2 ), 1 )
 600 *
 601                   ELSE
 602 *
 603 *                    2-by-2 diagonal block
 604 *
 605                      CALL SLALN2( .FALSE., 2, 2, SMIN, ONE,
 606      $                            T( J-1, J-1 ), LDT, ONE, ONE,
 607      $                            WORK( J-1+N ), N, WR, WI, X, 2, SCALE,
 608      $                            XNORM, IERR )
 609 *
 610 *                    Scale X to avoid overflow when updating
 611 *                    the right-hand side.
 612 *
 613                      IF( XNORM.GT.ONE ) THEN
 614                         BETA = MAX( WORK( J-1 ), WORK( J ) )
 615                         IF( BETA.GT.BIGNUM / XNORM ) THEN
 616                            REC = ONE / XNORM
 617                            X( 1, 1 ) = X( 1, 1 )*REC
 618                            X( 1, 2 ) = X( 1, 2 )*REC
 619                            X( 2, 1 ) = X( 2, 1 )*REC
 620                            X( 2, 2 ) = X( 2, 2 )*REC
 621                            SCALE = SCALE*REC
 622                         END IF
 623                      END IF
 624 *
 625 *                    Scale if necessary
 626 *
 627                      IF( SCALE.NE.ONE ) THEN
 628                         CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 )
 629                         CALL SSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
 630                      END IF
 631                      WORK( J-1+N ) = X( 1, 1 )
 632                      WORK( J+N ) = X( 2, 1 )
 633                      WORK( J-1+N2 ) = X( 1, 2 )
 634                      WORK( J+N2 ) = X( 2, 2 )
 635 *
 636 *                    Update the right-hand side
 637 *
 638                      CALL SAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
 639      $                           WORK( 1+N ), 1 )
 640                      CALL SAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
 641      $                           WORK( 1+N ), 1 )
 642                      CALL SAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
 643      $                           WORK( 1+N2 ), 1 )
 644                      CALL SAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1,
 645      $                           WORK( 1+N2 ), 1 )
 646                   END IF
 647    90          CONTINUE
 648 *
 649 *              Copy the vector x or Q*x to VR and normalize.
 650 *
 651                IF( .NOT.OVER ) THEN
 652                   CALL SCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 )
 653                   CALL SCOPY( KI, WORK( 1+N2 ), 1, VR( 1, IS ), 1 )
 654 *
 655                   EMAX = ZERO
 656                   DO 100 K = 1, KI
 657                      EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+
 658      $                      ABS( VR( K, IS ) ) )
 659   100             CONTINUE
 660 *
 661                   REMAX = ONE / EMAX
 662                   CALL SSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
 663                   CALL SSCAL( KI, REMAX, VR( 1, IS ), 1 )
 664 *
 665                   DO 110 K = KI + 1, N
 666                      VR( K, IS-1 ) = ZERO
 667                      VR( K, IS ) = ZERO
 668   110             CONTINUE
 669 *
 670                ELSE
 671 *
 672                   IF( KI.GT.2 ) THEN
 673                      CALL SGEMV( 'N', N, KI-2, ONE, VR, LDVR,
 674      $                           WORK( 1+N ), 1, WORK( KI-1+N ),
 675      $                           VR( 1, KI-1 ), 1 )
 676                      CALL SGEMV( 'N', N, KI-2, ONE, VR, LDVR,
 677      $                           WORK( 1+N2 ), 1, WORK( KI+N2 ),
 678      $                           VR( 1, KI ), 1 )
 679                   ELSE
 680                      CALL SSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 )
 681                      CALL SSCAL( N, WORK( KI+N2 ), VR( 1, KI ), 1 )
 682                   END IF
 683 *
 684                   EMAX = ZERO
 685                   DO 120 K = 1, N
 686                      EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+
 687      $                      ABS( VR( K, KI ) ) )
 688   120             CONTINUE
 689                   REMAX = ONE / EMAX
 690                   CALL SSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
 691                   CALL SSCAL( N, REMAX, VR( 1, KI ), 1 )
 692                END IF
 693             END IF
 694 *
 695             IS = IS - 1
 696             IF( IP.NE.0 )
 697      $         IS = IS - 1
 698   130       CONTINUE
 699             IF( IP.EQ.1 )
 700      $         IP = 0
 701             IF( IP.EQ.-1 )
 702      $         IP = 1
 703   140    CONTINUE
 704       END IF
 705 *
 706       IF( LEFTV ) THEN
 707 *
 708 *        Compute left eigenvectors.
 709 *
 710          IP = 0
 711          IS = 1
 712          DO 260 KI = 1, N
 713 *
 714             IF( IP.EQ.-1 )
 715      $         GO TO 250
 716             IF( KI.EQ.N )
 717      $         GO TO 150
 718             IF( T( KI+1, KI ).EQ.ZERO )
 719      $         GO TO 150
 720             IP = 1
 721 *
 722   150       CONTINUE
 723             IF( SOMEV ) THEN
 724                IF( .NOT.SELECT( KI ) )
 725      $            GO TO 250
 726             END IF
 727 *
 728 *           Compute the KI-th eigenvalue (WR,WI).
 729 *
 730             WR = T( KI, KI )
 731             WI = ZERO
 732             IF( IP.NE.0 )
 733      $         WI = SQRT( ABS( T( KI, KI+1 ) ) )*
 734      $              SQRT( ABS( T( KI+1, KI ) ) )
 735             SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
 736 *
 737             IF( IP.EQ.0 ) THEN
 738 *
 739 *              Real left eigenvector.
 740 *
 741                WORK( KI+N ) = ONE
 742 *
 743 *              Form right-hand side
 744 *
 745                DO 160 K = KI + 1, N
 746                   WORK( K+N ) = -T( KI, K )
 747   160          CONTINUE
 748 *
 749 *              Solve the quasi-triangular system:
 750 *                 (T(KI+1:N,KI+1:N) - WR)**T*X = SCALE*WORK
 751 *
 752                VMAX = ONE
 753                VCRIT = BIGNUM
 754 *
 755                JNXT = KI + 1
 756                DO 170 J = KI + 1, N
 757                   IF( J.LT.JNXT )
 758      $               GO TO 170
 759                   J1 = J
 760                   J2 = J
 761                   JNXT = J + 1
 762                   IF( J.LT.N ) THEN
 763                      IF( T( J+1, J ).NE.ZERO ) THEN
 764                         J2 = J + 1
 765                         JNXT = J + 2
 766                      END IF
 767                   END IF
 768 *
 769                   IF( J1.EQ.J2 ) THEN
 770 *
 771 *                    1-by-1 diagonal block
 772 *
 773 *                    Scale if necessary to avoid overflow when forming
 774 *                    the right-hand side.
 775 *
 776                      IF( WORK( J ).GT.VCRIT ) THEN
 777                         REC = ONE / VMAX
 778                         CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
 779                         VMAX = ONE
 780                         VCRIT = BIGNUM
 781                      END IF
 782 *
 783                      WORK( J+N ) = WORK( J+N ) -
 784      $                             SDOT( J-KI-1, T( KI+1, J ), 1,
 785      $                             WORK( KI+1+N ), 1 )
 786 *
 787 *                    Solve (T(J,J)-WR)**T*X = WORK
 788 *
 789                      CALL SLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
 790      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
 791      $                            ZERO, X, 2, SCALE, XNORM, IERR )
 792 *
 793 *                    Scale if necessary
 794 *
 795                      IF( SCALE.NE.ONE )
 796      $                  CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
 797                      WORK( J+N ) = X( 1, 1 )
 798                      VMAX = MAX( ABS( WORK( J+N ) ), VMAX )
 799                      VCRIT = BIGNUM / VMAX
 800 *
 801                   ELSE
 802 *
 803 *                    2-by-2 diagonal block
 804 *
 805 *                    Scale if necessary to avoid overflow when forming
 806 *                    the right-hand side.
 807 *
 808                      BETA = MAX( WORK( J ), WORK( J+1 ) )
 809                      IF( BETA.GT.VCRIT ) THEN
 810                         REC = ONE / VMAX
 811                         CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
 812                         VMAX = ONE
 813                         VCRIT = BIGNUM
 814                      END IF
 815 *
 816                      WORK( J+N ) = WORK( J+N ) -
 817      $                             SDOT( J-KI-1, T( KI+1, J ), 1,
 818      $                             WORK( KI+1+N ), 1 )
 819 *
 820                      WORK( J+1+N ) = WORK( J+1+N ) -
 821      $                               SDOT( J-KI-1, T( KI+1, J+1 ), 1,
 822      $                               WORK( KI+1+N ), 1 )
 823 *
 824 *                    Solve
 825 *                      [T(J,J)-WR   T(J,J+1)     ]**T* X = SCALE*( WORK1 )
 826 *                      [T(J+1,J)    T(J+1,J+1)-WR]               ( WORK2 )
 827 *
 828                      CALL SLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ),
 829      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
 830      $                            ZERO, X, 2, SCALE, XNORM, IERR )
 831 *
 832 *                    Scale if necessary
 833 *
 834                      IF( SCALE.NE.ONE )
 835      $                  CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
 836                      WORK( J+N ) = X( 1, 1 )
 837                      WORK( J+1+N ) = X( 2, 1 )
 838 *
 839                      VMAX = MAX( ABS( WORK( J+N ) ),
 840      $                      ABS( WORK( J+1+N ) ), VMAX )
 841                      VCRIT = BIGNUM / VMAX
 842 *
 843                   END IF
 844   170          CONTINUE
 845 *
 846 *              Copy the vector x or Q*x to VL and normalize.
 847 *
 848                IF( .NOT.OVER ) THEN
 849                   CALL SCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
 850 *
 851                   II = ISAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
 852                   REMAX = ONE / ABS( VL( II, IS ) )
 853                   CALL SSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
 854 *
 855                   DO 180 K = 1, KI - 1
 856                      VL( K, IS ) = ZERO
 857   180             CONTINUE
 858 *
 859                ELSE
 860 *
 861                   IF( KI.LT.N )
 862      $               CALL SGEMV( 'N', N, N-KI, ONE, VL( 1, KI+1 ), LDVL,
 863      $                           WORK( KI+1+N ), 1, WORK( KI+N ),
 864      $                           VL( 1, KI ), 1 )
 865 *
 866                   II = ISAMAX( N, VL( 1, KI ), 1 )
 867                   REMAX = ONE / ABS( VL( II, KI ) )
 868                   CALL SSCAL( N, REMAX, VL( 1, KI ), 1 )
 869 *
 870                END IF
 871 *
 872             ELSE
 873 *
 874 *              Complex left eigenvector.
 875 *
 876 *               Initial solve:
 877 *                 ((T(KI,KI)    T(KI,KI+1) )**T - (WR - I* WI))*X = 0.
 878 *                 ((T(KI+1,KI) T(KI+1,KI+1))                )
 879 *
 880                IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN
 881                   WORK( KI+N ) = WI / T( KI, KI+1 )
 882                   WORK( KI+1+N2 ) = ONE
 883                ELSE
 884                   WORK( KI+N ) = ONE
 885                   WORK( KI+1+N2 ) = -WI / T( KI+1, KI )
 886                END IF
 887                WORK( KI+1+N ) = ZERO
 888                WORK( KI+N2 ) = ZERO
 889 *
 890 *              Form right-hand side
 891 *
 892                DO 190 K = KI + 2, N
 893                   WORK( K+N ) = -WORK( KI+N )*T( KI, K )
 894                   WORK( K+N2 ) = -WORK( KI+1+N2 )*T( KI+1, K )
 895   190          CONTINUE
 896 *
 897 *              Solve complex quasi-triangular system:
 898 *              ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2
 899 *
 900                VMAX = ONE
 901                VCRIT = BIGNUM
 902 *
 903                JNXT = KI + 2
 904                DO 200 J = KI + 2, N
 905                   IF( J.LT.JNXT )
 906      $               GO TO 200
 907                   J1 = J
 908                   J2 = J
 909                   JNXT = J + 1
 910                   IF( J.LT.N ) THEN
 911                      IF( T( J+1, J ).NE.ZERO ) THEN
 912                         J2 = J + 1
 913                         JNXT = J + 2
 914                      END IF
 915                   END IF
 916 *
 917                   IF( J1.EQ.J2 ) THEN
 918 *
 919 *                    1-by-1 diagonal block
 920 *
 921 *                    Scale if necessary to avoid overflow when
 922 *                    forming the right-hand side elements.
 923 *
 924                      IF( WORK( J ).GT.VCRIT ) THEN
 925                         REC = ONE / VMAX
 926                         CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
 927                         CALL SSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
 928                         VMAX = ONE
 929                         VCRIT = BIGNUM
 930                      END IF
 931 *
 932                      WORK( J+N ) = WORK( J+N ) -
 933      $                             SDOT( J-KI-2, T( KI+2, J ), 1,
 934      $                             WORK( KI+2+N ), 1 )
 935                      WORK( J+N2 ) = WORK( J+N2 ) -
 936      $                              SDOT( J-KI-2, T( KI+2, J ), 1,
 937      $                              WORK( KI+2+N2 ), 1 )
 938 *
 939 *                    Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2
 940 *
 941                      CALL SLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
 942      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
 943      $                            -WI, X, 2, SCALE, XNORM, IERR )
 944 *
 945 *                    Scale if necessary
 946 *
 947                      IF( SCALE.NE.ONE ) THEN
 948                         CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
 949                         CALL SSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
 950                      END IF
 951                      WORK( J+N ) = X( 1, 1 )
 952                      WORK( J+N2 ) = X( 1, 2 )
 953                      VMAX = MAX( ABS( WORK( J+N ) ),
 954      $                      ABS( WORK( J+N2 ) ), VMAX )
 955                      VCRIT = BIGNUM / VMAX
 956 *
 957                   ELSE
 958 *
 959 *                    2-by-2 diagonal block
 960 *
 961 *                    Scale if necessary to avoid overflow when forming
 962 *                    the right-hand side elements.
 963 *
 964                      BETA = MAX( WORK( J ), WORK( J+1 ) )
 965                      IF( BETA.GT.VCRIT ) THEN
 966                         REC = ONE / VMAX
 967                         CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
 968                         CALL SSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
 969                         VMAX = ONE
 970                         VCRIT = BIGNUM
 971                      END IF
 972 *
 973                      WORK( J+N ) = WORK( J+N ) -
 974      $                             SDOT( J-KI-2, T( KI+2, J ), 1,
 975      $                             WORK( KI+2+N ), 1 )
 976 *
 977                      WORK( J+N2 ) = WORK( J+N2 ) -
 978      $                              SDOT( J-KI-2, T( KI+2, J ), 1,
 979      $                              WORK( KI+2+N2 ), 1 )
 980 *
 981                      WORK( J+1+N ) = WORK( J+1+N ) -
 982      $                               SDOT( J-KI-2, T( KI+2, J+1 ), 1,
 983      $                               WORK( KI+2+N ), 1 )
 984 *
 985                      WORK( J+1+N2 ) = WORK( J+1+N2 ) -
 986      $                                SDOT( J-KI-2, T( KI+2, J+1 ), 1,
 987      $                                WORK( KI+2+N2 ), 1 )
 988 *
 989 *                    Solve 2-by-2 complex linear equation
 990 *                      ([T(j,j)   T(j,j+1)  ]**T-(wr-i*wi)*I)*X = SCALE*B
 991 *                      ([T(j+1,j) T(j+1,j+1)]               )
 992 *
 993                      CALL SLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ),
 994      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
 995      $                            -WI, X, 2, SCALE, XNORM, IERR )
 996 *
 997 *                    Scale if necessary
 998 *
 999                      IF( SCALE.NE.ONE ) THEN
1000                         CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
1001                         CALL SSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
1002                      END IF
1003                      WORK( J+N ) = X( 1, 1 )
1004                      WORK( J+N2 ) = X( 1, 2 )
1005                      WORK( J+1+N ) = X( 2, 1 )
1006                      WORK( J+1+N2 ) = X( 2, 2 )
1007                      VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ),
1008      $                      ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ), VMAX )
1009                      VCRIT = BIGNUM / VMAX
1010 *
1011                   END IF
1012   200          CONTINUE
1013 *
1014 *              Copy the vector x or Q*x to VL and normalize.
1015 *
1016                IF( .NOT.OVER ) THEN
1017                   CALL SCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
1018                   CALL SCOPY( N-KI+1, WORK( KI+N2 ), 1, VL( KI, IS+1 ),
1019      $                        1 )
1020 *
1021                   EMAX = ZERO
1022                   DO 220 K = KI, N
1023                      EMAX = MAX( EMAX, ABS( VL( K, IS ) )+
1024      $                      ABS( VL( K, IS+1 ) ) )
1025   220             CONTINUE
1026                   REMAX = ONE / EMAX
1027                   CALL SSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
1028                   CALL SSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 )
1029 *
1030                   DO 230 K = 1, KI - 1
1031                      VL( K, IS ) = ZERO
1032                      VL( K, IS+1 ) = ZERO
1033   230             CONTINUE
1034                ELSE
1035                   IF( KI.LT.N-1 ) THEN
1036                      CALL SGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
1037      $                           LDVL, WORK( KI+2+N ), 1, WORK( KI+N ),
1038      $                           VL( 1, KI ), 1 )
1039                      CALL SGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
1040      $                           LDVL, WORK( KI+2+N2 ), 1,
1041      $                           WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
1042                   ELSE
1043                      CALL SSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 )
1044                      CALL SSCAL( N, WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
1045                   END IF
1046 *
1047                   EMAX = ZERO
1048                   DO 240 K = 1, N
1049                      EMAX = MAX( EMAX, ABS( VL( K, KI ) )+
1050      $                      ABS( VL( K, KI+1 ) ) )
1051   240             CONTINUE
1052                   REMAX = ONE / EMAX
1053                   CALL SSCAL( N, REMAX, VL( 1, KI ), 1 )
1054                   CALL SSCAL( N, REMAX, VL( 1, KI+1 ), 1 )
1055 *
1056                END IF
1057 *
1058             END IF
1059 *
1060             IS = IS + 1
1061             IF( IP.NE.0 )
1062      $         IS = IS + 1
1063   250       CONTINUE
1064             IF( IP.EQ.-1 )
1065      $         IP = 0
1066             IF( IP.EQ.1 )
1067      $         IP = -1
1068 *
1069   260    CONTINUE
1070 *
1071       END IF
1072 *
1073       RETURN
1074 *
1075 *     End of STREVC
1076 *
1077       END