From 8c8d2f518968ce7057eec6aa5cd5aec8faab861a Mon Sep 17 00:00:00 2001 From: jofret Date: Tue, 28 Apr 2009 07:17:00 +0000 Subject: Moving lapack to right place --- src/lib/lapack/dtrevc.f | 980 ------------------------------------------------ 1 file changed, 980 deletions(-) delete mode 100644 src/lib/lapack/dtrevc.f (limited to 'src/lib/lapack/dtrevc.f') 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 -- cgit