SUBROUTINE DTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, \$ LDVR, MM, M, WORK, INFO ) * * -- LAPACK routine (instrumented to count operations, version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * May 7, 2001 * * .. 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( * ) * .. * Common block to return operation count. * OPS is only incremented, OPST is used to accumulate small * contributions to OPS to avoid roundoff error * .. Common blocks .. COMMON / LATIME / OPS, ITCNT * .. * .. Scalars in Common .. DOUBLE PRECISION ITCNT, OPS * .. * * 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, OPST, 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, DLABAD, 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 *** * Initialize OPST = 0 *** * * 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 *** OPS = OPS + N*( N-1 ) / 2 *** * * 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 ) *** * Increment op count, ignoring the possible scaling OPST = OPST + ( 2*( J-1 )+6 ) *** * 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 ) *** * Increment op count, ignoring the possible scaling OPST = OPST + ( 4*( J-2 )+24 ) *** 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 ) *** OPST = OPST + ( 2*KI+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 ) *** OPS = OPS + ( 2*N*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 *** OPST = OPST + 2*( KI-2 ) *** * * 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 ) *** * Increment op count, ignoring the possible scaling OPST = OPST + ( 4*( J-1 )+24 ) *** * 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 ) *** * Increment op count, ignoring the possible scaling OPST = OPST + ( 8*( J-2 )+64 ) *** 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 ) *** OPST = OPST + ( 4*KI+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 ) *** OPS = OPS + ( 4*N*( KI-2 )+6*N+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 *** * Increment op count, ignoring the possible scaling OPST = OPST + ( 2*( J-KI-1 )+6 ) *** * 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 *** * Increment op count, ignoring the possible scaling OPST = OPST + ( 4*( J-KI-1 )+24 ) *** * 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 ) *** OPST = OPST + ( 2*( N-KI+1 )+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 ) *** OPS = OPS + ( 2*N*( N-KI+1 )+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 *** OPST = OPST + 2*( N-KI-1 ) *** * * 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 *** * Increment op count, ignoring the possible scaling OPST = OPST + ( 4*( J-KI-2 )+24 ) *** * 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 *** * Increment op count, ignoring the possible scaling OPST = OPST + ( 8*( J-KI-2 )+64 ) *** * END IF 200 CONTINUE * * Copy the vector x or Q*x to VL and normalize. * 210 CONTINUE 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 ) *** OPST = OPST + ( 4*( N-KI+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 ) *** OPS = OPS + ( 4*N*( N-KI-1 )+6*N+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 *** * Compute final op count OPS = OPS + OPST *** * RETURN * * End of DTREVC * END