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