Original Version

1 files changed, 980 insertions, 0 deletions

diff --git a/2.3-1/src/fortran/lapack/dtrevc.f b/2.3-1/src/fortran/lapack/dtrevc.f
new file mode 100644
index 00000000..a0215f02
--- /dev/null
+++ b/2.3-1/src/fortran/lapack/dtrevc.f
@@ -0,0 +1,980 @@
+      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
+      CALL DLABAD( UNFL, OVFL )
+      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