Logo Search packages:      
Sourcecode: octave2.9 version File versions

dtrevc.f

      SUBROUTINE DTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
     $                   LDVR, MM, M, WORK, INFO )
*
*  -- LAPACK routine (version 3.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     June 30, 1999
*
*     .. 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.
*
*  The right eigenvector x and the left eigenvector y of T corresponding
*  to an eigenvalue w are defined by:
*
*               T*x = w*x,     y'*T = w*y'
*
*  where y





' denotes the conjugate transpose of the vector y.**  If all eigenvectors are requested, the routine may either return the*  matrices X and/or Y of right or left eigenvectors of T, or the*  products Q*X and/or Q*Y, where Q is an input orthogonal*  matrix. If T was obtained from the real-Schur factorization of an*  original matrix A = Q*T*Q', then Q*X and Q*Y are the matrices of
*  right or left eigenvectors of A.
*
*  T must be in Schur canonical form (as returned by DHSEQR), that is,
*  block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
*  2-by-2 diagonal block has its diagonal elements equal and its
*  off-diagonal elements of opposite sign.  Corresponding to each 2-by-2
*  diagonal block is a complex conjugate pair of eigenvalues and
*  eigenvectors; only one eigenvector of the pair is computed, namely
*  the one corresponding to the eigenvalue with positive imaginary part.
*
*  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,
*                  and backtransform them using the input matrices
*                  supplied in VR and/or VL;
*          = 'S':  compute selected right and/or left eigenvectors,
*                  specified by the logical array SELECT.
*
*  SELECT  (input/output) LOGICAL array, dimension (N)
*          If HOWMNY = 'S', SELECT specifies the eigenvectors to be
*          computed.
*          If HOWMNY = 'A' or 'B', SELECT is not referenced.
*          To select the real eigenvector corresponding to a real
*          eigenvalue w(j), SELECT(j) must be set to .TRUE..  To select
*          the complex eigenvector corresponding to a complex conjugate
*          pair w(j) and w(j+1), either SELECT(j) or SELECT(j+1) must be
*          set to .TRUE.; then on exit SELECT(j) is .TRUE. and
*          SELECT(j+1) is .FALSE..
*
*  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;
*                           VL has the same quasi-lower triangular form
*                           as T







'. If T(i,i) is a real eigenvalue, then*                           the i-th column VL(i) of VL  is its*                           corresponding eigenvector. If T(i:i+1,i:i+1)*                           is a 2-by-2 block whose eigenvalues are*                           complex-conjugate eigenvalues of T, then*                           VL(i)+sqrt(-1)*VL(i+1) is the complex*                           eigenvector corresponding to the eigenvalue*                           with positive real part.*          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.*          If SIDE = 'R



', VL is not referenced.**  LDVL    (input) INTEGER*          The leading dimension of the array VL.  LDVL >= max(1,N) if*          SIDE = 'L' or 'B


'; LDVL >= 1 otherwise.**  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;*                           VR has the same quasi-upper triangular form*                           as T. If T(i,i) is a real eigenvalue, then*                           the i-th column VR(i) of VR  is its*                           corresponding eigenvector. If T(i:i+1,i:i+1)*                           is a 2-by-2 block whose eigenvalues are*                           complex-conjugate eigenvalues of T, then*                           VR(i)+sqrt(-1)*VR(i+1) is the complex*                           eigenvector corresponding to the eigenvalue*                           with positive real part.*          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.*          If SIDE = 'L



', VR is not referenced.**  LDVR    (input) INTEGER*          The leading dimension of the array VR.  LDVR >= max(1,N) if*          SIDE = 'R' or 'B







'; LDVR >= 1 otherwise.**  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.
*
  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 )
*
                  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  Doxygen 1.6.0   Back to index