summaryrefslogtreecommitdiff
path: root/src/lib/lapack/dtrevc.f
diff options
context:
space:
mode:
Diffstat (limited to 'src/lib/lapack/dtrevc.f')
-rw-r--r--src/lib/lapack/dtrevc.f980
1 files changed, 0 insertions, 980 deletions
diff --git a/src/lib/lapack/dtrevc.f b/src/lib/lapack/dtrevc.f
deleted file mode 100644
index a0215f02..00000000
--- a/src/lib/lapack/dtrevc.f
+++ /dev/null
@@ -1,980 +0,0 @@
- 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
generated by cgit (git 2.34.1) at 2025-03-19 07:47:06 +0000