LAPACK: SRC/dtrevc.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DTREVC( 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       DOUBLE PRECISION   T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
00016      $                   WORK( * )
00017 *     ..
00018 *
00019 *  Purpose
00020 *  =======
00021 *
00022 *  DTREVC 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 DHSEQR.
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) DOUBLE PRECISION 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) DOUBLE PRECISION 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 DHSEQR).
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) DOUBLE PRECISION 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 DHSEQR).
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) DOUBLE PRECISION 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       DOUBLE PRECISION   ZERO, ONE
00149       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+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       DOUBLE PRECISION   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            IDAMAX
00161       DOUBLE PRECISION   DDOT, DLAMCH
00162       EXTERNAL           LSAME, IDAMAX, DDOT, DLAMCH
00163 *     ..
00164 *     .. External Subroutines ..
00165       EXTERNAL           DAXPY, DCOPY, DGEMV, DLALN2, DSCAL, XERBLA
00166 *     ..
00167 *     .. Intrinsic Functions ..
00168       INTRINSIC          ABS, MAX, SQRT
00169 *     ..
00170 *     .. Local Arrays ..
00171       DOUBLE PRECISION   X( 2, 2 )
00172 *     ..
00173 *     .. Executable Statements ..
00174 *
00175 *     Decode and test the input parameters
00176 *
00177       BOTHV = LSAME( SIDE, 'B' )
00178       RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
00179       LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
00180 *
00181       ALLV = LSAME( HOWMNY, 'A' )
00182       OVER = LSAME( HOWMNY, 'B' )
00183       SOMEV = LSAME( HOWMNY, 'S' )
00184 *
00185       INFO = 0
00186       IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
00187          INFO = -1
00188       ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
00189          INFO = -2
00190       ELSE IF( N.LT.0 ) THEN
00191          INFO = -4
00192       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
00193          INFO = -6
00194       ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
00195          INFO = -8
00196       ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
00197          INFO = -10
00198       ELSE
00199 *
00200 *        Set M to the number of columns required to store the selected
00201 *        eigenvectors, standardize the array SELECT if necessary, and
00202 *        test MM.
00203 *
00204          IF( SOMEV ) THEN
00205             M = 0
00206             PAIR = .FALSE.
00207             DO 10 J = 1, N
00208                IF( PAIR ) THEN
00209                   PAIR = .FALSE.
00210                   SELECT( J ) = .FALSE.
00211                ELSE
00212                   IF( J.LT.N ) THEN
00213                      IF( T( J+1, J ).EQ.ZERO ) THEN
00214                         IF( SELECT( J ) )
00215      $                     M = M + 1
00216                      ELSE
00217                         PAIR = .TRUE.
00218                         IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN
00219                            SELECT( J ) = .TRUE.
00220                            M = M + 2
00221                         END IF
00222                      END IF
00223                   ELSE
00224                      IF( SELECT( N ) )
00225      $                  M = M + 1
00226                   END IF
00227                END IF
00228    10       CONTINUE
00229          ELSE
00230             M = N
00231          END IF
00232 *
00233          IF( MM.LT.M ) THEN
00234             INFO = -11
00235          END IF
00236       END IF
00237       IF( INFO.NE.0 ) THEN
00238          CALL XERBLA( 'DTREVC', -INFO )
00239          RETURN
00240       END IF
00241 *
00242 *     Quick return if possible.
00243 *
00244       IF( N.EQ.0 )
00245      $   RETURN
00246 *
00247 *     Set the constants to control overflow.
00248 *
00249       UNFL = DLAMCH( 'Safe minimum' )
00250       OVFL = ONE / UNFL
00251       CALL DLABAD( UNFL, OVFL )
00252       ULP = DLAMCH( 'Precision' )
00253       SMLNUM = UNFL*( N / ULP )
00254       BIGNUM = ( ONE-ULP ) / SMLNUM
00255 *
00256 *     Compute 1-norm of each column of strictly upper triangular
00257 *     part of T to control overflow in triangular solver.
00258 *
00259       WORK( 1 ) = ZERO
00260       DO 30 J = 2, N
00261          WORK( J ) = ZERO
00262          DO 20 I = 1, J - 1
00263             WORK( J ) = WORK( J ) + ABS( T( I, J ) )
00264    20    CONTINUE
00265    30 CONTINUE
00266 *
00267 *     Index IP is used to specify the real or complex eigenvalue:
00268 *       IP = 0, real eigenvalue,
00269 *            1, first of conjugate complex pair: (wr,wi)
00270 *           -1, second of conjugate complex pair: (wr,wi)
00271 *
00272       N2 = 2*N
00273 *
00274       IF( RIGHTV ) THEN
00275 *
00276 *        Compute right eigenvectors.
00277 *
00278          IP = 0
00279          IS = M
00280          DO 140 KI = N, 1, -1
00281 *
00282             IF( IP.EQ.1 )
00283      $         GO TO 130
00284             IF( KI.EQ.1 )
00285      $         GO TO 40
00286             IF( T( KI, KI-1 ).EQ.ZERO )
00287      $         GO TO 40
00288             IP = -1
00289 *
00290    40       CONTINUE
00291             IF( SOMEV ) THEN
00292                IF( IP.EQ.0 ) THEN
00293                   IF( .NOT.SELECT( KI ) )
00294      $               GO TO 130
00295                ELSE
00296                   IF( .NOT.SELECT( KI-1 ) )
00297      $               GO TO 130
00298                END IF
00299             END IF
00300 *
00301 *           Compute the KI-th eigenvalue (WR,WI).
00302 *
00303             WR = T( KI, KI )
00304             WI = ZERO
00305             IF( IP.NE.0 )
00306      $         WI = SQRT( ABS( T( KI, KI-1 ) ) )*
00307      $              SQRT( ABS( T( KI-1, KI ) ) )
00308             SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
00309 *
00310             IF( IP.EQ.0 ) THEN
00311 *
00312 *              Real right eigenvector
00313 *
00314                WORK( KI+N ) = ONE
00315 *
00316 *              Form right-hand side
00317 *
00318                DO 50 K = 1, KI - 1
00319                   WORK( K+N ) = -T( K, KI )
00320    50          CONTINUE
00321 *
00322 *              Solve the upper quasi-triangular system:
00323 *                 (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK.
00324 *
00325                JNXT = KI - 1
00326                DO 60 J = KI - 1, 1, -1
00327                   IF( J.GT.JNXT )
00328      $               GO TO 60
00329                   J1 = J
00330                   J2 = J
00331                   JNXT = J - 1
00332                   IF( J.GT.1 ) THEN
00333                      IF( T( J, J-1 ).NE.ZERO ) THEN
00334                         J1 = J - 1
00335                         JNXT = J - 2
00336                      END IF
00337                   END IF
00338 *
00339                   IF( J1.EQ.J2 ) THEN
00340 *
00341 *                    1-by-1 diagonal block
00342 *
00343                      CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
00344      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
00345      $                            ZERO, X, 2, SCALE, XNORM, IERR )
00346 *
00347 *                    Scale X(1,1) to avoid overflow when updating
00348 *                    the right-hand side.
00349 *
00350                      IF( XNORM.GT.ONE ) THEN
00351                         IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
00352                            X( 1, 1 ) = X( 1, 1 ) / XNORM
00353                            SCALE = SCALE / XNORM
00354                         END IF
00355                      END IF
00356 *
00357 *                    Scale if necessary
00358 *
00359                      IF( SCALE.NE.ONE )
00360      $                  CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
00361                      WORK( J+N ) = X( 1, 1 )
00362 *
00363 *                    Update right-hand side
00364 *
00365                      CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
00366      $                           WORK( 1+N ), 1 )
00367 *
00368                   ELSE
00369 *
00370 *                    2-by-2 diagonal block
00371 *
00372                      CALL DLALN2( .FALSE., 2, 1, SMIN, ONE,
00373      $                            T( J-1, J-1 ), LDT, ONE, ONE,
00374      $                            WORK( J-1+N ), N, WR, ZERO, X, 2,
00375      $                            SCALE, XNORM, IERR )
00376 *
00377 *                    Scale X(1,1) and X(2,1) to avoid overflow when
00378 *                    updating the right-hand side.
00379 *
00380                      IF( XNORM.GT.ONE ) THEN
00381                         BETA = MAX( WORK( J-1 ), WORK( J ) )
00382                         IF( BETA.GT.BIGNUM / XNORM ) THEN
00383                            X( 1, 1 ) = X( 1, 1 ) / XNORM
00384                            X( 2, 1 ) = X( 2, 1 ) / XNORM
00385                            SCALE = SCALE / XNORM
00386                         END IF
00387                      END IF
00388 *
00389 *                    Scale if necessary
00390 *
00391                      IF( SCALE.NE.ONE )
00392      $                  CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
00393                      WORK( J-1+N ) = X( 1, 1 )
00394                      WORK( J+N ) = X( 2, 1 )
00395 *
00396 *                    Update right-hand side
00397 *
00398                      CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
00399      $                           WORK( 1+N ), 1 )
00400                      CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
00401      $                           WORK( 1+N ), 1 )
00402                   END IF
00403    60          CONTINUE
00404 *
00405 *              Copy the vector x or Q*x to VR and normalize.
00406 *
00407                IF( .NOT.OVER ) THEN
00408                   CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 )
00409 *
00410                   II = IDAMAX( KI, VR( 1, IS ), 1 )
00411                   REMAX = ONE / ABS( VR( II, IS ) )
00412                   CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
00413 *
00414                   DO 70 K = KI + 1, N
00415                      VR( K, IS ) = ZERO
00416    70             CONTINUE
00417                ELSE
00418                   IF( KI.GT.1 )
00419      $               CALL DGEMV( 'N', N, KI-1, ONE, VR, LDVR,
00420      $                           WORK( 1+N ), 1, WORK( KI+N ),
00421      $                           VR( 1, KI ), 1 )
00422 *
00423                   II = IDAMAX( N, VR( 1, KI ), 1 )
00424                   REMAX = ONE / ABS( VR( II, KI ) )
00425                   CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
00426                END IF
00427 *
00428             ELSE
00429 *
00430 *              Complex right eigenvector.
00431 *
00432 *              Initial solve
00433 *                [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.
00434 *                [ (T(KI,KI-1)   T(KI,KI)   )               ]
00435 *
00436                IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN
00437                   WORK( KI-1+N ) = ONE
00438                   WORK( KI+N2 ) = WI / T( KI-1, KI )
00439                ELSE
00440                   WORK( KI-1+N ) = -WI / T( KI, KI-1 )
00441                   WORK( KI+N2 ) = ONE
00442                END IF
00443                WORK( KI+N ) = ZERO
00444                WORK( KI-1+N2 ) = ZERO
00445 *
00446 *              Form right-hand side
00447 *
00448                DO 80 K = 1, KI - 2
00449                   WORK( K+N ) = -WORK( KI-1+N )*T( K, KI-1 )
00450                   WORK( K+N2 ) = -WORK( KI+N2 )*T( K, KI )
00451    80          CONTINUE
00452 *
00453 *              Solve upper quasi-triangular system:
00454 *              (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2)
00455 *
00456                JNXT = KI - 2
00457                DO 90 J = KI - 2, 1, -1
00458                   IF( J.GT.JNXT )
00459      $               GO TO 90
00460                   J1 = J
00461                   J2 = J
00462                   JNXT = J - 1
00463                   IF( J.GT.1 ) THEN
00464                      IF( T( J, J-1 ).NE.ZERO ) THEN
00465                         J1 = J - 1
00466                         JNXT = J - 2
00467                      END IF
00468                   END IF
00469 *
00470                   IF( J1.EQ.J2 ) THEN
00471 *
00472 *                    1-by-1 diagonal block
00473 *
00474                      CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
00475      $                            LDT, ONE, ONE, WORK( J+N ), N, WR, WI,
00476      $                            X, 2, SCALE, XNORM, IERR )
00477 *
00478 *                    Scale X(1,1) and X(1,2) to avoid overflow when
00479 *                    updating the right-hand side.
00480 *
00481                      IF( XNORM.GT.ONE ) THEN
00482                         IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
00483                            X( 1, 1 ) = X( 1, 1 ) / XNORM
00484                            X( 1, 2 ) = X( 1, 2 ) / XNORM
00485                            SCALE = SCALE / XNORM
00486                         END IF
00487                      END IF
00488 *
00489 *                    Scale if necessary
00490 *
00491                      IF( SCALE.NE.ONE ) THEN
00492                         CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
00493                         CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
00494                      END IF
00495                      WORK( J+N ) = X( 1, 1 )
00496                      WORK( J+N2 ) = X( 1, 2 )
00497 *
00498 *                    Update the right-hand side
00499 *
00500                      CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
00501      $                           WORK( 1+N ), 1 )
00502                      CALL DAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1,
00503      $                           WORK( 1+N2 ), 1 )
00504 *
00505                   ELSE
00506 *
00507 *                    2-by-2 diagonal block
00508 *
00509                      CALL DLALN2( .FALSE., 2, 2, SMIN, ONE,
00510      $                            T( J-1, J-1 ), LDT, ONE, ONE,
00511      $                            WORK( J-1+N ), N, WR, WI, X, 2, SCALE,
00512      $                            XNORM, IERR )
00513 *
00514 *                    Scale X to avoid overflow when updating
00515 *                    the right-hand side.
00516 *
00517                      IF( XNORM.GT.ONE ) THEN
00518                         BETA = MAX( WORK( J-1 ), WORK( J ) )
00519                         IF( BETA.GT.BIGNUM / XNORM ) THEN
00520                            REC = ONE / XNORM
00521                            X( 1, 1 ) = X( 1, 1 )*REC
00522                            X( 1, 2 ) = X( 1, 2 )*REC
00523                            X( 2, 1 ) = X( 2, 1 )*REC
00524                            X( 2, 2 ) = X( 2, 2 )*REC
00525                            SCALE = SCALE*REC
00526                         END IF
00527                      END IF
00528 *
00529 *                    Scale if necessary
00530 *
00531                      IF( SCALE.NE.ONE ) THEN
00532                         CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
00533                         CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
00534                      END IF
00535                      WORK( J-1+N ) = X( 1, 1 )
00536                      WORK( J+N ) = X( 2, 1 )
00537                      WORK( J-1+N2 ) = X( 1, 2 )
00538                      WORK( J+N2 ) = X( 2, 2 )
00539 *
00540 *                    Update the right-hand side
00541 *
00542                      CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
00543      $                           WORK( 1+N ), 1 )
00544                      CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
00545      $                           WORK( 1+N ), 1 )
00546                      CALL DAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
00547      $                           WORK( 1+N2 ), 1 )
00548                      CALL DAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1,
00549      $                           WORK( 1+N2 ), 1 )
00550                   END IF
00551    90          CONTINUE
00552 *
00553 *              Copy the vector x or Q*x to VR and normalize.
00554 *
00555                IF( .NOT.OVER ) THEN
00556                   CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 )
00557                   CALL DCOPY( KI, WORK( 1+N2 ), 1, VR( 1, IS ), 1 )
00558 *
00559                   EMAX = ZERO
00560                   DO 100 K = 1, KI
00561                      EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+
00562      $                      ABS( VR( K, IS ) ) )
00563   100             CONTINUE
00564 *
00565                   REMAX = ONE / EMAX
00566                   CALL DSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
00567                   CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
00568 *
00569                   DO 110 K = KI + 1, N
00570                      VR( K, IS-1 ) = ZERO
00571                      VR( K, IS ) = ZERO
00572   110             CONTINUE
00573 *
00574                ELSE
00575 *
00576                   IF( KI.GT.2 ) THEN
00577                      CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
00578      $                           WORK( 1+N ), 1, WORK( KI-1+N ),
00579      $                           VR( 1, KI-1 ), 1 )
00580                      CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
00581      $                           WORK( 1+N2 ), 1, WORK( KI+N2 ),
00582      $                           VR( 1, KI ), 1 )
00583                   ELSE
00584                      CALL DSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 )
00585                      CALL DSCAL( N, WORK( KI+N2 ), VR( 1, KI ), 1 )
00586                   END IF
00587 *
00588                   EMAX = ZERO
00589                   DO 120 K = 1, N
00590                      EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+
00591      $                      ABS( VR( K, KI ) ) )
00592   120             CONTINUE
00593                   REMAX = ONE / EMAX
00594                   CALL DSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
00595                   CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
00596                END IF
00597             END IF
00598 *
00599             IS = IS - 1
00600             IF( IP.NE.0 )
00601      $         IS = IS - 1
00602   130       CONTINUE
00603             IF( IP.EQ.1 )
00604      $         IP = 0
00605             IF( IP.EQ.-1 )
00606      $         IP = 1
00607   140    CONTINUE
00608       END IF
00609 *
00610       IF( LEFTV ) THEN
00611 *
00612 *        Compute left eigenvectors.
00613 *
00614          IP = 0
00615          IS = 1
00616          DO 260 KI = 1, N
00617 *
00618             IF( IP.EQ.-1 )
00619      $         GO TO 250
00620             IF( KI.EQ.N )
00621      $         GO TO 150
00622             IF( T( KI+1, KI ).EQ.ZERO )
00623      $         GO TO 150
00624             IP = 1
00625 *
00626   150       CONTINUE
00627             IF( SOMEV ) THEN
00628                IF( .NOT.SELECT( KI ) )
00629      $            GO TO 250
00630             END IF
00631 *
00632 *           Compute the KI-th eigenvalue (WR,WI).
00633 *
00634             WR = T( KI, KI )
00635             WI = ZERO
00636             IF( IP.NE.0 )
00637      $         WI = SQRT( ABS( T( KI, KI+1 ) ) )*
00638      $              SQRT( ABS( T( KI+1, KI ) ) )
00639             SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
00640 *
00641             IF( IP.EQ.0 ) THEN
00642 *
00643 *              Real left eigenvector.
00644 *
00645                WORK( KI+N ) = ONE
00646 *
00647 *              Form right-hand side
00648 *
00649                DO 160 K = KI + 1, N
00650                   WORK( K+N ) = -T( KI, K )
00651   160          CONTINUE
00652 *
00653 *              Solve the quasi-triangular system:
00654 *                 (T(KI+1:N,KI+1:N) - WR)'*X = SCALE*WORK
00655 *
00656                VMAX = ONE
00657                VCRIT = BIGNUM
00658 *
00659                JNXT = KI + 1
00660                DO 170 J = KI + 1, N
00661                   IF( J.LT.JNXT )
00662      $               GO TO 170
00663                   J1 = J
00664                   J2 = J
00665                   JNXT = J + 1
00666                   IF( J.LT.N ) THEN
00667                      IF( T( J+1, J ).NE.ZERO ) THEN
00668                         J2 = J + 1
00669                         JNXT = J + 2
00670                      END IF
00671                   END IF
00672 *
00673                   IF( J1.EQ.J2 ) THEN
00674 *
00675 *                    1-by-1 diagonal block
00676 *
00677 *                    Scale if necessary to avoid overflow when forming
00678 *                    the right-hand side.
00679 *
00680                      IF( WORK( J ).GT.VCRIT ) THEN
00681                         REC = ONE / VMAX
00682                         CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
00683                         VMAX = ONE
00684                         VCRIT = BIGNUM
00685                      END IF
00686 *
00687                      WORK( J+N ) = WORK( J+N ) -
00688      $                             DDOT( J-KI-1, T( KI+1, J ), 1,
00689      $                             WORK( KI+1+N ), 1 )
00690 *
00691 *                    Solve (T(J,J)-WR)'*X = WORK
00692 *
00693                      CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
00694      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
00695      $                            ZERO, X, 2, SCALE, XNORM, IERR )
00696 *
00697 *                    Scale if necessary
00698 *
00699                      IF( SCALE.NE.ONE )
00700      $                  CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
00701                      WORK( J+N ) = X( 1, 1 )
00702                      VMAX = MAX( ABS( WORK( J+N ) ), VMAX )
00703                      VCRIT = BIGNUM / VMAX
00704 *
00705                   ELSE
00706 *
00707 *                    2-by-2 diagonal block
00708 *
00709 *                    Scale if necessary to avoid overflow when forming
00710 *                    the right-hand side.
00711 *
00712                      BETA = MAX( WORK( J ), WORK( J+1 ) )
00713                      IF( BETA.GT.VCRIT ) THEN
00714                         REC = ONE / VMAX
00715                         CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
00716                         VMAX = ONE
00717                         VCRIT = BIGNUM
00718                      END IF
00719 *
00720                      WORK( J+N ) = WORK( J+N ) -
00721      $                             DDOT( J-KI-1, T( KI+1, J ), 1,
00722      $                             WORK( KI+1+N ), 1 )
00723 *
00724                      WORK( J+1+N ) = WORK( J+1+N ) -
00725      $                               DDOT( J-KI-1, T( KI+1, J+1 ), 1,
00726      $                               WORK( KI+1+N ), 1 )
00727 *
00728 *                    Solve
00729 *                      [T(J,J)-WR   T(J,J+1)     ]'* X = SCALE*( WORK1 )
00730 *                      [T(J+1,J)    T(J+1,J+1)-WR]             ( WORK2 )
00731 *
00732                      CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ),
00733      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
00734      $                            ZERO, X, 2, SCALE, XNORM, IERR )
00735 *
00736 *                    Scale if necessary
00737 *
00738                      IF( SCALE.NE.ONE )
00739      $                  CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
00740                      WORK( J+N ) = X( 1, 1 )
00741                      WORK( J+1+N ) = X( 2, 1 )
00742 *
00743                      VMAX = MAX( ABS( WORK( J+N ) ),
00744      $                      ABS( WORK( J+1+N ) ), VMAX )
00745                      VCRIT = BIGNUM / VMAX
00746 *
00747                   END IF
00748   170          CONTINUE
00749 *
00750 *              Copy the vector x or Q*x to VL and normalize.
00751 *
00752                IF( .NOT.OVER ) THEN
00753                   CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
00754 *
00755                   II = IDAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
00756                   REMAX = ONE / ABS( VL( II, IS ) )
00757                   CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
00758 *
00759                   DO 180 K = 1, KI - 1
00760                      VL( K, IS ) = ZERO
00761   180             CONTINUE
00762 *
00763                ELSE
00764 *
00765                   IF( KI.LT.N )
00766      $               CALL DGEMV( 'N', N, N-KI, ONE, VL( 1, KI+1 ), LDVL,
00767      $                           WORK( KI+1+N ), 1, WORK( KI+N ),
00768      $                           VL( 1, KI ), 1 )
00769 *
00770                   II = IDAMAX( N, VL( 1, KI ), 1 )
00771                   REMAX = ONE / ABS( VL( II, KI ) )
00772                   CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
00773 *
00774                END IF
00775 *
00776             ELSE
00777 *
00778 *              Complex left eigenvector.
00779 *
00780 *               Initial solve:
00781 *                 ((T(KI,KI)    T(KI,KI+1) )' - (WR - I* WI))*X = 0.
00782 *                 ((T(KI+1,KI) T(KI+1,KI+1))                )
00783 *
00784                IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN
00785                   WORK( KI+N ) = WI / T( KI, KI+1 )
00786                   WORK( KI+1+N2 ) = ONE
00787                ELSE
00788                   WORK( KI+N ) = ONE
00789                   WORK( KI+1+N2 ) = -WI / T( KI+1, KI )
00790                END IF
00791                WORK( KI+1+N ) = ZERO
00792                WORK( KI+N2 ) = ZERO
00793 *
00794 *              Form right-hand side
00795 *
00796                DO 190 K = KI + 2, N
00797                   WORK( K+N ) = -WORK( KI+N )*T( KI, K )
00798                   WORK( K+N2 ) = -WORK( KI+1+N2 )*T( KI+1, K )
00799   190          CONTINUE
00800 *
00801 *              Solve complex quasi-triangular system:
00802 *              ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2
00803 *
00804                VMAX = ONE
00805                VCRIT = BIGNUM
00806 *
00807                JNXT = KI + 2
00808                DO 200 J = KI + 2, N
00809                   IF( J.LT.JNXT )
00810      $               GO TO 200
00811                   J1 = J
00812                   J2 = J
00813                   JNXT = J + 1
00814                   IF( J.LT.N ) THEN
00815                      IF( T( J+1, J ).NE.ZERO ) THEN
00816                         J2 = J + 1
00817                         JNXT = J + 2
00818                      END IF
00819                   END IF
00820 *
00821                   IF( J1.EQ.J2 ) THEN
00822 *
00823 *                    1-by-1 diagonal block
00824 *
00825 *                    Scale if necessary to avoid overflow when
00826 *                    forming the right-hand side elements.
00827 *
00828                      IF( WORK( J ).GT.VCRIT ) THEN
00829                         REC = ONE / VMAX
00830                         CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
00831                         CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
00832                         VMAX = ONE
00833                         VCRIT = BIGNUM
00834                      END IF
00835 *
00836                      WORK( J+N ) = WORK( J+N ) -
00837      $                             DDOT( J-KI-2, T( KI+2, J ), 1,
00838      $                             WORK( KI+2+N ), 1 )
00839                      WORK( J+N2 ) = WORK( J+N2 ) -
00840      $                              DDOT( J-KI-2, T( KI+2, J ), 1,
00841      $                              WORK( KI+2+N2 ), 1 )
00842 *
00843 *                    Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2
00844 *
00845                      CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
00846      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
00847      $                            -WI, X, 2, SCALE, XNORM, IERR )
00848 *
00849 *                    Scale if necessary
00850 *
00851                      IF( SCALE.NE.ONE ) THEN
00852                         CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
00853                         CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
00854                      END IF
00855                      WORK( J+N ) = X( 1, 1 )
00856                      WORK( J+N2 ) = X( 1, 2 )
00857                      VMAX = MAX( ABS( WORK( J+N ) ),
00858      $                      ABS( WORK( J+N2 ) ), VMAX )
00859                      VCRIT = BIGNUM / VMAX
00860 *
00861                   ELSE
00862 *
00863 *                    2-by-2 diagonal block
00864 *
00865 *                    Scale if necessary to avoid overflow when forming
00866 *                    the right-hand side elements.
00867 *
00868                      BETA = MAX( WORK( J ), WORK( J+1 ) )
00869                      IF( BETA.GT.VCRIT ) THEN
00870                         REC = ONE / VMAX
00871                         CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
00872                         CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
00873                         VMAX = ONE
00874                         VCRIT = BIGNUM
00875                      END IF
00876 *
00877                      WORK( J+N ) = WORK( J+N ) -
00878      $                             DDOT( J-KI-2, T( KI+2, J ), 1,
00879      $                             WORK( KI+2+N ), 1 )
00880 *
00881                      WORK( J+N2 ) = WORK( J+N2 ) -
00882      $                              DDOT( J-KI-2, T( KI+2, J ), 1,
00883      $                              WORK( KI+2+N2 ), 1 )
00884 *
00885                      WORK( J+1+N ) = WORK( J+1+N ) -
00886      $                               DDOT( J-KI-2, T( KI+2, J+1 ), 1,
00887      $                               WORK( KI+2+N ), 1 )
00888 *
00889                      WORK( J+1+N2 ) = WORK( J+1+N2 ) -
00890      $                                DDOT( J-KI-2, T( KI+2, J+1 ), 1,
00891      $                                WORK( KI+2+N2 ), 1 )
00892 *
00893 *                    Solve 2-by-2 complex linear equation
00894 *                      ([T(j,j)   T(j,j+1)  ]'-(wr-i*wi)*I)*X = SCALE*B
00895 *                      ([T(j+1,j) T(j+1,j+1)]             )
00896 *
00897                      CALL DLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ),
00898      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
00899      $                            -WI, X, 2, SCALE, XNORM, IERR )
00900 *
00901 *                    Scale if necessary
00902 *
00903                      IF( SCALE.NE.ONE ) THEN
00904                         CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
00905                         CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
00906                      END IF
00907                      WORK( J+N ) = X( 1, 1 )
00908                      WORK( J+N2 ) = X( 1, 2 )
00909                      WORK( J+1+N ) = X( 2, 1 )
00910                      WORK( J+1+N2 ) = X( 2, 2 )
00911                      VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ),
00912      $                      ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ), VMAX )
00913                      VCRIT = BIGNUM / VMAX
00914 *
00915                   END IF
00916   200          CONTINUE
00917 *
00918 *              Copy the vector x or Q*x to VL and normalize.
00919 *
00920                IF( .NOT.OVER ) THEN
00921                   CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
00922                   CALL DCOPY( N-KI+1, WORK( KI+N2 ), 1, VL( KI, IS+1 ),
00923      $                        1 )
00924 *
00925                   EMAX = ZERO
00926                   DO 220 K = KI, N
00927                      EMAX = MAX( EMAX, ABS( VL( K, IS ) )+
00928      $                      ABS( VL( K, IS+1 ) ) )
00929   220             CONTINUE
00930                   REMAX = ONE / EMAX
00931                   CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
00932                   CALL DSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 )
00933 *
00934                   DO 230 K = 1, KI - 1
00935                      VL( K, IS ) = ZERO
00936                      VL( K, IS+1 ) = ZERO
00937   230             CONTINUE
00938                ELSE
00939                   IF( KI.LT.N-1 ) THEN
00940                      CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
00941      $                           LDVL, WORK( KI+2+N ), 1, WORK( KI+N ),
00942      $                           VL( 1, KI ), 1 )
00943                      CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
00944      $                           LDVL, WORK( KI+2+N2 ), 1,
00945      $                           WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
00946                   ELSE
00947                      CALL DSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 )
00948                      CALL DSCAL( N, WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
00949                   END IF
00950 *
00951                   EMAX = ZERO
00952                   DO 240 K = 1, N
00953                      EMAX = MAX( EMAX, ABS( VL( K, KI ) )+
00954      $                      ABS( VL( K, KI+1 ) ) )
00955   240             CONTINUE
00956                   REMAX = ONE / EMAX
00957                   CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
00958                   CALL DSCAL( N, REMAX, VL( 1, KI+1 ), 1 )
00959 *
00960                END IF
00961 *
00962             END IF
00963 *
00964             IS = IS + 1
00965             IF( IP.NE.0 )
00966      $         IS = IS + 1
00967   250       CONTINUE
00968             IF( IP.EQ.-1 )
00969      $         IP = 0
00970             IF( IP.EQ.1 )
00971      $         IP = -1
00972 *
00973   260    CONTINUE
00974 *
00975       END IF
00976 *
00977       RETURN
00978 *
00979 *     End of DTREVC
00980 *
00981       END