SRC/dtrevc.f

   1 *> \brief \b DTREVC
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DTREVC + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrevc.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrevc.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrevc.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
  22 *                          LDVR, MM, M, WORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          HOWMNY, SIDE
  26 *       INTEGER            INFO, LDT, LDVL, LDVR, M, MM, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       LOGICAL            SELECT( * )
  30 *       DOUBLE PRECISION   T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
  31 *      $                   WORK( * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> DTREVC computes some or all of the right and/or left eigenvectors of
  41 *> a real upper quasi-triangular matrix T.
  42 *> Matrices of this type are produced by the Schur factorization of
  43 *> a real general matrix:  A = Q*T*Q**T, as computed by DHSEQR.
  44 *>
  45 *> The right eigenvector x and the left eigenvector y of T corresponding
  46 *> to an eigenvalue w are defined by:
  47 *>
  48 *>    T*x = w*x,     (y**T)*T = w*(y**T)
  49 *>
  50 *> where y**T denotes the transpose of y.
  51 *> The eigenvalues are not input to this routine, but are read directly
  52 *> from the diagonal blocks of T.
  53 *>
  54 *> This routine returns the matrices X and/or Y of right and left
  55 *> eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
  56 *> input matrix.  If Q is the orthogonal factor that reduces a matrix
  57 *> A to Schur form T, then Q*X and Q*Y are the matrices of right and
  58 *> left eigenvectors of A.
  59 *> \endverbatim
  60 *
  61 *  Arguments:
  62 *  ==========
  63 *
  64 *> \param[in] SIDE
  65 *> \verbatim
  66 *>          SIDE is CHARACTER*1
  67 *>          = 'R':  compute right eigenvectors only;
  68 *>          = 'L':  compute left eigenvectors only;
  69 *>          = 'B':  compute both right and left eigenvectors.
  70 *> \endverbatim
  71 *>
  72 *> \param[in] HOWMNY
  73 *> \verbatim
  74 *>          HOWMNY is CHARACTER*1
  75 *>          = 'A':  compute all right and/or left eigenvectors;
  76 *>          = 'B':  compute all right and/or left eigenvectors,
  77 *>                  backtransformed by the matrices in VR and/or VL;
  78 *>          = 'S':  compute selected right and/or left eigenvectors,
  79 *>                  as indicated by the logical array SELECT.
  80 *> \endverbatim
  81 *>
  82 *> \param[in,out] SELECT
  83 *> \verbatim
  84 *>          SELECT is LOGICAL array, dimension (N)
  85 *>          If HOWMNY = 'S', SELECT specifies the eigenvectors to be
  86 *>          computed.
  87 *>          If w(j) is a real eigenvalue, the corresponding real
  88 *>          eigenvector is computed if SELECT(j) is .TRUE..
  89 *>          If w(j) and w(j+1) are the real and imaginary parts of a
  90 *>          complex eigenvalue, the corresponding complex eigenvector is
  91 *>          computed if either SELECT(j) or SELECT(j+1) is .TRUE., and
  92 *>          on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to
  93 *>          .FALSE..
  94 *>          Not referenced if HOWMNY = 'A' or 'B'.
  95 *> \endverbatim
  96 *>
  97 *> \param[in] N
  98 *> \verbatim
  99 *>          N is INTEGER
 100 *>          The order of the matrix T. N >= 0.
 101 *> \endverbatim
 102 *>
 103 *> \param[in] T
 104 *> \verbatim
 105 *>          T is DOUBLE PRECISION array, dimension (LDT,N)
 106 *>          The upper quasi-triangular matrix T in Schur canonical form.
 107 *> \endverbatim
 108 *>
 109 *> \param[in] LDT
 110 *> \verbatim
 111 *>          LDT is INTEGER
 112 *>          The leading dimension of the array T. LDT >= max(1,N).
 113 *> \endverbatim
 114 *>
 115 *> \param[in,out] VL
 116 *> \verbatim
 117 *>          VL is DOUBLE PRECISION array, dimension (LDVL,MM)
 118 *>          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
 119 *>          contain an N-by-N matrix Q (usually the orthogonal matrix Q
 120 *>          of Schur vectors returned by DHSEQR).
 121 *>          On exit, if SIDE = 'L' or 'B', VL contains:
 122 *>          if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
 123 *>          if HOWMNY = 'B', the matrix Q*Y;
 124 *>          if HOWMNY = 'S', the left eigenvectors of T specified by
 125 *>                           SELECT, stored consecutively in the columns
 126 *>                           of VL, in the same order as their
 127 *>                           eigenvalues.
 128 *>          A complex eigenvector corresponding to a complex eigenvalue
 129 *>          is stored in two consecutive columns, the first holding the
 130 *>          real part, and the second the imaginary part.
 131 *>          Not referenced if SIDE = 'R'.
 132 *> \endverbatim
 133 *>
 134 *> \param[in] LDVL
 135 *> \verbatim
 136 *>          LDVL is INTEGER
 137 *>          The leading dimension of the array VL.  LDVL >= 1, and if
 138 *>          SIDE = 'L' or 'B', LDVL >= N.
 139 *> \endverbatim
 140 *>
 141 *> \param[in,out] VR
 142 *> \verbatim
 143 *>          VR is DOUBLE PRECISION array, dimension (LDVR,MM)
 144 *>          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
 145 *>          contain an N-by-N matrix Q (usually the orthogonal matrix Q
 146 *>          of Schur vectors returned by DHSEQR).
 147 *>          On exit, if SIDE = 'R' or 'B', VR contains:
 148 *>          if HOWMNY = 'A', the matrix X of right eigenvectors of T;
 149 *>          if HOWMNY = 'B', the matrix Q*X;
 150 *>          if HOWMNY = 'S', the right eigenvectors of T specified by
 151 *>                           SELECT, stored consecutively in the columns
 152 *>                           of VR, in the same order as their
 153 *>                           eigenvalues.
 154 *>          A complex eigenvector corresponding to a complex eigenvalue
 155 *>          is stored in two consecutive columns, the first holding the
 156 *>          real part and the second the imaginary part.
 157 *>          Not referenced if SIDE = 'L'.
 158 *> \endverbatim
 159 *>
 160 *> \param[in] LDVR
 161 *> \verbatim
 162 *>          LDVR is INTEGER
 163 *>          The leading dimension of the array VR.  LDVR >= 1, and if
 164 *>          SIDE = 'R' or 'B', LDVR >= N.
 165 *> \endverbatim
 166 *>
 167 *> \param[in] MM
 168 *> \verbatim
 169 *>          MM is INTEGER
 170 *>          The number of columns in the arrays VL and/or VR. MM >= M.
 171 *> \endverbatim
 172 *>
 173 *> \param[out] M
 174 *> \verbatim
 175 *>          M is INTEGER
 176 *>          The number of columns in the arrays VL and/or VR actually
 177 *>          used to store the eigenvectors.
 178 *>          If HOWMNY = 'A' or 'B', M is set to N.
 179 *>          Each selected real eigenvector occupies one column and each
 180 *>          selected complex eigenvector occupies two columns.
 181 *> \endverbatim
 182 *>
 183 *> \param[out] WORK
 184 *> \verbatim
 185 *>          WORK is DOUBLE PRECISION array, dimension (3*N)
 186 *> \endverbatim
 187 *>
 188 *> \param[out] INFO
 189 *> \verbatim
 190 *>          INFO is INTEGER
 191 *>          = 0:  successful exit
 192 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 193 *> \endverbatim
 194 *
 195 *  Authors:
 196 *  ========
 197 *
 198 *> \author Univ. of Tennessee
 199 *> \author Univ. of California Berkeley
 200 *> \author Univ. of Colorado Denver
 201 *> \author NAG Ltd.
 202 *
 203 *> \date November 2011
 204 *
 205 *> \ingroup doubleOTHERcomputational
 206 *
 207 *> \par Further Details:
 208 *  =====================
 209 *>
 210 *> \verbatim
 211 *>
 212 *>  The algorithm used in this program is basically backward (forward)
 213 *>  substitution, with scaling to make the the code robust against
 214 *>  possible overflow.
 215 *>
 216 *>  Each eigenvector is normalized so that the element of largest
 217 *>  magnitude has magnitude 1; here the magnitude of a complex number
 218 *>  (x,y) is taken to be |x| + |y|.
 219 *> \endverbatim
 220 *>
 221 *  =====================================================================
 222       SUBROUTINE DTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
 223      $                   LDVR, MM, M, WORK, INFO )
 224 *
 225 *  -- LAPACK computational routine (version 3.4.0) --
 226 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 227 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 228 *     November 2011
 229 *
 230 *     .. Scalar Arguments ..
 231       CHARACTER          HOWMNY, SIDE
 232       INTEGER            INFO, LDT, LDVL, LDVR, M, MM, N
 233 *     ..
 234 *     .. Array Arguments ..
 235       LOGICAL            SELECT( * )
 236       DOUBLE PRECISION   T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
 237      $                   WORK( * )
 238 *     ..
 239 *
 240 *  =====================================================================
 241 *
 242 *     .. Parameters ..
 243       DOUBLE PRECISION   ZERO, ONE
 244       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 245 *     ..
 246 *     .. Local Scalars ..
 247       LOGICAL            ALLV, BOTHV, LEFTV, OVER, PAIR, RIGHTV, SOMEV
 248       INTEGER            I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI, N2
 249       DOUBLE PRECISION   BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
 250      $                   SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR,
 251      $                   XNORM
 252 *     ..
 253 *     .. External Functions ..
 254       LOGICAL            LSAME
 255       INTEGER            IDAMAX
 256       DOUBLE PRECISION   DDOT, DLAMCH
 257       EXTERNAL           LSAME, IDAMAX, DDOT, DLAMCH
 258 *     ..
 259 *     .. External Subroutines ..
 260       EXTERNAL           DAXPY, DCOPY, DGEMV, DLALN2, DSCAL, XERBLA
 261 *     ..
 262 *     .. Intrinsic Functions ..
 263       INTRINSIC          ABS, MAX, SQRT
 264 *     ..
 265 *     .. Local Arrays ..
 266       DOUBLE PRECISION   X( 2, 2 )
 267 *     ..
 268 *     .. Executable Statements ..
 269 *
 270 *     Decode and test the input parameters
 271 *
 272       BOTHV = LSAME( SIDE, 'B' )
 273       RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
 274       LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
 275 *
 276       ALLV = LSAME( HOWMNY, 'A' )
 277       OVER = LSAME( HOWMNY, 'B' )
 278       SOMEV = LSAME( HOWMNY, 'S' )
 279 *
 280       INFO = 0
 281       IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
 282          INFO = -1
 283       ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
 284          INFO = -2
 285       ELSE IF( N.LT.0 ) THEN
 286          INFO = -4
 287       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
 288          INFO = -6
 289       ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
 290          INFO = -8
 291       ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
 292          INFO = -10
 293       ELSE
 294 *
 295 *        Set M to the number of columns required to store the selected
 296 *        eigenvectors, standardize the array SELECT if necessary, and
 297 *        test MM.
 298 *
 299          IF( SOMEV ) THEN
 300             M = 0
 301             PAIR = .FALSE.
 302             DO 10 J = 1, N
 303                IF( PAIR ) THEN
 304                   PAIR = .FALSE.
 305                   SELECT( J ) = .FALSE.
 306                ELSE
 307                   IF( J.LT.N ) THEN
 308                      IF( T( J+1, J ).EQ.ZERO ) THEN
 309                         IF( SELECT( J ) )
 310      $                     M = M + 1
 311                      ELSE
 312                         PAIR = .TRUE.
 313                         IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN
 314                            SELECT( J ) = .TRUE.
 315                            M = M + 2
 316                         END IF
 317                      END IF
 318                   ELSE
 319                      IF( SELECT( N ) )
 320      $                  M = M + 1
 321                   END IF
 322                END IF
 323    10       CONTINUE
 324          ELSE
 325             M = N
 326          END IF
 327 *
 328          IF( MM.LT.M ) THEN
 329             INFO = -11
 330          END IF
 331       END IF
 332       IF( INFO.NE.0 ) THEN
 333          CALL XERBLA( 'DTREVC', -INFO )
 334          RETURN
 335       END IF
 336 *
 337 *     Quick return if possible.
 338 *
 339       IF( N.EQ.0 )
 340      $   RETURN
 341 *
 342 *     Set the constants to control overflow.
 343 *
 344       UNFL = DLAMCH( 'Safe minimum' )
 345       OVFL = ONE / UNFL
 346       CALL DLABAD( UNFL, OVFL )
 347       ULP = DLAMCH( 'Precision' )
 348       SMLNUM = UNFL*( N / ULP )
 349       BIGNUM = ( ONE-ULP ) / SMLNUM
 350 *
 351 *     Compute 1-norm of each column of strictly upper triangular
 352 *     part of T to control overflow in triangular solver.
 353 *
 354       WORK( 1 ) = ZERO
 355       DO 30 J = 2, N
 356          WORK( J ) = ZERO
 357          DO 20 I = 1, J - 1
 358             WORK( J ) = WORK( J ) + ABS( T( I, J ) )
 359    20    CONTINUE
 360    30 CONTINUE
 361 *
 362 *     Index IP is used to specify the real or complex eigenvalue:
 363 *       IP = 0, real eigenvalue,
 364 *            1, first of conjugate complex pair: (wr,wi)
 365 *           -1, second of conjugate complex pair: (wr,wi)
 366 *
 367       N2 = 2*N
 368 *
 369       IF( RIGHTV ) THEN
 370 *
 371 *        Compute right eigenvectors.
 372 *
 373          IP = 0
 374          IS = M
 375          DO 140 KI = N, 1, -1
 376 *
 377             IF( IP.EQ.1 )
 378      $         GO TO 130
 379             IF( KI.EQ.1 )
 380      $         GO TO 40
 381             IF( T( KI, KI-1 ).EQ.ZERO )
 382      $         GO TO 40
 383             IP = -1
 384 *
 385    40       CONTINUE
 386             IF( SOMEV ) THEN
 387                IF( IP.EQ.0 ) THEN
 388                   IF( .NOT.SELECT( KI ) )
 389      $               GO TO 130
 390                ELSE
 391                   IF( .NOT.SELECT( KI-1 ) )
 392      $               GO TO 130
 393                END IF
 394             END IF
 395 *
 396 *           Compute the KI-th eigenvalue (WR,WI).
 397 *
 398             WR = T( KI, KI )
 399             WI = ZERO
 400             IF( IP.NE.0 )
 401      $         WI = SQRT( ABS( T( KI, KI-1 ) ) )*
 402      $              SQRT( ABS( T( KI-1, KI ) ) )
 403             SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
 404 *
 405             IF( IP.EQ.0 ) THEN
 406 *
 407 *              Real right eigenvector
 408 *
 409                WORK( KI+N ) = ONE
 410 *
 411 *              Form right-hand side
 412 *
 413                DO 50 K = 1, KI - 1
 414                   WORK( K+N ) = -T( K, KI )
 415    50          CONTINUE
 416 *
 417 *              Solve the upper quasi-triangular system:
 418 *                 (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK.
 419 *
 420                JNXT = KI - 1
 421                DO 60 J = KI - 1, 1, -1
 422                   IF( J.GT.JNXT )
 423      $               GO TO 60
 424                   J1 = J
 425                   J2 = J
 426                   JNXT = J - 1
 427                   IF( J.GT.1 ) THEN
 428                      IF( T( J, J-1 ).NE.ZERO ) THEN
 429                         J1 = J - 1
 430                         JNXT = J - 2
 431                      END IF
 432                   END IF
 433 *
 434                   IF( J1.EQ.J2 ) THEN
 435 *
 436 *                    1-by-1 diagonal block
 437 *
 438                      CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
 439      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
 440      $                            ZERO, X, 2, SCALE, XNORM, IERR )
 441 *
 442 *                    Scale X(1,1) to avoid overflow when updating
 443 *                    the right-hand side.
 444 *
 445                      IF( XNORM.GT.ONE ) THEN
 446                         IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
 447                            X( 1, 1 ) = X( 1, 1 ) / XNORM
 448                            SCALE = SCALE / XNORM
 449                         END IF
 450                      END IF
 451 *
 452 *                    Scale if necessary
 453 *
 454                      IF( SCALE.NE.ONE )
 455      $                  CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
 456                      WORK( J+N ) = X( 1, 1 )
 457 *
 458 *                    Update right-hand side
 459 *
 460                      CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
 461      $                           WORK( 1+N ), 1 )
 462 *
 463                   ELSE
 464 *
 465 *                    2-by-2 diagonal block
 466 *
 467                      CALL DLALN2( .FALSE., 2, 1, SMIN, ONE,
 468      $                            T( J-1, J-1 ), LDT, ONE, ONE,
 469      $                            WORK( J-1+N ), N, WR, ZERO, X, 2,
 470      $                            SCALE, XNORM, IERR )
 471 *
 472 *                    Scale X(1,1) and X(2,1) to avoid overflow when
 473 *                    updating the right-hand side.
 474 *
 475                      IF( XNORM.GT.ONE ) THEN
 476                         BETA = MAX( WORK( J-1 ), WORK( J ) )
 477                         IF( BETA.GT.BIGNUM / XNORM ) THEN
 478                            X( 1, 1 ) = X( 1, 1 ) / XNORM
 479                            X( 2, 1 ) = X( 2, 1 ) / XNORM
 480                            SCALE = SCALE / XNORM
 481                         END IF
 482                      END IF
 483 *
 484 *                    Scale if necessary
 485 *
 486                      IF( SCALE.NE.ONE )
 487      $                  CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
 488                      WORK( J-1+N ) = X( 1, 1 )
 489                      WORK( J+N ) = X( 2, 1 )
 490 *
 491 *                    Update right-hand side
 492 *
 493                      CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
 494      $                           WORK( 1+N ), 1 )
 495                      CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
 496      $                           WORK( 1+N ), 1 )
 497                   END IF
 498    60          CONTINUE
 499 *
 500 *              Copy the vector x or Q*x to VR and normalize.
 501 *
 502                IF( .NOT.OVER ) THEN
 503                   CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 )
 504 *
 505                   II = IDAMAX( KI, VR( 1, IS ), 1 )
 506                   REMAX = ONE / ABS( VR( II, IS ) )
 507                   CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
 508 *
 509                   DO 70 K = KI + 1, N
 510                      VR( K, IS ) = ZERO
 511    70             CONTINUE
 512                ELSE
 513                   IF( KI.GT.1 )
 514      $               CALL DGEMV( 'N', N, KI-1, ONE, VR, LDVR,
 515      $                           WORK( 1+N ), 1, WORK( KI+N ),
 516      $                           VR( 1, KI ), 1 )
 517 *
 518                   II = IDAMAX( N, VR( 1, KI ), 1 )
 519                   REMAX = ONE / ABS( VR( II, KI ) )
 520                   CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
 521                END IF
 522 *
 523             ELSE
 524 *
 525 *              Complex right eigenvector.
 526 *
 527 *              Initial solve
 528 *                [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.
 529 *                [ (T(KI,KI-1)   T(KI,KI)   )               ]
 530 *
 531                IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN
 532                   WORK( KI-1+N ) = ONE
 533                   WORK( KI+N2 ) = WI / T( KI-1, KI )
 534                ELSE
 535                   WORK( KI-1+N ) = -WI / T( KI, KI-1 )
 536                   WORK( KI+N2 ) = ONE
 537                END IF
 538                WORK( KI+N ) = ZERO
 539                WORK( KI-1+N2 ) = ZERO
 540 *
 541 *              Form right-hand side
 542 *
 543                DO 80 K = 1, KI - 2
 544                   WORK( K+N ) = -WORK( KI-1+N )*T( K, KI-1 )
 545                   WORK( K+N2 ) = -WORK( KI+N2 )*T( K, KI )
 546    80          CONTINUE
 547 *
 548 *              Solve upper quasi-triangular system:
 549 *              (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2)
 550 *
 551                JNXT = KI - 2
 552                DO 90 J = KI - 2, 1, -1
 553                   IF( J.GT.JNXT )
 554      $               GO TO 90
 555                   J1 = J
 556                   J2 = J
 557                   JNXT = J - 1
 558                   IF( J.GT.1 ) THEN
 559                      IF( T( J, J-1 ).NE.ZERO ) THEN
 560                         J1 = J - 1
 561                         JNXT = J - 2
 562                      END IF
 563                   END IF
 564 *
 565                   IF( J1.EQ.J2 ) THEN
 566 *
 567 *                    1-by-1 diagonal block
 568 *
 569                      CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
 570      $                            LDT, ONE, ONE, WORK( J+N ), N, WR, WI,
 571      $                            X, 2, SCALE, XNORM, IERR )
 572 *
 573 *                    Scale X(1,1) and X(1,2) to avoid overflow when
 574 *                    updating the right-hand side.
 575 *
 576                      IF( XNORM.GT.ONE ) THEN
 577                         IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
 578                            X( 1, 1 ) = X( 1, 1 ) / XNORM
 579                            X( 1, 2 ) = X( 1, 2 ) / XNORM
 580                            SCALE = SCALE / XNORM
 581                         END IF
 582                      END IF
 583 *
 584 *                    Scale if necessary
 585 *
 586                      IF( SCALE.NE.ONE ) THEN
 587                         CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
 588                         CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
 589                      END IF
 590                      WORK( J+N ) = X( 1, 1 )
 591                      WORK( J+N2 ) = X( 1, 2 )
 592 *
 593 *                    Update the right-hand side
 594 *
 595                      CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
 596      $                           WORK( 1+N ), 1 )
 597                      CALL DAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1,
 598      $                           WORK( 1+N2 ), 1 )
 599 *
 600                   ELSE
 601 *
 602 *                    2-by-2 diagonal block
 603 *
 604                      CALL DLALN2( .FALSE., 2, 2, SMIN, ONE,
 605      $                            T( J-1, J-1 ), LDT, ONE, ONE,
 606      $                            WORK( J-1+N ), N, WR, WI, X, 2, SCALE,
 607      $                            XNORM, IERR )
 608 *
 609 *                    Scale X to avoid overflow when updating
 610 *                    the right-hand side.
 611 *
 612                      IF( XNORM.GT.ONE ) THEN
 613                         BETA = MAX( WORK( J-1 ), WORK( J ) )
 614                         IF( BETA.GT.BIGNUM / XNORM ) THEN
 615                            REC = ONE / XNORM
 616                            X( 1, 1 ) = X( 1, 1 )*REC
 617                            X( 1, 2 ) = X( 1, 2 )*REC
 618                            X( 2, 1 ) = X( 2, 1 )*REC
 619                            X( 2, 2 ) = X( 2, 2 )*REC
 620                            SCALE = SCALE*REC
 621                         END IF
 622                      END IF
 623 *
 624 *                    Scale if necessary
 625 *
 626                      IF( SCALE.NE.ONE ) THEN
 627                         CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
 628                         CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
 629                      END IF
 630                      WORK( J-1+N ) = X( 1, 1 )
 631                      WORK( J+N ) = X( 2, 1 )
 632                      WORK( J-1+N2 ) = X( 1, 2 )
 633                      WORK( J+N2 ) = X( 2, 2 )
 634 *
 635 *                    Update the right-hand side
 636 *
 637                      CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
 638      $                           WORK( 1+N ), 1 )
 639                      CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
 640      $                           WORK( 1+N ), 1 )
 641                      CALL DAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
 642      $                           WORK( 1+N2 ), 1 )
 643                      CALL DAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1,
 644      $                           WORK( 1+N2 ), 1 )
 645                   END IF
 646    90          CONTINUE
 647 *
 648 *              Copy the vector x or Q*x to VR and normalize.
 649 *
 650                IF( .NOT.OVER ) THEN
 651                   CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 )
 652                   CALL DCOPY( KI, WORK( 1+N2 ), 1, VR( 1, IS ), 1 )
 653 *
 654                   EMAX = ZERO
 655                   DO 100 K = 1, KI
 656                      EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+
 657      $                      ABS( VR( K, IS ) ) )
 658   100             CONTINUE
 659 *
 660                   REMAX = ONE / EMAX
 661                   CALL DSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
 662                   CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
 663 *
 664                   DO 110 K = KI + 1, N
 665                      VR( K, IS-1 ) = ZERO
 666                      VR( K, IS ) = ZERO
 667   110             CONTINUE
 668 *
 669                ELSE
 670 *
 671                   IF( KI.GT.2 ) THEN
 672                      CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
 673      $                           WORK( 1+N ), 1, WORK( KI-1+N ),
 674      $                           VR( 1, KI-1 ), 1 )
 675                      CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
 676      $                           WORK( 1+N2 ), 1, WORK( KI+N2 ),
 677      $                           VR( 1, KI ), 1 )
 678                   ELSE
 679                      CALL DSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 )
 680                      CALL DSCAL( N, WORK( KI+N2 ), VR( 1, KI ), 1 )
 681                   END IF
 682 *
 683                   EMAX = ZERO
 684                   DO 120 K = 1, N
 685                      EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+
 686      $                      ABS( VR( K, KI ) ) )
 687   120             CONTINUE
 688                   REMAX = ONE / EMAX
 689                   CALL DSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
 690                   CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
 691                END IF
 692             END IF
 693 *
 694             IS = IS - 1
 695             IF( IP.NE.0 )
 696      $         IS = IS - 1
 697   130       CONTINUE
 698             IF( IP.EQ.1 )
 699      $         IP = 0
 700             IF( IP.EQ.-1 )
 701      $         IP = 1
 702   140    CONTINUE
 703       END IF
 704 *
 705       IF( LEFTV ) THEN
 706 *
 707 *        Compute left eigenvectors.
 708 *
 709          IP = 0
 710          IS = 1
 711          DO 260 KI = 1, N
 712 *
 713             IF( IP.EQ.-1 )
 714      $         GO TO 250
 715             IF( KI.EQ.N )
 716      $         GO TO 150
 717             IF( T( KI+1, KI ).EQ.ZERO )
 718      $         GO TO 150
 719             IP = 1
 720 *
 721   150       CONTINUE
 722             IF( SOMEV ) THEN
 723                IF( .NOT.SELECT( KI ) )
 724      $            GO TO 250
 725             END IF
 726 *
 727 *           Compute the KI-th eigenvalue (WR,WI).
 728 *
 729             WR = T( KI, KI )
 730             WI = ZERO
 731             IF( IP.NE.0 )
 732      $         WI = SQRT( ABS( T( KI, KI+1 ) ) )*
 733      $              SQRT( ABS( T( KI+1, KI ) ) )
 734             SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
 735 *
 736             IF( IP.EQ.0 ) THEN
 737 *
 738 *              Real left eigenvector.
 739 *
 740                WORK( KI+N ) = ONE
 741 *
 742 *              Form right-hand side
 743 *
 744                DO 160 K = KI + 1, N
 745                   WORK( K+N ) = -T( KI, K )
 746   160          CONTINUE
 747 *
 748 *              Solve the quasi-triangular system:
 749 *                 (T(KI+1:N,KI+1:N) - WR)**T*X = SCALE*WORK
 750 *
 751                VMAX = ONE
 752                VCRIT = BIGNUM
 753 *
 754                JNXT = KI + 1
 755                DO 170 J = KI + 1, N
 756                   IF( J.LT.JNXT )
 757      $               GO TO 170
 758                   J1 = J
 759                   J2 = J
 760                   JNXT = J + 1
 761                   IF( J.LT.N ) THEN
 762                      IF( T( J+1, J ).NE.ZERO ) THEN
 763                         J2 = J + 1
 764                         JNXT = J + 2
 765                      END IF
 766                   END IF
 767 *
 768                   IF( J1.EQ.J2 ) THEN
 769 *
 770 *                    1-by-1 diagonal block
 771 *
 772 *                    Scale if necessary to avoid overflow when forming
 773 *                    the right-hand side.
 774 *
 775                      IF( WORK( J ).GT.VCRIT ) THEN
 776                         REC = ONE / VMAX
 777                         CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
 778                         VMAX = ONE
 779                         VCRIT = BIGNUM
 780                      END IF
 781 *
 782                      WORK( J+N ) = WORK( J+N ) -
 783      $                             DDOT( J-KI-1, T( KI+1, J ), 1,
 784      $                             WORK( KI+1+N ), 1 )
 785 *
 786 *                    Solve (T(J,J)-WR)**T*X = WORK
 787 *
 788                      CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
 789      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
 790      $                            ZERO, X, 2, SCALE, XNORM, IERR )
 791 *
 792 *                    Scale if necessary
 793 *
 794                      IF( SCALE.NE.ONE )
 795      $                  CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
 796                      WORK( J+N ) = X( 1, 1 )
 797                      VMAX = MAX( ABS( WORK( J+N ) ), VMAX )
 798                      VCRIT = BIGNUM / VMAX
 799 *
 800                   ELSE
 801 *
 802 *                    2-by-2 diagonal block
 803 *
 804 *                    Scale if necessary to avoid overflow when forming
 805 *                    the right-hand side.
 806 *
 807                      BETA = MAX( WORK( J ), WORK( J+1 ) )
 808                      IF( BETA.GT.VCRIT ) THEN
 809                         REC = ONE / VMAX
 810                         CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
 811                         VMAX = ONE
 812                         VCRIT = BIGNUM
 813                      END IF
 814 *
 815                      WORK( J+N ) = WORK( J+N ) -
 816      $                             DDOT( J-KI-1, T( KI+1, J ), 1,
 817      $                             WORK( KI+1+N ), 1 )
 818 *
 819                      WORK( J+1+N ) = WORK( J+1+N ) -
 820      $                               DDOT( J-KI-1, T( KI+1, J+1 ), 1,
 821      $                               WORK( KI+1+N ), 1 )
 822 *
 823 *                    Solve
 824 *                      [T(J,J)-WR   T(J,J+1)     ]**T * X = SCALE*( WORK1 )
 825 *                      [T(J+1,J)    T(J+1,J+1)-WR]                ( WORK2 )
 826 *
 827                      CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ),
 828      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
 829      $                            ZERO, X, 2, SCALE, XNORM, IERR )
 830 *
 831 *                    Scale if necessary
 832 *
 833                      IF( SCALE.NE.ONE )
 834      $                  CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
 835                      WORK( J+N ) = X( 1, 1 )
 836                      WORK( J+1+N ) = X( 2, 1 )
 837 *
 838                      VMAX = MAX( ABS( WORK( J+N ) ),
 839      $                      ABS( WORK( J+1+N ) ), VMAX )
 840                      VCRIT = BIGNUM / VMAX
 841 *
 842                   END IF
 843   170          CONTINUE
 844 *
 845 *              Copy the vector x or Q*x to VL and normalize.
 846 *
 847                IF( .NOT.OVER ) THEN
 848                   CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
 849 *
 850                   II = IDAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
 851                   REMAX = ONE / ABS( VL( II, IS ) )
 852                   CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
 853 *
 854                   DO 180 K = 1, KI - 1
 855                      VL( K, IS ) = ZERO
 856   180             CONTINUE
 857 *
 858                ELSE
 859 *
 860                   IF( KI.LT.N )
 861      $               CALL DGEMV( 'N', N, N-KI, ONE, VL( 1, KI+1 ), LDVL,
 862      $                           WORK( KI+1+N ), 1, WORK( KI+N ),
 863      $                           VL( 1, KI ), 1 )
 864 *
 865                   II = IDAMAX( N, VL( 1, KI ), 1 )
 866                   REMAX = ONE / ABS( VL( II, KI ) )
 867                   CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
 868 *
 869                END IF
 870 *
 871             ELSE
 872 *
 873 *              Complex left eigenvector.
 874 *
 875 *               Initial solve:
 876 *                 ((T(KI,KI)    T(KI,KI+1) )**T - (WR - I* WI))*X = 0.
 877 *                 ((T(KI+1,KI) T(KI+1,KI+1))                )
 878 *
 879                IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN
 880                   WORK( KI+N ) = WI / T( KI, KI+1 )
 881                   WORK( KI+1+N2 ) = ONE
 882                ELSE
 883                   WORK( KI+N ) = ONE
 884                   WORK( KI+1+N2 ) = -WI / T( KI+1, KI )
 885                END IF
 886                WORK( KI+1+N ) = ZERO
 887                WORK( KI+N2 ) = ZERO
 888 *
 889 *              Form right-hand side
 890 *
 891                DO 190 K = KI + 2, N
 892                   WORK( K+N ) = -WORK( KI+N )*T( KI, K )
 893                   WORK( K+N2 ) = -WORK( KI+1+N2 )*T( KI+1, K )
 894   190          CONTINUE
 895 *
 896 *              Solve complex quasi-triangular system:
 897 *              ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2
 898 *
 899                VMAX = ONE
 900                VCRIT = BIGNUM
 901 *
 902                JNXT = KI + 2
 903                DO 200 J = KI + 2, N
 904                   IF( J.LT.JNXT )
 905      $               GO TO 200
 906                   J1 = J
 907                   J2 = J
 908                   JNXT = J + 1
 909                   IF( J.LT.N ) THEN
 910                      IF( T( J+1, J ).NE.ZERO ) THEN
 911                         J2 = J + 1
 912                         JNXT = J + 2
 913                      END IF
 914                   END IF
 915 *
 916                   IF( J1.EQ.J2 ) THEN
 917 *
 918 *                    1-by-1 diagonal block
 919 *
 920 *                    Scale if necessary to avoid overflow when
 921 *                    forming the right-hand side elements.
 922 *
 923                      IF( WORK( J ).GT.VCRIT ) THEN
 924                         REC = ONE / VMAX
 925                         CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
 926                         CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
 927                         VMAX = ONE
 928                         VCRIT = BIGNUM
 929                      END IF
 930 *
 931                      WORK( J+N ) = WORK( J+N ) -
 932      $                             DDOT( J-KI-2, T( KI+2, J ), 1,
 933      $                             WORK( KI+2+N ), 1 )
 934                      WORK( J+N2 ) = WORK( J+N2 ) -
 935      $                              DDOT( J-KI-2, T( KI+2, J ), 1,
 936      $                              WORK( KI+2+N2 ), 1 )
 937 *
 938 *                    Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2
 939 *
 940                      CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
 941      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
 942      $                            -WI, X, 2, SCALE, XNORM, IERR )
 943 *
 944 *                    Scale if necessary
 945 *
 946                      IF( SCALE.NE.ONE ) THEN
 947                         CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
 948                         CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
 949                      END IF
 950                      WORK( J+N ) = X( 1, 1 )
 951                      WORK( J+N2 ) = X( 1, 2 )
 952                      VMAX = MAX( ABS( WORK( J+N ) ),
 953      $                      ABS( WORK( J+N2 ) ), VMAX )
 954                      VCRIT = BIGNUM / VMAX
 955 *
 956                   ELSE
 957 *
 958 *                    2-by-2 diagonal block
 959 *
 960 *                    Scale if necessary to avoid overflow when forming
 961 *                    the right-hand side elements.
 962 *
 963                      BETA = MAX( WORK( J ), WORK( J+1 ) )
 964                      IF( BETA.GT.VCRIT ) THEN
 965                         REC = ONE / VMAX
 966                         CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
 967                         CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
 968                         VMAX = ONE
 969                         VCRIT = BIGNUM
 970                      END IF
 971 *
 972                      WORK( J+N ) = WORK( J+N ) -
 973      $                             DDOT( J-KI-2, T( KI+2, J ), 1,
 974      $                             WORK( KI+2+N ), 1 )
 975 *
 976                      WORK( J+N2 ) = WORK( J+N2 ) -
 977      $                              DDOT( J-KI-2, T( KI+2, J ), 1,
 978      $                              WORK( KI+2+N2 ), 1 )
 979 *
 980                      WORK( J+1+N ) = WORK( J+1+N ) -
 981      $                               DDOT( J-KI-2, T( KI+2, J+1 ), 1,
 982      $                               WORK( KI+2+N ), 1 )
 983 *
 984                      WORK( J+1+N2 ) = WORK( J+1+N2 ) -
 985      $                                DDOT( J-KI-2, T( KI+2, J+1 ), 1,
 986      $                                WORK( KI+2+N2 ), 1 )
 987 *
 988 *                    Solve 2-by-2 complex linear equation
 989 *                      ([T(j,j)   T(j,j+1)  ]**T-(wr-i*wi)*I)*X = SCALE*B
 990 *                      ([T(j+1,j) T(j+1,j+1)]               )
 991 *
 992                      CALL DLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ),
 993      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
 994      $                            -WI, X, 2, SCALE, XNORM, IERR )
 995 *
 996 *                    Scale if necessary
 997 *
 998                      IF( SCALE.NE.ONE ) THEN
 999                         CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
