```      SUBROUTINE DTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
\$                   LDVR, MM, M, WORK, INFO )
*
*  -- LAPACK routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
CHARACTER          HOWMNY, SIDE
INTEGER            INFO, LDT, LDVL, LDVR, M, MM, N
*     ..
*     .. Array Arguments ..
LOGICAL            SELECT( * )
DOUBLE PRECISION   T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
\$                   WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  DTREVC computes some or all of the right and/or left eigenvectors of
*  a real upper quasi-triangular matrix T.
*  Matrices of this type are produced by the Schur factorization of
*  a real general matrix:  A = Q*T*Q**T, as computed by DHSEQR.
*
*  The right eigenvector x and the left eigenvector y of T corresponding
*  to an eigenvalue w are defined by:
*
*     T*x = w*x,     (y**H)*T = w*(y**H)
*
*  where y**H denotes the conjugate transpose of y.
*  The eigenvalues are not input to this routine, but are read directly
*  from the diagonal blocks of T.
*
*  This routine returns the matrices X and/or Y of right and left
*  eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
*  input matrix.  If Q is the orthogonal factor that reduces a matrix
*  A to Schur form T, then Q*X and Q*Y are the matrices of right and
*  left eigenvectors of A.
*
*  Arguments
*  =========
*
*  SIDE    (input) CHARACTER*1
*          = 'R':  compute right eigenvectors only;
*          = 'L':  compute left eigenvectors only;
*          = 'B':  compute both right and left eigenvectors.
*
*  HOWMNY  (input) CHARACTER*1
*          = 'A':  compute all right and/or left eigenvectors;
*          = 'B':  compute all right and/or left eigenvectors,
*                  backtransformed by the matrices in VR and/or VL;
*          = 'S':  compute selected right and/or left eigenvectors,
*                  as indicated by the logical array SELECT.
*
*  SELECT  (input/output) LOGICAL array, dimension (N)
*          If HOWMNY = 'S', SELECT specifies the eigenvectors to be
*          computed.
*          If w(j) is a real eigenvalue, the corresponding real
*          eigenvector is computed if SELECT(j) is .TRUE..
*          If w(j) and w(j+1) are the real and imaginary parts of a
*          complex eigenvalue, the corresponding complex eigenvector is
*          computed if either SELECT(j) or SELECT(j+1) is .TRUE., and
*          on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to
*          .FALSE..
*          Not referenced if HOWMNY = 'A' or 'B'.
*
*  N       (input) INTEGER
*          The order of the matrix T. N >= 0.
*
*  T       (input) DOUBLE PRECISION array, dimension (LDT,N)
*          The upper quasi-triangular matrix T in Schur canonical form.
*
*  LDT     (input) INTEGER
*          The leading dimension of the array T. LDT >= max(1,N).
*
*  VL      (input/output) DOUBLE PRECISION array, dimension (LDVL,MM)
*          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
*          contain an N-by-N matrix Q (usually the orthogonal matrix Q
*          of Schur vectors returned by DHSEQR).
*          On exit, if SIDE = 'L' or 'B', VL contains:
*          if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
*          if HOWMNY = 'B', the matrix Q*Y;
*          if HOWMNY = 'S', the left eigenvectors of T specified by
*                           SELECT, stored consecutively in the columns
*                           of VL, in the same order as their
*                           eigenvalues.
*          A complex eigenvector corresponding to a complex eigenvalue
*          is stored in two consecutive columns, the first holding the
*          real part, and the second the imaginary part.
*          Not referenced if SIDE = 'R'.
*
*  LDVL    (input) INTEGER
*          The leading dimension of the array VL.  LDVL >= 1, and if
*          SIDE = 'L' or 'B', LDVL >= N.
*
*  VR      (input/output) DOUBLE PRECISION array, dimension (LDVR,MM)
*          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
*          contain an N-by-N matrix Q (usually the orthogonal matrix Q
*          of Schur vectors returned by DHSEQR).
*          On exit, if SIDE = 'R' or 'B', VR contains:
*          if HOWMNY = 'A', the matrix X of right eigenvectors of T;
*          if HOWMNY = 'B', the matrix Q*X;
*          if HOWMNY = 'S', the right eigenvectors of T specified by
*                           SELECT, stored consecutively in the columns
*                           of VR, in the same order as their
*                           eigenvalues.
*          A complex eigenvector corresponding to a complex eigenvalue
*          is stored in two consecutive columns, the first holding the
*          real part and the second the imaginary part.
*          Not referenced if SIDE = 'L'.
*
*  LDVR    (input) INTEGER
*          The leading dimension of the array VR.  LDVR >= 1, and if
*          SIDE = 'R' or 'B', LDVR >= N.
*
*  MM      (input) INTEGER
*          The number of columns in the arrays VL and/or VR. MM >= M.
*
*  M       (output) INTEGER
*          The number of columns in the arrays VL and/or VR actually
*          used to store the eigenvectors.
*          If HOWMNY = 'A' or 'B', M is set to N.
*          Each selected real eigenvector occupies one column and each
*          selected complex eigenvector occupies two columns.
*
*  WORK    (workspace) DOUBLE PRECISION array, dimension (3*N)
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value
*
*  Further Details
*  ===============
*
*  The algorithm used in this program is basically backward (forward)
*  substitution, with scaling to make the the code robust against
*  possible overflow.
*
*  Each eigenvector is normalized so that the element of largest
*  magnitude has magnitude 1; here the magnitude of a complex number
*  (x,y) is taken to be |x| + |y|.
*
*  =====================================================================
*
*     .. Parameters ..
DOUBLE PRECISION   ZERO, ONE
PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
*     ..
*     .. Local Scalars ..
LOGICAL            ALLV, BOTHV, LEFTV, OVER, PAIR, RIGHTV, SOMEV
INTEGER            I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI, N2
DOUBLE PRECISION   BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
\$                   SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR,
\$                   XNORM
*     ..
*     .. External Functions ..
LOGICAL            LSAME
INTEGER            IDAMAX
DOUBLE PRECISION   DDOT, DLAMCH
EXTERNAL           LSAME, IDAMAX, DDOT, DLAMCH
*     ..
*     .. External Subroutines ..
EXTERNAL           DAXPY, DCOPY, DGEMV, DLALN2, DSCAL, XERBLA
*     ..
*     .. Intrinsic Functions ..
INTRINSIC          ABS, MAX, SQRT
*     ..
*     .. Local Arrays ..
DOUBLE PRECISION   X( 2, 2 )
*     ..
*     .. Executable Statements ..
*
*     Decode and test the input parameters
*
BOTHV = LSAME( SIDE, 'B' )
RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
*
ALLV = LSAME( HOWMNY, 'A' )
OVER = LSAME( HOWMNY, 'B' )
SOMEV = LSAME( HOWMNY, 'S' )
*
INFO = 0
IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
INFO = -1
ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
INFO = -2
ELSE IF( N.LT.0 ) THEN
INFO = -4
ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
INFO = -6
ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
INFO = -8
ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
INFO = -10
ELSE
*
*        Set M to the number of columns required to store the selected
*        eigenvectors, standardize the array SELECT if necessary, and
*        test MM.
*
IF( SOMEV ) THEN
M = 0
PAIR = .FALSE.
DO 10 J = 1, N
IF( PAIR ) THEN
PAIR = .FALSE.
SELECT( J ) = .FALSE.
ELSE
IF( J.LT.N ) THEN
IF( T( J+1, J ).EQ.ZERO ) THEN
IF( SELECT( J ) )
\$                     M = M + 1
ELSE
PAIR = .TRUE.
IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN
SELECT( J ) = .TRUE.
M = M + 2
END IF
END IF
ELSE
IF( SELECT( N ) )
\$                  M = M + 1
END IF
END IF
10       CONTINUE
ELSE
M = N
END IF
*
IF( MM.LT.M ) THEN
INFO = -11
END IF
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DTREVC', -INFO )
RETURN
END IF
*
*     Quick return if possible.
*
IF( N.EQ.0 )
\$   RETURN
*
*     Set the constants to control overflow.
*
UNFL = DLAMCH( 'Safe minimum' )
OVFL = ONE / UNFL
ULP = DLAMCH( 'Precision' )
SMLNUM = UNFL*( N / ULP )
BIGNUM = ( ONE-ULP ) / SMLNUM
*
*     Compute 1-norm of each column of strictly upper triangular
*     part of T to control overflow in triangular solver.
*
WORK( 1 ) = ZERO
DO 30 J = 2, N
WORK( J ) = ZERO
DO 20 I = 1, J - 1
WORK( J ) = WORK( J ) + ABS( T( I, J ) )
20    CONTINUE
30 CONTINUE
*
*     Index IP is used to specify the real or complex eigenvalue:
*       IP = 0, real eigenvalue,
*            1, first of conjugate complex pair: (wr,wi)
*           -1, second of conjugate complex pair: (wr,wi)
*
N2 = 2*N
*
IF( RIGHTV ) THEN
*
*        Compute right eigenvectors.
*
IP = 0
IS = M
DO 140 KI = N, 1, -1
*
IF( IP.EQ.1 )
\$         GO TO 130
IF( KI.EQ.1 )
\$         GO TO 40
IF( T( KI, KI-1 ).EQ.ZERO )
\$         GO TO 40
IP = -1
*
40       CONTINUE
IF( SOMEV ) THEN
IF( IP.EQ.0 ) THEN
IF( .NOT.SELECT( KI ) )
\$               GO TO 130
ELSE
IF( .NOT.SELECT( KI-1 ) )
\$               GO TO 130
END IF
END IF
*
*           Compute the KI-th eigenvalue (WR,WI).
*
WR = T( KI, KI )
WI = ZERO
IF( IP.NE.0 )
\$         WI = SQRT( ABS( T( KI, KI-1 ) ) )*
\$              SQRT( ABS( T( KI-1, KI ) ) )
SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
*
IF( IP.EQ.0 ) THEN
*
*              Real right eigenvector
*
WORK( KI+N ) = ONE
*
*              Form right-hand side
*
DO 50 K = 1, KI - 1
WORK( K+N ) = -T( K, KI )
50          CONTINUE
*
*              Solve the upper quasi-triangular system:
*                 (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK.
*
JNXT = KI - 1
DO 60 J = KI - 1, 1, -1
IF( J.GT.JNXT )
\$               GO TO 60
J1 = J
J2 = J
JNXT = J - 1
IF( J.GT.1 ) THEN
IF( T( J, J-1 ).NE.ZERO ) THEN
J1 = J - 1
JNXT = J - 2
END IF
END IF
*
IF( J1.EQ.J2 ) THEN
*
*                    1-by-1 diagonal block
*
CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
\$                            LDT, ONE, ONE, WORK( J+N ), N, WR,
\$                            ZERO, X, 2, SCALE, XNORM, IERR )
*
*                    Scale X(1,1) to avoid overflow when updating
*                    the right-hand side.
*
IF( XNORM.GT.ONE ) THEN
IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
X( 1, 1 ) = X( 1, 1 ) / XNORM
SCALE = SCALE / XNORM
END IF
END IF
*
*                    Scale if necessary
*
IF( SCALE.NE.ONE )
\$                  CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
WORK( J+N ) = X( 1, 1 )
*
*                    Update right-hand side
*
CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
\$                           WORK( 1+N ), 1 )
*
ELSE
*
*                    2-by-2 diagonal block
*
CALL DLALN2( .FALSE., 2, 1, SMIN, ONE,
\$                            T( J-1, J-1 ), LDT, ONE, ONE,
\$                            WORK( J-1+N ), N, WR, ZERO, X, 2,
\$                            SCALE, XNORM, IERR )
*
*                    Scale X(1,1) and X(2,1) to avoid overflow when
*                    updating the right-hand side.
*
IF( XNORM.GT.ONE ) THEN
BETA = MAX( WORK( J-1 ), WORK( J ) )
IF( BETA.GT.BIGNUM / XNORM ) THEN
X( 1, 1 ) = X( 1, 1 ) / XNORM
X( 2, 1 ) = X( 2, 1 ) / XNORM
SCALE = SCALE / XNORM
END IF
END IF
*
*                    Scale if necessary
*
IF( SCALE.NE.ONE )
\$                  CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
WORK( J-1+N ) = X( 1, 1 )
WORK( J+N ) = X( 2, 1 )
*
*                    Update right-hand side
*
CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
\$                           WORK( 1+N ), 1 )
CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
\$                           WORK( 1+N ), 1 )
END IF
60          CONTINUE
*
*              Copy the vector x or Q*x to VR and normalize.
*
IF( .NOT.OVER ) THEN
CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 )
*
II = IDAMAX( KI, VR( 1, IS ), 1 )
REMAX = ONE / ABS( VR( II, IS ) )
CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
*
DO 70 K = KI + 1, N
VR( K, IS ) = ZERO
70             CONTINUE
ELSE
IF( KI.GT.1 )
\$               CALL DGEMV( 'N', N, KI-1, ONE, VR, LDVR,
\$                           WORK( 1+N ), 1, WORK( KI+N ),
\$                           VR( 1, KI ), 1 )
*
II = IDAMAX( N, VR( 1, KI ), 1 )
REMAX = ONE / ABS( VR( II, KI ) )
CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
END IF
*
ELSE
*
*              Complex right eigenvector.
*
*              Initial solve
*                [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.
*                [ (T(KI,KI-1)   T(KI,KI)   )               ]
*
IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN
WORK( KI-1+N ) = ONE
WORK( KI+N2 ) = WI / T( KI-1, KI )
ELSE
WORK( KI-1+N ) = -WI / T( KI, KI-1 )
WORK( KI+N2 ) = ONE
END IF
WORK( KI+N ) = ZERO
WORK( KI-1+N2 ) = ZERO
*
*              Form right-hand side
*
DO 80 K = 1, KI - 2
WORK( K+N ) = -WORK( KI-1+N )*T( K, KI-1 )
WORK( K+N2 ) = -WORK( KI+N2 )*T( K, KI )
80          CONTINUE
*
*              Solve upper quasi-triangular system:
*              (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2)
*
JNXT = KI - 2
DO 90 J = KI - 2, 1, -1
IF( J.GT.JNXT )
\$               GO TO 90
J1 = J
J2 = J
JNXT = J - 1
IF( J.GT.1 ) THEN
IF( T( J, J-1 ).NE.ZERO ) THEN
J1 = J - 1
JNXT = J - 2
END IF
END IF
*
IF( J1.EQ.J2 ) THEN
*
*                    1-by-1 diagonal block
*
CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
\$                            LDT, ONE, ONE, WORK( J+N ), N, WR, WI,
\$                            X, 2, SCALE, XNORM, IERR )
*
*                    Scale X(1,1) and X(1,2) to avoid overflow when
*                    updating the right-hand side.
*
IF( XNORM.GT.ONE ) THEN
IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
X( 1, 1 ) = X( 1, 1 ) / XNORM
X( 1, 2 ) = X( 1, 2 ) / XNORM
SCALE = SCALE / XNORM
END IF
END IF
*
*                    Scale if necessary
*
IF( SCALE.NE.ONE ) THEN
CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
END IF
WORK( J+N ) = X( 1, 1 )
WORK( J+N2 ) = X( 1, 2 )
*
*                    Update the right-hand side
*
CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
\$                           WORK( 1+N ), 1 )
CALL DAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1,
\$                           WORK( 1+N2 ), 1 )
*
ELSE
*
*                    2-by-2 diagonal block
*
CALL DLALN2( .FALSE., 2, 2, SMIN, ONE,
\$                            T( J-1, J-1 ), LDT, ONE, ONE,
\$                            WORK( J-1+N ), N, WR, WI, X, 2, SCALE,
\$                            XNORM, IERR )
*
*                    Scale X to avoid overflow when updating
*                    the right-hand side.
*
IF( XNORM.GT.ONE ) THEN
BETA = MAX( WORK( J-1 ), WORK( J ) )
IF( BETA.GT.BIGNUM / XNORM ) THEN
REC = ONE / XNORM
X( 1, 1 ) = X( 1, 1 )*REC
X( 1, 2 ) = X( 1, 2 )*REC
X( 2, 1 ) = X( 2, 1 )*REC
X( 2, 2 ) = X( 2, 2 )*REC
SCALE = SCALE*REC
END IF
END IF
*
*                    Scale if necessary
*
IF( SCALE.NE.ONE ) THEN
CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
END IF
WORK( J-1+N ) = X( 1, 1 )
WORK( J+N ) = X( 2, 1 )
WORK( J-1+N2 ) = X( 1, 2 )
WORK( J+N2 ) = X( 2, 2 )
*
*                    Update the right-hand side
*
CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
\$                           WORK( 1+N ), 1 )
CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
\$                           WORK( 1+N ), 1 )
CALL DAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
\$                           WORK( 1+N2 ), 1 )
CALL DAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1,
\$                           WORK( 1+N2 ), 1 )
END IF
90          CONTINUE
*
*              Copy the vector x or Q*x to VR and normalize.
*
IF( .NOT.OVER ) THEN
CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 )
CALL DCOPY( KI, WORK( 1+N2 ), 1, VR( 1, IS ), 1 )
*
EMAX = ZERO
DO 100 K = 1, KI
EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+
\$                      ABS( VR( K, IS ) ) )
100             CONTINUE
*
REMAX = ONE / EMAX
CALL DSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
*
DO 110 K = KI + 1, N
VR( K, IS-1 ) = ZERO
VR( K, IS ) = ZERO
110             CONTINUE
*
ELSE
*
IF( KI.GT.2 ) THEN
CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
\$                           WORK( 1+N ), 1, WORK( KI-1+N ),
\$                           VR( 1, KI-1 ), 1 )
CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
\$                           WORK( 1+N2 ), 1, WORK( KI+N2 ),
\$                           VR( 1, KI ), 1 )
ELSE
CALL DSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 )
CALL DSCAL( N, WORK( KI+N2 ), VR( 1, KI ), 1 )
END IF
*
EMAX = ZERO
DO 120 K = 1, N
EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+
\$                      ABS( VR( K, KI ) ) )
120             CONTINUE
REMAX = ONE / EMAX
CALL DSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
END IF
END IF
*
IS = IS - 1
IF( IP.NE.0 )
\$         IS = IS - 1
130       CONTINUE
IF( IP.EQ.1 )
\$         IP = 0
IF( IP.EQ.-1 )
\$         IP = 1
140    CONTINUE
END IF
*
IF( LEFTV ) THEN
*
*        Compute left eigenvectors.
*
IP = 0
IS = 1
DO 260 KI = 1, N
*
IF( IP.EQ.-1 )
\$         GO TO 250
IF( KI.EQ.N )
\$         GO TO 150
IF( T( KI+1, KI ).EQ.ZERO )
\$         GO TO 150
IP = 1
*
150       CONTINUE
IF( SOMEV ) THEN
IF( .NOT.SELECT( KI ) )
\$            GO TO 250
END IF
*
*           Compute the KI-th eigenvalue (WR,WI).
*
WR = T( KI, KI )
WI = ZERO
IF( IP.NE.0 )
\$         WI = SQRT( ABS( T( KI, KI+1 ) ) )*
\$              SQRT( ABS( T( KI+1, KI ) ) )
SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
*
IF( IP.EQ.0 ) THEN
*
*              Real left eigenvector.
*
WORK( KI+N ) = ONE
*
*              Form right-hand side
*
DO 160 K = KI + 1, N
WORK( K+N ) = -T( KI, K )
160          CONTINUE
*
*              Solve the quasi-triangular system:
*                 (T(KI+1:N,KI+1:N) - WR)'*X = SCALE*WORK
*
VMAX = ONE
VCRIT = BIGNUM
*
JNXT = KI + 1
DO 170 J = KI + 1, N
IF( J.LT.JNXT )
\$               GO TO 170
J1 = J
J2 = J
JNXT = J + 1
IF( J.LT.N ) THEN
IF( T( J+1, J ).NE.ZERO ) THEN
J2 = J + 1
JNXT = J + 2
END IF
END IF
*
IF( J1.EQ.J2 ) THEN
*
*                    1-by-1 diagonal block
*
*                    Scale if necessary to avoid overflow when forming
*                    the right-hand side.
*
IF( WORK( J ).GT.VCRIT ) THEN
REC = ONE / VMAX
CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
VMAX = ONE
VCRIT = BIGNUM
END IF
*
WORK( J+N ) = WORK( J+N ) -
\$                             DDOT( J-KI-1, T( KI+1, J ), 1,
\$                             WORK( KI+1+N ), 1 )
*
*                    Solve (T(J,J)-WR)'*X = WORK
*
CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
\$                            LDT, ONE, ONE, WORK( J+N ), N, WR,
\$                            ZERO, X, 2, SCALE, XNORM, IERR )
*
*                    Scale if necessary
*
IF( SCALE.NE.ONE )
\$                  CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
WORK( J+N ) = X( 1, 1 )
VMAX = MAX( ABS( WORK( J+N ) ), VMAX )
VCRIT = BIGNUM / VMAX
*
ELSE
*
*                    2-by-2 diagonal block
*
*                    Scale if necessary to avoid overflow when forming
*                    the right-hand side.
*
BETA = MAX( WORK( J ), WORK( J+1 ) )
IF( BETA.GT.VCRIT ) THEN
REC = ONE / VMAX
CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
VMAX = ONE
VCRIT = BIGNUM
END IF
*
WORK( J+N ) = WORK( J+N ) -
\$                             DDOT( J-KI-1, T( KI+1, J ), 1,
\$                             WORK( KI+1+N ), 1 )
*
WORK( J+1+N ) = WORK( J+1+N ) -
\$                               DDOT( J-KI-1, T( KI+1, J+1 ), 1,
\$                               WORK( KI+1+N ), 1 )
*
*                    Solve
*                      [T(J,J)-WR   T(J,J+1)     ]'* X = SCALE*( WORK1 )
*                      [T(J+1,J)    T(J+1,J+1)-WR]             ( WORK2 )
*
CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ),
\$                            LDT, ONE, ONE, WORK( J+N ), N, WR,
\$                            ZERO, X, 2, SCALE, XNORM, IERR )
*
*                    Scale if necessary
*
IF( SCALE.NE.ONE )
\$                  CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
WORK( J+N ) = X( 1, 1 )
WORK( J+1+N ) = X( 2, 1 )
*
VMAX = MAX( ABS( WORK( J+N ) ),
\$                      ABS( WORK( J+1+N ) ), VMAX )
VCRIT = BIGNUM / VMAX
*
END IF
170          CONTINUE
*
*              Copy the vector x or Q*x to VL and normalize.
*
IF( .NOT.OVER ) THEN
CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
*
II = IDAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
REMAX = ONE / ABS( VL( II, IS ) )
CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
*
DO 180 K = 1, KI - 1
VL( K, IS ) = ZERO
180             CONTINUE
*
ELSE
*
IF( KI.LT.N )
\$               CALL DGEMV( 'N', N, N-KI, ONE, VL( 1, KI+1 ), LDVL,
\$                           WORK( KI+1+N ), 1, WORK( KI+N ),
\$                           VL( 1, KI ), 1 )
*
II = IDAMAX( N, VL( 1, KI ), 1 )
REMAX = ONE / ABS( VL( II, KI ) )
CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
*
END IF
*
ELSE
*
*              Complex left eigenvector.
*
*               Initial solve:
*                 ((T(KI,KI)    T(KI,KI+1) )' - (WR - I* WI))*X = 0.
*                 ((T(KI+1,KI) T(KI+1,KI+1))                )
*
IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN
WORK( KI+N ) = WI / T( KI, KI+1 )
WORK( KI+1+N2 ) = ONE
ELSE
WORK( KI+N ) = ONE
WORK( KI+1+N2 ) = -WI / T( KI+1, KI )
END IF
WORK( KI+1+N ) = ZERO
WORK( KI+N2 ) = ZERO
*
*              Form right-hand side
*
DO 190 K = KI + 2, N
WORK( K+N ) = -WORK( KI+N )*T( KI, K )
WORK( K+N2 ) = -WORK( KI+1+N2 )*T( KI+1, K )
190          CONTINUE
*
*              Solve complex quasi-triangular system:
*              ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2
*
VMAX = ONE
VCRIT = BIGNUM
*
JNXT = KI + 2
DO 200 J = KI + 2, N
IF( J.LT.JNXT )
\$               GO TO 200
J1 = J
J2 = J
JNXT = J + 1
IF( J.LT.N ) THEN
IF( T( J+1, J ).NE.ZERO ) THEN
J2 = J + 1
JNXT = J + 2
END IF
END IF
*
IF( J1.EQ.J2 ) THEN
*
*                    1-by-1 diagonal block
*
*                    Scale if necessary to avoid overflow when
*                    forming the right-hand side elements.
*
IF( WORK( J ).GT.VCRIT ) THEN
REC = ONE / VMAX
CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
VMAX = ONE
VCRIT = BIGNUM
END IF
*
WORK( J+N ) = WORK( J+N ) -
\$                             DDOT( J-KI-2, T( KI+2, J ), 1,
\$                             WORK( KI+2+N ), 1 )
WORK( J+N2 ) = WORK( J+N2 ) -
\$                              DDOT( J-KI-2, T( KI+2, J ), 1,
\$                              WORK( KI+2+N2 ), 1 )
*
*                    Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2
*
CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
\$                            LDT, ONE, ONE, WORK( J+N ), N, WR,
\$                            -WI, X, 2, SCALE, XNORM, IERR )
*
*                    Scale if necessary
*
IF( SCALE.NE.ONE ) THEN
CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
END IF
WORK( J+N ) = X( 1, 1 )
WORK( J+N2 ) = X( 1, 2 )
VMAX = MAX( ABS( WORK( J+N ) ),
\$                      ABS( WORK( J+N2 ) ), VMAX )
VCRIT = BIGNUM / VMAX
*
ELSE
*
*                    2-by-2 diagonal block
*
*                    Scale if necessary to avoid overflow when forming
*                    the right-hand side elements.
*
BETA = MAX( WORK( J ), WORK( J+1 ) )
IF( BETA.GT.VCRIT ) THEN
REC = ONE / VMAX
CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
VMAX = ONE
VCRIT = BIGNUM
END IF
*
WORK( J+N ) = WORK( J+N ) -
\$                             DDOT( J-KI-2, T( KI+2, J ), 1,
\$                             WORK( KI+2+N ), 1 )
*
WORK( J+N2 ) = WORK( J+N2 ) -
\$                              DDOT( J-KI-2, T( KI+2, J ), 1,
\$                              WORK( KI+2+N2 ), 1 )
*
WORK( J+1+N ) = WORK( J+1+N ) -
\$                               DDOT( J-KI-2, T( KI+2, J+1 ), 1,
\$                               WORK( KI+2+N ), 1 )
*
WORK( J+1+N2 ) = WORK( J+1+N2 ) -
\$                                DDOT( J-KI-2, T( KI+2, J+1 ), 1,
\$                                WORK( KI+2+N2 ), 1 )
*
*                    Solve 2-by-2 complex linear equation
*                      ([T(j,j)   T(j,j+1)  ]'-(wr-i*wi)*I)*X = SCALE*B
*                      ([T(j+1,j) T(j+1,j+1)]             )
*
CALL DLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ),
\$                            LDT, ONE, ONE, WORK( J+N ), N, WR,
\$                            -WI, X, 2, SCALE, XNORM, IERR )
*
*                    Scale if necessary
*
IF( SCALE.NE.ONE ) THEN
CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
END IF
WORK( J+N ) = X( 1, 1 )
WORK( J+N2 ) = X( 1, 2 )
WORK( J+1+N ) = X( 2, 1 )
WORK( J+1+N2 ) = X( 2, 2 )
VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ),
\$                      ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ), VMAX )
VCRIT = BIGNUM / VMAX
*
END IF
200          CONTINUE
*
*              Copy the vector x or Q*x to VL and normalize.
*
IF( .NOT.OVER ) THEN
CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
CALL DCOPY( N-KI+1, WORK( KI+N2 ), 1, VL( KI, IS+1 ),
\$                        1 )
*
EMAX = ZERO
DO 220 K = KI, N
EMAX = MAX( EMAX, ABS( VL( K, IS ) )+
\$                      ABS( VL( K, IS+1 ) ) )
220             CONTINUE
REMAX = ONE / EMAX
CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
CALL DSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 )
*
DO 230 K = 1, KI - 1
VL( K, IS ) = ZERO
VL( K, IS+1 ) = ZERO
230             CONTINUE
ELSE
IF( KI.LT.N-1 ) THEN
CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
\$                           LDVL, WORK( KI+2+N ), 1, WORK( KI+N ),
\$                           VL( 1, KI ), 1 )
CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
\$                           LDVL, WORK( KI+2+N2 ), 1,
\$                           WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
ELSE
CALL DSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 )
CALL DSCAL( N, WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
END IF
*
EMAX = ZERO
DO 240 K = 1, N
EMAX = MAX( EMAX, ABS( VL( K, KI ) )+
\$                      ABS( VL( K, KI+1 ) ) )
240             CONTINUE
REMAX = ONE / EMAX
CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
CALL DSCAL( N, REMAX, VL( 1, KI+1 ), 1 )
*
END IF
*
END IF
*
IS = IS + 1
IF( IP.NE.0 )
\$         IS = IS + 1
250       CONTINUE
IF( IP.EQ.-1 )
\$         IP = 0
IF( IP.EQ.1 )
\$         IP = -1
*
260    CONTINUE
*
END IF
*
RETURN
*
*     End of DTREVC
*
END

```