LAPACK: SRC/strevc.f Source File

Go to the documentation of this file.
00001       SUBROUTINE STREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
00002      $                   LDVR, MM, M, WORK, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          HOWMNY, SIDE
00011       INTEGER            INFO, LDT, LDVL, LDVR, M, MM, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       LOGICAL            SELECT( * )
00015       REAL               T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
00016      $                   WORK( * )
00017 *     ..
00018 *
00019 *  Purpose
00020 *  =======
00021 *
00022 *  STREVC computes some or all of the right and/or left eigenvectors of
00023 *  a real upper quasi-triangular matrix T.
00024 *  Matrices of this type are produced by the Schur factorization of
00025 *  a real general matrix:  A = Q*T*Q**T, as computed by SHSEQR.
00026 *  
00027 *  The right eigenvector x and the left eigenvector y of T corresponding
00028 *  to an eigenvalue w are defined by:
00029 *  
00030 *     T*x = w*x,     (y**H)*T = w*(y**H)
00031 *  
00032 *  where y**H denotes the conjugate transpose of y.
00033 *  The eigenvalues are not input to this routine, but are read directly
00034 *  from the diagonal blocks of T.
00035 *  
00036 *  This routine returns the matrices X and/or Y of right and left
00037 *  eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
00038 *  input matrix.  If Q is the orthogonal factor that reduces a matrix
00039 *  A to Schur form T, then Q*X and Q*Y are the matrices of right and
00040 *  left eigenvectors of A.
00041 *
00042 *  Arguments
00043 *  =========
00044 *
00045 *  SIDE    (input) CHARACTER*1
00046 *          = 'R':  compute right eigenvectors only;
00047 *          = 'L':  compute left eigenvectors only;
00048 *          = 'B':  compute both right and left eigenvectors.
00049 *
00050 *  HOWMNY  (input) CHARACTER*1
00051 *          = 'A':  compute all right and/or left eigenvectors;
00052 *          = 'B':  compute all right and/or left eigenvectors,
00053 *                  backtransformed by the matrices in VR and/or VL;
00054 *          = 'S':  compute selected right and/or left eigenvectors,
00055 *                  as indicated by the logical array SELECT.
00056 *
00057 *  SELECT  (input/output) LOGICAL array, dimension (N)
00058 *          If HOWMNY = 'S', SELECT specifies the eigenvectors to be
00059 *          computed.
00060 *          If w(j) is a real eigenvalue, the corresponding real
00061 *          eigenvector is computed if SELECT(j) is .TRUE..
00062 *          If w(j) and w(j+1) are the real and imaginary parts of a
00063 *          complex eigenvalue, the corresponding complex eigenvector is
00064 *          computed if either SELECT(j) or SELECT(j+1) is .TRUE., and
00065 *          on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to
00066 *          .FALSE..
00067 *          Not referenced if HOWMNY = 'A' or 'B'.
00068 *
00069 *  N       (input) INTEGER
00070 *          The order of the matrix T. N >= 0.
00071 *
00072 *  T       (input) REAL array, dimension (LDT,N)
00073 *          The upper quasi-triangular matrix T in Schur canonical form.
00074 *
00075 *  LDT     (input) INTEGER
00076 *          The leading dimension of the array T. LDT >= max(1,N).
00077 *
00078 *  VL      (input/output) REAL array, dimension (LDVL,MM)
00079 *          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
00080 *          contain an N-by-N matrix Q (usually the orthogonal matrix Q
00081 *          of Schur vectors returned by SHSEQR).
00082 *          On exit, if SIDE = 'L' or 'B', VL contains:
00083 *          if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
00084 *          if HOWMNY = 'B', the matrix Q*Y;
00085 *          if HOWMNY = 'S', the left eigenvectors of T specified by
00086 *                           SELECT, stored consecutively in the columns
00087 *                           of VL, in the same order as their
00088 *                           eigenvalues.
00089 *          A complex eigenvector corresponding to a complex eigenvalue
00090 *          is stored in two consecutive columns, the first holding the
00091 *          real part, and the second the imaginary part.
00092 *          Not referenced if SIDE = 'R'.
00093 *
00094 *  LDVL    (input) INTEGER
00095 *          The leading dimension of the array VL.  LDVL >= 1, and if
00096 *          SIDE = 'L' or 'B', LDVL >= N.
00097 *
00098 *  VR      (input/output) REAL array, dimension (LDVR,MM)
00099 *          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
00100 *          contain an N-by-N matrix Q (usually the orthogonal matrix Q
00101 *          of Schur vectors returned by SHSEQR).
00102 *          On exit, if SIDE = 'R' or 'B', VR contains:
00103 *          if HOWMNY = 'A', the matrix X of right eigenvectors of T;
00104 *          if HOWMNY = 'B', the matrix Q*X;
00105 *          if HOWMNY = 'S', the right eigenvectors of T specified by
00106 *                           SELECT, stored consecutively in the columns
00107 *                           of VR, in the same order as their
00108 *                           eigenvalues.
00109 *          A complex eigenvector corresponding to a complex eigenvalue
00110 *          is stored in two consecutive columns, the first holding the
00111 *          real part and the second the imaginary part.
00112 *          Not referenced if SIDE = 'L'.
00113 *
00114 *  LDVR    (input) INTEGER
00115 *          The leading dimension of the array VR.  LDVR >= 1, and if
00116 *          SIDE = 'R' or 'B', LDVR >= N.
00117 *
00118 *  MM      (input) INTEGER
00119 *          The number of columns in the arrays VL and/or VR. MM >= M.
00120 *
00121 *  M       (output) INTEGER
00122 *          The number of columns in the arrays VL and/or VR actually
00123 *          used to store the eigenvectors.
00124 *          If HOWMNY = 'A' or 'B', M is set to N.
00125 *          Each selected real eigenvector occupies one column and each
00126 *          selected complex eigenvector occupies two columns.
00127 *
00128 *  WORK    (workspace) REAL array, dimension (3*N)
00129 *
00130 *  INFO    (output) INTEGER
00131 *          = 0:  successful exit
00132 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00133 *
00134 *  Further Details
00135 *  ===============
00136 *
00137 *  The algorithm used in this program is basically backward (forward)
00138 *  substitution, with scaling to make the the code robust against
00139 *  possible overflow.
00140 *
00141 *  Each eigenvector is normalized so that the element of largest
00142 *  magnitude has magnitude 1; here the magnitude of a complex number
00143 *  (x,y) is taken to be |x| + |y|.
00144 *
00145 *  =====================================================================
00146 *
00147 *     .. Parameters ..
00148       REAL               ZERO, ONE
00149       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00150 *     ..
00151 *     .. Local Scalars ..
00152       LOGICAL            ALLV, BOTHV, LEFTV, OVER, PAIR, RIGHTV, SOMEV
00153       INTEGER            I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI, N2
00154       REAL               BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
00155      $                   SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR,
00156      $                   XNORM
00157 *     ..
00158 *     .. External Functions ..
00159       LOGICAL            LSAME
00160       INTEGER            ISAMAX
00161       REAL               SDOT, SLAMCH
00162       EXTERNAL           LSAME, ISAMAX, SDOT, SLAMCH
00163 *     ..
00164 *     .. External Subroutines ..
00165       EXTERNAL           SAXPY, SCOPY, SGEMV, SLABAD, SLALN2, SSCAL,
00166      $                   XERBLA
00167 *     ..
00168 *     .. Intrinsic Functions ..
00169       INTRINSIC          ABS, MAX, SQRT
00170 *     ..
00171 *     .. Local Arrays ..
00172       REAL               X( 2, 2 )
00173 *     ..
00174 *     .. Executable Statements ..
00175 *
00176 *     Decode and test the input parameters
00177 *
00178       BOTHV = LSAME( SIDE, 'B' )
00179       RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
00180       LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
00181 *
00182       ALLV = LSAME( HOWMNY, 'A' )
00183       OVER = LSAME( HOWMNY, 'B' )
00184       SOMEV = LSAME( HOWMNY, 'S' )
00185 *
00186       INFO = 0
00187       IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
00188          INFO = -1
00189       ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
00190          INFO = -2
00191       ELSE IF( N.LT.0 ) THEN
00192          INFO = -4
00193       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
00194          INFO = -6
00195       ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
00196          INFO = -8
00197       ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
00198          INFO = -10
00199       ELSE
00200 *
00201 *        Set M to the number of columns required to store the selected
00202 *        eigenvectors, standardize the array SELECT if necessary, and
00203 *        test MM.
00204 *
00205          IF( SOMEV ) THEN
00206             M = 0
00207             PAIR = .FALSE.
00208             DO 10 J = 1, N
00209                IF( PAIR ) THEN
00210                   PAIR = .FALSE.
00211                   SELECT( J ) = .FALSE.
00212                ELSE
00213                   IF( J.LT.N ) THEN
00214                      IF( T( J+1, J ).EQ.ZERO ) THEN
00215                         IF( SELECT( J ) )
00216      $                     M = M + 1
00217                      ELSE
00218                         PAIR = .TRUE.
00219                         IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN
00220                            SELECT( J ) = .TRUE.
00221                            M = M + 2
00222                         END IF
00223                      END IF
00224                   ELSE
00225                      IF( SELECT( N ) )
00226      $                  M = M + 1
00227                   END IF
00228                END IF
00229    10       CONTINUE
00230          ELSE
00231             M = N
00232          END IF
00233 *
00234          IF( MM.LT.M ) THEN
00235             INFO = -11
00236          END IF
00237       END IF
00238       IF( INFO.NE.0 ) THEN
00239          CALL XERBLA( 'STREVC', -INFO )
00240          RETURN
00241       END IF
00242 *
00243 *     Quick return if possible.
00244 *
00245       IF( N.EQ.0 )
00246      $   RETURN
00247 *
00248 *     Set the constants to control overflow.
00249 *
00250       UNFL = SLAMCH( 'Safe minimum' )
00251       OVFL = ONE / UNFL
00252       CALL SLABAD( UNFL, OVFL )
00253       ULP = SLAMCH( 'Precision' )
00254       SMLNUM = UNFL*( N / ULP )
00255       BIGNUM = ( ONE-ULP ) / SMLNUM
00256 *
00257 *     Compute 1-norm of each column of strictly upper triangular
00258 *     part of T to control overflow in triangular solver.
00259 *
00260       WORK( 1 ) = ZERO
00261       DO 30 J = 2, N
00262          WORK( J ) = ZERO
00263          DO 20 I = 1, J - 1
00264             WORK( J ) = WORK( J ) + ABS( T( I, J ) )
00265    20    CONTINUE
00266    30 CONTINUE
00267 *
00268 *     Index IP is used to specify the real or complex eigenvalue:
00269 *       IP = 0, real eigenvalue,
00270 *            1, first of conjugate complex pair: (wr,wi)
00271 *           -1, second of conjugate complex pair: (wr,wi)
00272 *
00273       N2 = 2*N
00274 *
00275       IF( RIGHTV ) THEN
00276 *
00277 *        Compute right eigenvectors.
00278 *
00279          IP = 0
00280          IS = M
00281          DO 140 KI = N, 1, -1
00282 *
00283             IF( IP.EQ.1 )
00284      $         GO TO 130
00285             IF( KI.EQ.1 )
00286      $         GO TO 40
00287             IF( T( KI, KI-1 ).EQ.ZERO )
00288      $         GO TO 40
00289             IP = -1
00290 *
00291    40       CONTINUE
00292             IF( SOMEV ) THEN
00293                IF( IP.EQ.0 ) THEN
00294                   IF( .NOT.SELECT( KI ) )
00295      $               GO TO 130
00296                ELSE
00297                   IF( .NOT.SELECT( KI-1 ) )
00298      $               GO TO 130
00299                END IF
00300             END IF
00301 *
00302 *           Compute the KI-th eigenvalue (WR,WI).
00303 *
00304             WR = T( KI, KI )
00305             WI = ZERO
00306             IF( IP.NE.0 )
00307      $         WI = SQRT( ABS( T( KI, KI-1 ) ) )*
00308      $              SQRT( ABS( T( KI-1, KI ) ) )
00309             SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
00310 *
00311             IF( IP.EQ.0 ) THEN
00312 *
00313 *              Real right eigenvector
00314 *
00315                WORK( KI+N ) = ONE
00316 *
00317 *              Form right-hand side
00318 *
00319                DO 50 K = 1, KI - 1
00320                   WORK( K+N ) = -T( K, KI )
00321    50          CONTINUE
00322 *
00323 *              Solve the upper quasi-triangular system:
00324 *                 (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK.
00325 *
00326                JNXT = KI - 1
00327                DO 60 J = KI - 1, 1, -1
00328                   IF( J.GT.JNXT )
00329      $               GO TO 60
00330                   J1 = J
00331                   J2 = J
00332                   JNXT = J - 1
00333                   IF( J.GT.1 ) THEN
00334                      IF( T( J, J-1 ).NE.ZERO ) THEN
00335                         J1 = J - 1
00336                         JNXT = J - 2
00337                      END IF
00338                   END IF
00339 *
00340                   IF( J1.EQ.J2 ) THEN
00341 *
00342 *                    1-by-1 diagonal block
00343 *
00344                      CALL SLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
00345      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
00346      $                            ZERO, X, 2, SCALE, XNORM, IERR )
00347 *
00348 *                    Scale X(1,1) to avoid overflow when updating
00349 *                    the right-hand side.
00350 *
00351                      IF( XNORM.GT.ONE ) THEN
00352                         IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
00353                            X( 1, 1 ) = X( 1, 1 ) / XNORM
00354                            SCALE = SCALE / XNORM
00355                         END IF
00356                      END IF
00357 *
00358 *                    Scale if necessary
00359 *
00360                      IF( SCALE.NE.ONE )
00361      $                  CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 )
00362                      WORK( J+N ) = X( 1, 1 )
00363 *
00364 *                    Update right-hand side
00365 *
00366                      CALL SAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
00367      $                           WORK( 1+N ), 1 )
00368 *
00369                   ELSE
00370 *
00371 *                    2-by-2 diagonal block
00372 *
00373                      CALL SLALN2( .FALSE., 2, 1, SMIN, ONE,
00374      $                            T( J-1, J-1 ), LDT, ONE, ONE,
00375      $                            WORK( J-1+N ), N, WR, ZERO, X, 2,
00376      $                            SCALE, XNORM, IERR )
00377 *
00378 *                    Scale X(1,1) and X(2,1) to avoid overflow when
00379 *                    updating the right-hand side.
00380 *
00381                      IF( XNORM.GT.ONE ) THEN
00382                         BETA = MAX( WORK( J-1 ), WORK( J ) )
00383                         IF( BETA.GT.BIGNUM / XNORM ) THEN
00384                            X( 1, 1 ) = X( 1, 1 ) / XNORM
00385                            X( 2, 1 ) = X( 2, 1 ) / XNORM
00386                            SCALE = SCALE / XNORM
00387                         END IF
00388                      END IF
00389 *
00390 *                    Scale if necessary
00391 *
00392                      IF( SCALE.NE.ONE )
00393      $                  CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 )
00394                      WORK( J-1+N ) = X( 1, 1 )
00395                      WORK( J+N ) = X( 2, 1 )
00396 *
00397 *                    Update right-hand side
00398 *
00399                      CALL SAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
00400      $                           WORK( 1+N ), 1 )
00401                      CALL SAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
00402      $                           WORK( 1+N ), 1 )
00403                   END IF
00404    60          CONTINUE
00405 *
00406 *              Copy the vector x or Q*x to VR and normalize.
00407 *
00408                IF( .NOT.OVER ) THEN
00409                   CALL SCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 )
00410 *
00411                   II = ISAMAX( KI, VR( 1, IS ), 1 )
00412                   REMAX = ONE / ABS( VR( II, IS ) )
00413                   CALL SSCAL( KI, REMAX, VR( 1, IS ), 1 )
00414 *
00415                   DO 70 K = KI + 1, N
00416                      VR( K, IS ) = ZERO
00417    70             CONTINUE
00418                ELSE
00419                   IF( KI.GT.1 )
00420      $               CALL SGEMV( 'N', N, KI-1, ONE, VR, LDVR,
00421      $                           WORK( 1+N ), 1, WORK( KI+N ),
00422      $                           VR( 1, KI ), 1 )
00423 *
00424                   II = ISAMAX( N, VR( 1, KI ), 1 )
00425                   REMAX = ONE / ABS( VR( II, KI ) )
00426                   CALL SSCAL( N, REMAX, VR( 1, KI ), 1 )
00427                END IF
00428 *
00429             ELSE
00430 *
00431 *              Complex right eigenvector.
00432 *
00433 *              Initial solve
00434 *                [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.
00435 *                [ (T(KI,KI-1)   T(KI,KI)   )               ]
00436 *
00437                IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN
00438                   WORK( KI-1+N ) = ONE
00439                   WORK( KI+N2 ) = WI / T( KI-1, KI )
00440                ELSE
00441                   WORK( KI-1+N ) = -WI / T( KI, KI-1 )
00442                   WORK( KI+N2 ) = ONE
00443                END IF
00444                WORK( KI+N ) = ZERO
00445                WORK( KI-1+N2 ) = ZERO
00446 *
00447 *              Form right-hand side
00448 *
00449                DO 80 K = 1, KI - 2
00450                   WORK( K+N ) = -WORK( KI-1+N )*T( K, KI-1 )
00451                   WORK( K+N2 ) = -WORK( KI+N2 )*T( K, KI )
00452    80          CONTINUE
00453 *
00454 *              Solve upper quasi-triangular system:
00455 *              (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2)
00456 *
00457                JNXT = KI - 2
00458                DO 90 J = KI - 2, 1, -1
00459                   IF( J.GT.JNXT )
00460      $               GO TO 90
00461                   J1 = J
00462                   J2 = J
00463                   JNXT = J - 1
00464                   IF( J.GT.1 ) THEN
00465                      IF( T( J, J-1 ).NE.ZERO ) THEN
00466                         J1 = J - 1
00467                         JNXT = J - 2
00468                      END IF
00469                   END IF
00470 *
00471                   IF( J1.EQ.J2 ) THEN
00472 *
00473 *                    1-by-1 diagonal block
00474 *
00475                      CALL SLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
00476      $                            LDT, ONE, ONE, WORK( J+N ), N, WR, WI,
00477      $                            X, 2, SCALE, XNORM, IERR )
00478 *
00479 *                    Scale X(1,1) and X(1,2) to avoid overflow when
00480 *                    updating the right-hand side.
00481 *
00482                      IF( XNORM.GT.ONE ) THEN
00483                         IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
00484                            X( 1, 1 ) = X( 1, 1 ) / XNORM
00485                            X( 1, 2 ) = X( 1, 2 ) / XNORM
00486                            SCALE = SCALE / XNORM
00487                         END IF
00488                      END IF
00489 *
00490 *                    Scale if necessary
00491 *
00492                      IF( SCALE.NE.ONE ) THEN
00493                         CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 )
00494                         CALL SSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
00495                      END IF
00496                      WORK( J+N ) = X( 1, 1 )
00497                      WORK( J+N2 ) = X( 1, 2 )
00498 *
00499 *                    Update the right-hand side
00500 *
00501                      CALL SAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
00502      $                           WORK( 1+N ), 1 )
00503                      CALL SAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1,
00504      $                           WORK( 1+N2 ), 1 )
00505 *
00506                   ELSE
00507 *
00508 *                    2-by-2 diagonal block
00509 *
00510                      CALL SLALN2( .FALSE., 2, 2, SMIN, ONE,
00511      $                            T( J-1, J-1 ), LDT, ONE, ONE,
00512      $                            WORK( J-1+N ), N, WR, WI, X, 2, SCALE,
00513      $                            XNORM, IERR )
00514 *
00515 *                    Scale X to avoid overflow when updating
00516 *                    the right-hand side.
00517 *
00518                      IF( XNORM.GT.ONE ) THEN
00519                         BETA = MAX( WORK( J-1 ), WORK( J ) )
00520                         IF( BETA.GT.BIGNUM / XNORM ) THEN
00521                            REC = ONE / XNORM
00522                            X( 1, 1 ) = X( 1, 1 )*REC
00523                            X( 1, 2 ) = X( 1, 2 )*REC
00524                            X( 2, 1 ) = X( 2, 1 )*REC
00525                            X( 2, 2 ) = X( 2, 2 )*REC
00526                            SCALE = SCALE*REC
00527                         END IF
00528                      END IF
00529 *
00530 *                    Scale if necessary
00531 *
00532                      IF( SCALE.NE.ONE ) THEN
00533                         CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 )
00534                         CALL SSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
00535                      END IF
00536                      WORK( J-1+N ) = X( 1, 1 )
00537                      WORK( J+N ) = X( 2, 1 )
00538                      WORK( J-1+N2 ) = X( 1, 2 )
00539                      WORK( J+N2 ) = X( 2, 2 )
00540 *
00541 *                    Update the right-hand side
00542 *
00543                      CALL SAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
00544      $                           WORK( 1+N ), 1 )
00545                      CALL SAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
00546      $                           WORK( 1+N ), 1 )
00547                      CALL SAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
00548      $                           WORK( 1+N2 ), 1 )
00549                      CALL SAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1,
00550      $                           WORK( 1+N2 ), 1 )
00551                   END IF
00552    90          CONTINUE
00553 *
00554 *              Copy the vector x or Q*x to VR and normalize.
00555 *
00556                IF( .NOT.OVER ) THEN
00557                   CALL SCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 )
00558                   CALL SCOPY( KI, WORK( 1+N2 ), 1, VR( 1, IS ), 1 )
00559 *
00560                   EMAX = ZERO
00561                   DO 100 K = 1, KI
00562                      EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+
00563      $                      ABS( VR( K, IS ) ) )
00564   100             CONTINUE
00565 *
00566                   REMAX = ONE / EMAX
00567                   CALL SSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
00568                   CALL SSCAL( KI, REMAX, VR( 1, IS ), 1 )
00569 *
00570                   DO 110 K = KI + 1, N
00571                      VR( K, IS-1 ) = ZERO
00572                      VR( K, IS ) = ZERO
00573   110             CONTINUE
00574 *
00575                ELSE
00576 *
00577                   IF( KI.GT.2 ) THEN
00578                      CALL SGEMV( 'N', N, KI-2, ONE, VR, LDVR,
00579      $                           WORK( 1+N ), 1, WORK( KI-1+N ),
00580      $                           VR( 1, KI-1 ), 1 )
00581                      CALL SGEMV( 'N', N, KI-2, ONE, VR, LDVR,
00582      $                           WORK( 1+N2 ), 1, WORK( KI+N2 ),
00583      $                           VR( 1, KI ), 1 )
00584                   ELSE
00585                      CALL SSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 )
00586                      CALL SSCAL( N, WORK( KI+N2 ), VR( 1, KI ), 1 )
00587                   END IF
00588 *
00589                   EMAX = ZERO
00590                   DO 120 K = 1, N
00591                      EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+
00592      $                      ABS( VR( K, KI ) ) )
00593   120             CONTINUE
00594                   REMAX = ONE / EMAX
00595                   CALL SSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
00596                   CALL SSCAL( N, REMAX, VR( 1, KI ), 1 )
00597                END IF
00598             END IF
00599 *
00600             IS = IS - 1
00601             IF( IP.NE.0 )
00602      $         IS = IS - 1
00603   130       CONTINUE
00604             IF( IP.EQ.1 )
00605      $         IP = 0
00606             IF( IP.EQ.-1 )
00607      $         IP = 1
00608   140    CONTINUE
00609       END IF
00610 *
00611       IF( LEFTV ) THEN
00612 *
00613 *        Compute left eigenvectors.
00614 *
00615          IP = 0
00616          IS = 1
00617          DO 260 KI = 1, N
00618 *
00619             IF( IP.EQ.-1 )
00620      $         GO TO 250
00621             IF( KI.EQ.N )
00622      $         GO TO 150
00623             IF( T( KI+1, KI ).EQ.ZERO )
00624      $         GO TO 150
00625             IP = 1
00626 *
00627   150       CONTINUE
00628             IF( SOMEV ) THEN
00629                IF( .NOT.SELECT( KI ) )
00630      $            GO TO 250
00631             END IF
00632 *
00633 *           Compute the KI-th eigenvalue (WR,WI).
00634 *
00635             WR = T( KI, KI )
00636             WI = ZERO
00637             IF( IP.NE.0 )
00638      $         WI = SQRT( ABS( T( KI, KI+1 ) ) )*
00639      $              SQRT( ABS( T( KI+1, KI ) ) )
00640             SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
00641 *
00642             IF( IP.EQ.0 ) THEN
00643 *
00644 *              Real left eigenvector.
00645 *
00646                WORK( KI+N ) = ONE
00647 *
00648 *              Form right-hand side
00649 *
00650                DO 160 K = KI + 1, N
00651                   WORK( K+N ) = -T( KI, K )
00652   160          CONTINUE
00653 *
00654 *              Solve the quasi-triangular system:
00655 *                 (T(KI+1:N,KI+1:N) - WR)'*X = SCALE*WORK
00656 *
00657                VMAX = ONE
00658                VCRIT = BIGNUM
00659 *
00660                JNXT = KI + 1
00661                DO 170 J = KI + 1, N
00662                   IF( J.LT.JNXT )
00663      $               GO TO 170
00664                   J1 = J
00665                   J2 = J
00666                   JNXT = J + 1
00667                   IF( J.LT.N ) THEN
00668                      IF( T( J+1, J ).NE.ZERO ) THEN
00669                         J2 = J + 1
00670                         JNXT = J + 2
00671                      END IF
00672                   END IF
00673 *
00674                   IF( J1.EQ.J2 ) THEN
00675 *
00676 *                    1-by-1 diagonal block
00677 *
00678 *                    Scale if necessary to avoid overflow when forming
00679 *                    the right-hand side.
00680 *
00681                      IF( WORK( J ).GT.VCRIT ) THEN
00682                         REC = ONE / VMAX
00683                         CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
00684                         VMAX = ONE
00685                         VCRIT = BIGNUM
00686                      END IF
00687 *
00688                      WORK( J+N ) = WORK( J+N ) -
00689      $                             SDOT( J-KI-1, T( KI+1, J ), 1,
00690      $                             WORK( KI+1+N ), 1 )
00691 *
00692 *                    Solve (T(J,J)-WR)'*X = WORK
00693 *
00694                      CALL SLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
00695      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
00696      $                            ZERO, X, 2, SCALE, XNORM, IERR )
00697 *
00698 *                    Scale if necessary
00699 *
00700                      IF( SCALE.NE.ONE )
00701      $                  CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
00702                      WORK( J+N ) = X( 1, 1 )
00703                      VMAX = MAX( ABS( WORK( J+N ) ), VMAX )
00704                      VCRIT = BIGNUM / VMAX
00705 *
00706                   ELSE
00707 *
00708 *                    2-by-2 diagonal block
00709 *
00710 *                    Scale if necessary to avoid overflow when forming
00711 *                    the right-hand side.
00712 *
00713                      BETA = MAX( WORK( J ), WORK( J+1 ) )
00714                      IF( BETA.GT.VCRIT ) THEN
00715                         REC = ONE / VMAX
00716                         CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
00717                         VMAX = ONE
00718                         VCRIT = BIGNUM
00719                      END IF
00720 *
00721                      WORK( J+N ) = WORK( J+N ) -
00722      $                             SDOT( J-KI-1, T( KI+1, J ), 1,
00723      $                             WORK( KI+1+N ), 1 )
00724 *
00725                      WORK( J+1+N ) = WORK( J+1+N ) -
00726      $                               SDOT( J-KI-1, T( KI+1, J+1 ), 1,
00727      $                               WORK( KI+1+N ), 1 )
00728 *
00729 *                    Solve
00730 *                      [T(J,J)-WR   T(J,J+1)     ]'* X = SCALE*( WORK1 )
00731 *                      [T(J+1,J)    T(J+1,J+1)-WR]             ( WORK2 )
00732 *
00733                      CALL SLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ),
00734      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
00735      $                            ZERO, X, 2, SCALE, XNORM, IERR )
00736 *
00737 *                    Scale if necessary
00738 *
00739                      IF( SCALE.NE.ONE )
00740      $                  CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
00741                      WORK( J+N ) = X( 1, 1 )
00742                      WORK( J+1+N ) = X( 2, 1 )
00743 *
00744                      VMAX = MAX( ABS( WORK( J+N ) ),
00745      $                      ABS( WORK( J+1+N ) ), VMAX )
00746                      VCRIT = BIGNUM / VMAX
00747 *
00748                   END IF
00749   170          CONTINUE
00750 *
00751 *              Copy the vector x or Q*x to VL and normalize.
00752 *
00753                IF( .NOT.OVER ) THEN
00754                   CALL SCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
00755 *
00756                   II = ISAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
00757                   REMAX = ONE / ABS( VL( II, IS ) )
00758                   CALL SSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
00759 *
00760                   DO 180 K = 1, KI - 1
00761                      VL( K, IS ) = ZERO
00762   180             CONTINUE
00763 *
00764                ELSE
00765 *
00766                   IF( KI.LT.N )
00767      $               CALL SGEMV( 'N', N, N-KI, ONE, VL( 1, KI+1 ), LDVL,
00768      $                           WORK( KI+1+N ), 1, WORK( KI+N ),
00769      $                           VL( 1, KI ), 1 )
00770 *
00771                   II = ISAMAX( N, VL( 1, KI ), 1 )
00772                   REMAX = ONE / ABS( VL( II, KI ) )
00773                   CALL SSCAL( N, REMAX, VL( 1, KI ), 1 )
00774 *
00775                END IF
00776 *
00777             ELSE
00778 *
00779 *              Complex left eigenvector.
00780 *
00781 *               Initial solve:
00782 *                 ((T(KI,KI)    T(KI,KI+1) )' - (WR - I* WI))*X = 0.
00783 *                 ((T(KI+1,KI) T(KI+1,KI+1))                )
00784 *
00785                IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN
00786                   WORK( KI+N ) = WI / T( KI, KI+1 )
00787                   WORK( KI+1+N2 ) = ONE
00788                ELSE
00789                   WORK( KI+N ) = ONE
00790                   WORK( KI+1+N2 ) = -WI / T( KI+1, KI )
00791                END IF
00792                WORK( KI+1+N ) = ZERO
00793                WORK( KI+N2 ) = ZERO
00794 *
00795 *              Form right-hand side
00796 *
00797                DO 190 K = KI + 2, N
00798                   WORK( K+N ) = -WORK( KI+N )*T( KI, K )
00799                   WORK( K+N2 ) = -WORK( KI+1+N2 )*T( KI+1, K )
00800   190          CONTINUE
00801 *
00802 *              Solve complex quasi-triangular system:
00803 *              ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2
00804 *
00805                VMAX = ONE
00806                VCRIT = BIGNUM
00807 *
00808                JNXT = KI + 2
00809                DO 200 J = KI + 2, N
00810                   IF( J.LT.JNXT )
00811      $               GO TO 200
00812                   J1 = J
00813                   J2 = J
00814                   JNXT = J + 1
00815                   IF( J.LT.N ) THEN
00816                      IF( T( J+1, J ).NE.ZERO ) THEN
00817                         J2 = J + 1
00818                         JNXT = J + 2
00819                      END IF
00820                   END IF
00821 *
00822                   IF( J1.EQ.J2 ) THEN
00823 *
00824 *                    1-by-1 diagonal block
00825 *
00826 *                    Scale if necessary to avoid overflow when
00827 *                    forming the right-hand side elements.
00828 *
00829                      IF( WORK( J ).GT.VCRIT ) THEN
00830                         REC = ONE / VMAX
00831                         CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
00832                         CALL SSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
00833                         VMAX = ONE
00834                         VCRIT = BIGNUM
00835                      END IF
00836 *
00837                      WORK( J+N ) = WORK( J+N ) -
00838      $                             SDOT( J-KI-2, T( KI+2, J ), 1,
00839      $                             WORK( KI+2+N ), 1 )
00840                      WORK( J+N2 ) = WORK( J+N2 ) -
00841      $                              SDOT( J-KI-2, T( KI+2, J ), 1,
00842      $                              WORK( KI+2+N2 ), 1 )
00843 *
00844 *                    Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2
00845 *
00846                      CALL SLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
00847      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
00848      $                            -WI, X, 2, SCALE, XNORM, IERR )
00849 *
00850 *                    Scale if necessary
00851 *
00852                      IF( SCALE.NE.ONE ) THEN
00853                         CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
00854                         CALL SSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
00855                      END IF
00856                      WORK( J+N ) = X( 1, 1 )
00857                      WORK( J+N2 ) = X( 1, 2 )
00858                      VMAX = MAX( ABS( WORK( J+N ) ),
00859      $                      ABS( WORK( J+N2 ) ), VMAX )
00860                      VCRIT = BIGNUM / VMAX
00861 *
00862                   ELSE
00863 *
00864 *                    2-by-2 diagonal block
00865 *
00866 *                    Scale if necessary to avoid overflow when forming
00867 *                    the right-hand side elements.
00868 *
00869                      BETA = MAX( WORK( J ), WORK( J+1 ) )
00870                      IF( BETA.GT.VCRIT ) THEN
00871                         REC = ONE / VMAX
00872                         CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
00873                         CALL SSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
00874                         VMAX = ONE
00875                         VCRIT = BIGNUM
00876                      END IF
00877 *
00878                      WORK( J+N ) = WORK( J+N ) -
00879      $                             SDOT( J-KI-2, T( KI+2, J ), 1,
00880      $                             WORK( KI+2+N ), 1 )
00881 *
00882                      WORK( J+N2 ) = WORK( J+N2 ) -
00883      $                              SDOT( J-KI-2, T( KI+2, J ), 1,
00884      $                              WORK( KI+2+N2 ), 1 )
00885 *
00886                      WORK( J+1+N ) = WORK( J+1+N ) -
00887      $                               SDOT( J-KI-2, T( KI+2, J+1 ), 1,
00888      $                               WORK( KI+2+N ), 1 )
00889 *
00890                      WORK( J+1+N2 ) = WORK( J+1+N2 ) -
00891      $                                SDOT( J-KI-2, T( KI+2, J+1 ), 1,
00892      $                                WORK( KI+2+N2 ), 1 )
00893 *
00894 *                    Solve 2-by-2 complex linear equation
00895 *                      ([T(j,j)   T(j,j+1)  ]'-(wr-i*wi)*I)*X = SCALE*B
00896 *                      ([T(j+1,j) T(j+1,j+1)]             )
00897 *
00898                      CALL SLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ),
00899      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
00900      $                            -WI, X, 2, SCALE, XNORM, IERR )
00901 *
00902 *                    Scale if necessary
00903 *
00904                      IF( SCALE.NE.ONE ) THEN
00905                         CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
00906                         CALL SSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
00907                      END IF
00908                      WORK( J+N ) = X( 1, 1 )
00909                      WORK( J+N2 ) = X( 1, 2 )
00910                      WORK( J+1+N ) = X( 2, 1 )
00911                      WORK( J+1+N2 ) = X( 2, 2 )
00912                      VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ),
00913      $                      ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ), VMAX )
00914                      VCRIT = BIGNUM / VMAX
00915 *
00916                   END IF
00917   200          CONTINUE
00918 *
00919 *              Copy the vector x or Q*x to VL and normalize.
00920 *
00921                IF( .NOT.OVER ) THEN
00922                   CALL SCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
00923                   CALL SCOPY( N-KI+1, WORK( KI+N2 ), 1, VL( KI, IS+1 ),
00924      $                        1 )
00925 *
00926                   EMAX = ZERO
00927                   DO 220 K = KI, N
00928                      EMAX = MAX( EMAX, ABS( VL( K, IS ) )+
00929      $                      ABS( VL( K, IS+1 ) ) )
00930   220             CONTINUE
00931                   REMAX = ONE / EMAX
00932                   CALL SSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
00933                   CALL SSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 )
00934 *
00935                   DO 230 K = 1, KI - 1
00936                      VL( K, IS ) = ZERO
00937                      VL( K, IS+1 ) = ZERO
00938   230             CONTINUE
00939                ELSE
00940                   IF( KI.LT.N-1 ) THEN
00941                      CALL SGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
00942      $                           LDVL, WORK( KI+2+N ), 1, WORK( KI+N ),
00943      $                           VL( 1, KI ), 1 )
00944                      CALL SGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
00945      $                           LDVL, WORK( KI+2+N2 ), 1,
00946      $                           WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
00947                   ELSE
00948                      CALL SSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 )
00949                      CALL SSCAL( N, WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
00950                   END IF
00951 *
00952                   EMAX = ZERO
00953                   DO 240 K = 1, N
00954                      EMAX = MAX( EMAX, ABS( VL( K, KI ) )+
00955      $                      ABS( VL( K, KI+1 ) ) )
00956   240             CONTINUE
00957                   REMAX = ONE / EMAX
00958                   CALL SSCAL( N, REMAX, VL( 1, KI ), 1 )
00959                   CALL SSCAL( N, REMAX, VL( 1, KI+1 ), 1 )
00960 *
00961                END IF
00962 *
00963             END IF
00964 *
00965             IS = IS + 1
00966             IF( IP.NE.0 )
00967      $         IS = IS + 1
00968   250       CONTINUE
00969             IF( IP.EQ.-1 )
00970      $         IP = 0
00971             IF( IP.EQ.1 )
00972      $         IP = -1
00973 *
00974   260    CONTINUE
00975 *
00976       END IF
00977 *
00978       RETURN
00979 *
00980 *     End of STREVC
00981 *
00982       END