1000                         CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
1001                      END IF
1002                      WORK( J+N ) = X( 1, 1 )
1003                      WORK( J+N2 ) = X( 1, 2 )
1004                      WORK( J+1+N ) = X( 2, 1 )
1005                      WORK( J+1+N2 ) = X( 2, 2 )
1006                      VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ),
1007      $                      ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ), VMAX )
1008                      VCRIT = BIGNUM / VMAX
1009 *
1010                   END IF
1011   200          CONTINUE
1012 *
1013 *              Copy the vector x or Q*x to VL and normalize.
1014 *
1015                IF( .NOT.OVER ) THEN
1016                   CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
1017                   CALL DCOPY( N-KI+1, WORK( KI+N2 ), 1, VL( KI, IS+1 ),
1018      $                        1 )
1019 *
1020                   EMAX = ZERO
1021                   DO 220 K = KI, N
1022                      EMAX = MAX( EMAX, ABS( VL( K, IS ) )+
1023      $                      ABS( VL( K, IS+1 ) ) )
1024   220             CONTINUE
1025                   REMAX = ONE / EMAX
1026                   CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
1027                   CALL DSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 )
1028 *
1029                   DO 230 K = 1, KI - 1
1030                      VL( K, IS ) = ZERO
1031                      VL( K, IS+1 ) = ZERO
1032   230             CONTINUE
1033                ELSE
1034                   IF( KI.LT.N-1 ) THEN
1035                      CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
1036      $                           LDVL, WORK( KI+2+N ), 1, WORK( KI+N ),
1037      $                           VL( 1, KI ), 1 )
1038                      CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
1039      $                           LDVL, WORK( KI+2+N2 ), 1,
1040      $                           WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
1041                   ELSE
1042                      CALL DSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 )
1043                      CALL DSCAL( N, WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
1044                   END IF
1045 *
1046                   EMAX = ZERO
1047                   DO 240 K = 1, N
1048                      EMAX = MAX( EMAX, ABS( VL( K, KI ) )+
1049      $                      ABS( VL( K, KI+1 ) ) )
1050   240             CONTINUE
1051                   REMAX = ONE / EMAX
1052                   CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
1053                   CALL DSCAL( N, REMAX, VL( 1, KI+1 ), 1 )
1054 *
1055                END IF
1056 *
1057             END IF
1058 *
1059             IS = IS + 1
1060             IF( IP.NE.0 )
1061      $         IS = IS + 1
1062   250       CONTINUE
1063             IF( IP.EQ.-1 )
1064      $         IP = 0
1065             IF( IP.EQ.1 )
1066      $         IP = -1
1067 *
1068   260    CONTINUE
1069 *
1070       END IF
1071 *
1072       RETURN
1073 *
1074 *     End of DTREVC
1075 *
1076       END