SRC/dtrevc3.f

   1 *> \brief \b DTREVC3
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DTREVC3 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrevc3.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrevc3.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrevc3.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
  22 *                           VR, LDVR, MM, M, WORK, LWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          HOWMNY, SIDE
  26 *       INTEGER            INFO, LDT, LDVL, LDVR, LWORK, 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 *> DTREVC3 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 the vector 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 *>
  60 *> This uses a Level 3 BLAS version of the back transformation.
  61 *> \endverbatim
  62 *
  63 *  Arguments:
  64 *  ==========
  65 *
  66 *> \param[in] SIDE
  67 *> \verbatim
  68 *>          SIDE is CHARACTER*1
  69 *>          = 'R':  compute right eigenvectors only;
  70 *>          = 'L':  compute left eigenvectors only;
  71 *>          = 'B':  compute both right and left eigenvectors.
  72 *> \endverbatim
  73 *>
  74 *> \param[in] HOWMNY
  75 *> \verbatim
  76 *>          HOWMNY is CHARACTER*1
  77 *>          = 'A':  compute all right and/or left eigenvectors;
  78 *>          = 'B':  compute all right and/or left eigenvectors,
  79 *>                  backtransformed by the matrices in VR and/or VL;
  80 *>          = 'S':  compute selected right and/or left eigenvectors,
  81 *>                  as indicated by the logical array SELECT.
  82 *> \endverbatim
  83 *>
  84 *> \param[in,out] SELECT
  85 *> \verbatim
  86 *>          SELECT is LOGICAL array, dimension (N)
  87 *>          If HOWMNY = 'S', SELECT specifies the eigenvectors to be
  88 *>          computed.
  89 *>          If w(j) is a real eigenvalue, the corresponding real
  90 *>          eigenvector is computed if SELECT(j) is .TRUE..
  91 *>          If w(j) and w(j+1) are the real and imaginary parts of a
  92 *>          complex eigenvalue, the corresponding complex eigenvector is
  93 *>          computed if either SELECT(j) or SELECT(j+1) is .TRUE., and
  94 *>          on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to
  95 *>          .FALSE..
  96 *>          Not referenced if HOWMNY = 'A' or 'B'.
  97 *> \endverbatim
  98 *>
  99 *> \param[in] N
 100 *> \verbatim
 101 *>          N is INTEGER
 102 *>          The order of the matrix T. N >= 0.
 103 *> \endverbatim
 104 *>
 105 *> \param[in] T
 106 *> \verbatim
 107 *>          T is DOUBLE PRECISION array, dimension (LDT,N)
 108 *>          The upper quasi-triangular matrix T in Schur canonical form.
 109 *> \endverbatim
 110 *>
 111 *> \param[in] LDT
 112 *> \verbatim
 113 *>          LDT is INTEGER
 114 *>          The leading dimension of the array T. LDT >= max(1,N).
 115 *> \endverbatim
 116 *>
 117 *> \param[in,out] VL
 118 *> \verbatim
 119 *>          VL is DOUBLE PRECISION array, dimension (LDVL,MM)
 120 *>          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
 121 *>          contain an N-by-N matrix Q (usually the orthogonal matrix Q
 122 *>          of Schur vectors returned by DHSEQR).
 123 *>          On exit, if SIDE = 'L' or 'B', VL contains:
 124 *>          if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
 125 *>          if HOWMNY = 'B', the matrix Q*Y;
 126 *>          if HOWMNY = 'S', the left eigenvectors of T specified by
 127 *>                           SELECT, stored consecutively in the columns
 128 *>                           of VL, in the same order as their
 129 *>                           eigenvalues.
 130 *>          A complex eigenvector corresponding to a complex eigenvalue
 131 *>          is stored in two consecutive columns, the first holding the
 132 *>          real part, and the second the imaginary part.
 133 *>          Not referenced if SIDE = 'R'.
 134 *> \endverbatim
 135 *>
 136 *> \param[in] LDVL
 137 *> \verbatim
 138 *>          LDVL is INTEGER
 139 *>          The leading dimension of the array VL.
 140 *>          LDVL >= 1, and if SIDE = 'L' or 'B', LDVL >= N.
 141 *> \endverbatim
 142 *>
 143 *> \param[in,out] VR
 144 *> \verbatim
 145 *>          VR is DOUBLE PRECISION array, dimension (LDVR,MM)
 146 *>          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
 147 *>          contain an N-by-N matrix Q (usually the orthogonal matrix Q
 148 *>          of Schur vectors returned by DHSEQR).
 149 *>          On exit, if SIDE = 'R' or 'B', VR contains:
 150 *>          if HOWMNY = 'A', the matrix X of right eigenvectors of T;
 151 *>          if HOWMNY = 'B', the matrix Q*X;
 152 *>          if HOWMNY = 'S', the right eigenvectors of T specified by
 153 *>                           SELECT, stored consecutively in the columns
 154 *>                           of VR, in the same order as their
 155 *>                           eigenvalues.
 156 *>          A complex eigenvector corresponding to a complex eigenvalue
 157 *>          is stored in two consecutive columns, the first holding the
 158 *>          real part and the second the imaginary part.
 159 *>          Not referenced if SIDE = 'L'.
 160 *> \endverbatim
 161 *>
 162 *> \param[in] LDVR
 163 *> \verbatim
 164 *>          LDVR is INTEGER
 165 *>          The leading dimension of the array VR.
 166 *>          LDVR >= 1, and if SIDE = 'R' or 'B', LDVR >= N.
 167 *> \endverbatim
 168 *>
 169 *> \param[in] MM
 170 *> \verbatim
 171 *>          MM is INTEGER
 172 *>          The number of columns in the arrays VL and/or VR. MM >= M.
 173 *> \endverbatim
 174 *>
 175 *> \param[out] M
 176 *> \verbatim
 177 *>          M is INTEGER
 178 *>          The number of columns in the arrays VL and/or VR actually
 179 *>          used to store the eigenvectors.
 180 *>          If HOWMNY = 'A' or 'B', M is set to N.
 181 *>          Each selected real eigenvector occupies one column and each
 182 *>          selected complex eigenvector occupies two columns.
 183 *> \endverbatim
 184 *>
 185 *> \param[out] WORK
 186 *> \verbatim
 187 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 188 *> \endverbatim
 189 *>
 190 *> \param[in] LWORK
 191 *> \verbatim
 192 *>          LWORK is INTEGER
 193 *>          The dimension of array WORK. LWORK >= max(1,3*N).
 194 *>          For optimum performance, LWORK >= N + 2*N*NB, where NB is
 195 *>          the optimal blocksize.
 196 *>
 197 *>          If LWORK = -1, then a workspace query is assumed; the routine
 198 *>          only calculates the optimal size of the WORK array, returns
 199 *>          this value as the first entry of the WORK array, and no error
 200 *>          message related to LWORK is issued by XERBLA.
 201 *> \endverbatim
 202 *>
 203 *> \param[out] INFO
 204 *> \verbatim
 205 *>          INFO is INTEGER
 206 *>          = 0:  successful exit
 207 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 208 *> \endverbatim
 209 *
 210 *  Authors:
 211 *  ========
 212 *
 213 *> \author Univ. of Tennessee
 214 *> \author Univ. of California Berkeley
 215 *> \author Univ. of Colorado Denver
 216 *> \author NAG Ltd.
 217 *
 218 *> \date November 2011
 219 *
 220 *  @precisions fortran d -> s
 221 *
 222 *> \ingroup doubleOTHERcomputational
 223 *
 224 *> \par Further Details:
 225 *  =====================
 226 *>
 227 *> \verbatim
 228 *>
 229 *>  The algorithm used in this program is basically backward (forward)
 230 *>  substitution, with scaling to make the the code robust against
 231 *>  possible overflow.
 232 *>
 233 *>  Each eigenvector is normalized so that the element of largest
 234 *>  magnitude has magnitude 1; here the magnitude of a complex number
 235 *>  (x,y) is taken to be |x| + |y|.
 236 *> \endverbatim
 237 *>
 238 *  =====================================================================
 239       SUBROUTINE DTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
 240      $                    VR, LDVR, MM, M, WORK, LWORK, INFO )
 241       IMPLICIT NONE
 242 *
 243 *  -- LAPACK computational routine (version 3.4.0) --
 244 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 245 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 246 *     November 2011
 247 *
 248 *     .. Scalar Arguments ..
 249       CHARACTER          HOWMNY, SIDE
 250       INTEGER            INFO, LDT, LDVL, LDVR, LWORK, M, MM, N
 251 *     ..
 252 *     .. Array Arguments ..
 253       LOGICAL            SELECT( * )
 254       DOUBLE PRECISION   T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
 255      $                   WORK( * )
 256 *     ..
 257 *
 258 *  =====================================================================
 259 *
 260 *     .. Parameters ..
 261       DOUBLE PRECISION   ZERO, ONE
 262       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 263       INTEGER            NBMIN, NBMAX
 264       PARAMETER          ( NBMIN = 8, NBMAX = 128 )
 265 *     ..
 266 *     .. Local Scalars ..
 267       LOGICAL            ALLV, BOTHV, LEFTV, LQUERY, OVER, PAIR,
 268      $                   RIGHTV, SOMEV
 269       INTEGER            I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI,
 270      $                   IV, MAXWRK, NB, KI2
 271       DOUBLE PRECISION   BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
 272      $                   SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR,
 273      $                   XNORM
 274 *     ..
 275 *     .. External Functions ..
 276       LOGICAL            LSAME
 277       INTEGER            IDAMAX, ILAENV
 278       DOUBLE PRECISION   DDOT, DLAMCH
 279       EXTERNAL           LSAME, IDAMAX, ILAENV, DDOT, DLAMCH
 280 *     ..
 281 *     .. External Subroutines ..
 282       EXTERNAL           DAXPY, DCOPY, DGEMV, DLALN2, DSCAL, XERBLA
 283 *     ..
 284 *     .. Intrinsic Functions ..
 285       INTRINSIC          ABS, MAX, SQRT
 286 *     ..
 287 *     .. Local Arrays ..
 288       DOUBLE PRECISION   X( 2, 2 )
 289       INTEGER            ISCOMPLEX( NBMAX )
 290 *     ..
 291 *     .. Executable Statements ..
 292 *
 293 *     Decode and test the input parameters
 294 *
 295       BOTHV  = LSAME( SIDE, 'B' )
 296       RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
 297       LEFTV  = LSAME( SIDE, 'L' ) .OR. BOTHV
 298 *
 299       ALLV  = LSAME( HOWMNY, 'A' )
 300       OVER  = LSAME( HOWMNY, 'B' )
 301       SOMEV = LSAME( HOWMNY, 'S' )
 302 *
 303       INFO = 0
 304       NB = ILAENV( 1, 'DTREVC', SIDE // HOWMNY, N, -1, -1, -1 )
 305       MAXWRK = N + 2*N*NB
 306       WORK(1) = MAXWRK
 307       LQUERY = ( LWORK.EQ.-1 )
 308       IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
 309          INFO = -1
 310       ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
 311          INFO = -2
 312       ELSE IF( N.LT.0 ) THEN
 313          INFO = -4
 314       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
 315          INFO = -6
 316       ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
 317          INFO = -8
 318       ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
 319          INFO = -10
 320       ELSE IF( LWORK.LT.MAX( 1, 3*N ) .AND. .NOT.LQUERY ) THEN
 321          INFO = -14
 322       ELSE
 323 *
 324 *        Set M to the number of columns required to store the selected
 325 *        eigenvectors, standardize the array SELECT if necessary, and
 326 *        test MM.
 327 *
 328          IF( SOMEV ) THEN
 329             M = 0
 330             PAIR = .FALSE.
 331             DO 10 J = 1, N
 332                IF( PAIR ) THEN
 333                   PAIR = .FALSE.
 334                   SELECT( J ) = .FALSE.
 335                ELSE
 336                   IF( J.LT.N ) THEN
 337                      IF( T( J+1, J ).EQ.ZERO ) THEN
 338                         IF( SELECT( J ) )
 339      $                     M = M + 1
 340                      ELSE
 341                         PAIR = .TRUE.
 342                         IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN
 343                            SELECT( J ) = .TRUE.
 344                            M = M + 2
 345                         END IF
 346                      END IF
 347                   ELSE
 348                      IF( SELECT( N ) )
 349      $                  M = M + 1
 350                   END IF
 351                END IF
 352    10       CONTINUE
 353          ELSE
 354             M = N
 355          END IF
 356 *
 357          IF( MM.LT.M ) THEN
 358             INFO = -11
 359          END IF
 360       END IF
 361       IF( INFO.NE.0 ) THEN
 362          CALL XERBLA( 'DTREVC3', -INFO )
 363          RETURN
 364       ELSE IF( LQUERY ) THEN
 365          RETURN
 366       END IF
 367 *
 368 *     Quick return if possible.
 369 *
 370       IF( N.EQ.0 )
 371      $   RETURN
 372 *
 373 *     Use blocked version of back-transformation if sufficient workspace.
 374 *     Zero-out the workspace to avoid potential NaN propagation.
 375 *
 376       IF( OVER .AND. LWORK .GE. N + 2*N*NBMIN ) THEN
 377          NB = (LWORK - N) / (2*N)
 378          NB = MIN( NB, NBMAX )
 379          CALL DLASET( 'F', N, 1+2*NB, ZERO, ZERO, WORK, N )
 380       ELSE
 381          NB = 1
 382       END IF
 383 *
 384 *     Set the constants to control overflow.
 385 *
 386       UNFL = DLAMCH( 'Safe minimum' )
 387       OVFL = ONE / UNFL
 388       CALL DLABAD( UNFL, OVFL )
 389       ULP = DLAMCH( 'Precision' )
 390       SMLNUM = UNFL*( N / ULP )
 391       BIGNUM = ( ONE-ULP ) / SMLNUM
 392 *
 393 *     Compute 1-norm of each column of strictly upper triangular
 394 *     part of T to control overflow in triangular solver.
 395 *
 396       WORK( 1 ) = ZERO
 397       DO 30 J = 2, N
 398          WORK( J ) = ZERO
 399          DO 20 I = 1, J - 1
 400             WORK( J ) = WORK( J ) + ABS( T( I, J ) )
 401    20    CONTINUE
 402    30 CONTINUE
 403 *
 404 *     Index IP is used to specify the real or complex eigenvalue:
 405 *       IP = 0, real eigenvalue,
 406 *            1, first  of conjugate complex pair: (wr,wi)
 407 *           -1, second of conjugate complex pair: (wr,wi)
 408 *       ISCOMPLEX array stores IP for each column in current block.
 409 *
 410       IF( RIGHTV ) THEN
 411 *
 412 *        ============================================================
 413 *        Compute right eigenvectors.
 414 *
 415 *        IV is index of column in current block.
 416 *        For complex right vector, uses IV-1 for real part and IV for complex part.
 417 *        Non-blocked version always uses IV=2;
 418 *        blocked     version starts with IV=NB, goes down to 1 or 2.
 419 *        (Note the "0-th" column is used for 1-norms computed above.)
 420          IV = 2
 421          IF( NB.GT.2 ) THEN
 422             IV = NB
 423          END IF
 424
 425          IP = 0
 426          IS = M
 427          DO 140 KI = N, 1, -1
 428             IF( IP.EQ.-1 ) THEN
 429 *              previous iteration (ki+1) was second of conjugate pair,
 430 *              so this ki is first of conjugate pair; skip to end of loop
 431                IP = 1
 432                GO TO 140
 433             ELSE IF( KI.EQ.1 ) THEN
 434 *              last column, so this ki must be real eigenvalue
 435                IP = 0
 436             ELSE IF( T( KI, KI-1 ).EQ.ZERO ) THEN
 437 *              zero on sub-diagonal, so this ki is real eigenvalue
 438                IP = 0
 439             ELSE
 440 *              non-zero on sub-diagonal, so this ki is second of conjugate pair
 441                IP = -1
 442             END IF
 443
 444             IF( SOMEV ) THEN
 445                IF( IP.EQ.0 ) THEN
 446                   IF( .NOT.SELECT( KI ) )
 447      $               GO TO 140
 448                ELSE
 449                   IF( .NOT.SELECT( KI-1 ) )
 450      $               GO TO 140
 451                END IF
 452             END IF
 453 *
 454 *           Compute the KI-th eigenvalue (WR,WI).
 455 *
 456             WR = T( KI, KI )
 457             WI = ZERO
 458             IF( IP.NE.0 )
 459      $         WI = SQRT( ABS( T( KI, KI-1 ) ) )*
 460      $              SQRT( ABS( T( KI-1, KI ) ) )
 461             SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
 462 *
 463             IF( IP.EQ.0 ) THEN
 464 *
 465 *              --------------------------------------------------------
 466 *              Real right eigenvector
 467 *
 468                WORK( KI + IV*N ) = ONE
 469 *
 470 *              Form right-hand side.
 471 *
 472                DO 50 K = 1, KI - 1
 473                   WORK( K + IV*N ) = -T( K, KI )
 474    50          CONTINUE
 475 *
 476 *              Solve upper quasi-triangular system:
 477 *              [ T(1:KI-1,1:KI-1) - WR ]*X = SCALE*WORK.
 478 *
 479                JNXT = KI - 1
 480                DO 60 J = KI - 1, 1, -1
 481                   IF( J.GT.JNXT )
 482      $               GO TO 60
 483                   J1 = J
 484                   J2 = J
 485                   JNXT = J - 1
 486                   IF( J.GT.1 ) THEN
 487                      IF( T( J, J-1 ).NE.ZERO ) THEN
 488                         J1   = J - 1
 489                         JNXT = J - 2
 490                      END IF
 491                   END IF
 492 *
 493                   IF( J1.EQ.J2 ) THEN
 494 *
 495 *                    1-by-1 diagonal block
 496 *
 497                      CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
 498      $                            LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
 499      $                            ZERO, X, 2, SCALE, XNORM, IERR )
 500 *
 501 *                    Scale X(1,1) to avoid overflow when updating
 502 *                    the right-hand side.
 503 *
 504                      IF( XNORM.GT.ONE ) THEN
 505                         IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
 506                            X( 1, 1 ) = X( 1, 1 ) / XNORM
 507                            SCALE = SCALE / XNORM
 508                         END IF
 509                      END IF
 510 *
 511 *                    Scale if necessary
 512 *
 513                      IF( SCALE.NE.ONE )
 514      $                  CALL DSCAL( KI, SCALE, WORK( 1+IV*N ), 1 )
 515                      WORK( J+IV*N ) = X( 1, 1 )
 516 *
 517 *                    Update right-hand side
 518 *
 519                      CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
 520      $                           WORK( 1+IV*N ), 1 )
 521 *
 522                   ELSE
 523 *
 524 *                    2-by-2 diagonal block
 525 *
 526                      CALL DLALN2( .FALSE., 2, 1, SMIN, ONE,
 527      $                            T( J-1, J-1 ), LDT, ONE, ONE,
 528      $                            WORK( J-1+IV*N ), N, WR, ZERO, X, 2,
 529      $                            SCALE, XNORM, IERR )
 530 *
 531 *                    Scale X(1,1) and X(2,1) to avoid overflow when
 532 *                    updating the right-hand side.
 533 *
 534                      IF( XNORM.GT.ONE ) THEN
 535                         BETA = MAX( WORK( J-1 ), WORK( J ) )
 536                         IF( BETA.GT.BIGNUM / XNORM ) THEN
 537                            X( 1, 1 ) = X( 1, 1 ) / XNORM
 538                            X( 2, 1 ) = X( 2, 1 ) / XNORM
 539                            SCALE = SCALE / XNORM
 540                         END IF
 541                      END IF
 542 *
 543 *                    Scale if necessary
 544 *
 545                      IF( SCALE.NE.ONE )
 546      $                  CALL DSCAL( KI, SCALE, WORK( 1+IV*N ), 1 )
 547                      WORK( J-1+IV*N ) = X( 1, 1 )
 548                      WORK( J  +IV*N ) = X( 2, 1 )
 549 *
 550 *                    Update right-hand side
 551 *
 552                      CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
 553      $                           WORK( 1+IV*N ), 1 )
 554                      CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
 555      $                           WORK( 1+IV*N ), 1 )
 556                   END IF
 557    60          CONTINUE
 558 *
 559 *              Copy the vector x or Q*x to VR and normalize.
 560 *
 561                IF( .NOT.OVER ) THEN
 562 *                 ------------------------------
 563 *                 no back-transform: copy x to VR and normalize.
 564                   CALL DCOPY( KI, WORK( 1 + IV*N ), 1, VR( 1, IS ), 1 )
 565 *
 566                   II = IDAMAX( KI, VR( 1, IS ), 1 )
 567                   REMAX = ONE / ABS( VR( II, IS ) )
 568                   CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
 569 *
 570                   DO 70 K = KI + 1, N
 571                      VR( K, IS ) = ZERO
 572    70             CONTINUE
 573 *
 574                ELSE IF( NB.EQ.1 ) THEN
 575 *                 ------------------------------
 576 *                 version 1: back-transform each vector with GEMV, Q*x.
 577                   IF( KI.GT.1 )
 578      $               CALL DGEMV( 'N', N, KI-1, ONE, VR, LDVR,
 579      $                           WORK( 1 + IV*N ), 1, WORK( KI + IV*N ),
 580      $                           VR( 1, KI ), 1 )
 581 *
 582                   II = IDAMAX( N, VR( 1, KI ), 1 )
 583                   REMAX = ONE / ABS( VR( II, KI ) )
 584                   CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
 585 *
 586                ELSE
 587 *                 ------------------------------
 588 *                 version 2: back-transform block of vectors with GEMM
 589 *                 zero out below vector
 590                   DO K = KI + 1, N
 591                      WORK( K + IV*N ) = ZERO
 592                   END DO
 593                   ISCOMPLEX( IV ) = IP
 594 *                 back-transform and normalization is done below
 595                END IF
 596             ELSE
 597 *
 598 *              --------------------------------------------------------
 599 *              Complex right eigenvector.
 600 *
 601 *              Initial solve
 602 *              [ ( T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I*WI) ]*X = 0.
 603 *              [ ( T(KI,  KI-1) T(KI,  KI) )               ]
 604 *
 605                IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN
 606                   WORK( KI-1 + (IV-1)*N ) = ONE
 607                   WORK( KI   + (IV  )*N ) = WI / T( KI-1, KI )
 608                ELSE
 609                   WORK( KI-1 + (IV-1)*N ) = -WI / T( KI, KI-1 )
 610                   WORK( KI   + (IV  )*N ) = ONE
 611                END IF
 612                WORK( KI   + (IV-1)*N ) = ZERO
 613                WORK( KI-1 + (IV  )*N ) = ZERO
 614 *
 615 *              Form right-hand side.
 616 *
 617                DO 80 K = 1, KI - 2
 618                   WORK( K+(IV-1)*N ) = -WORK( KI-1+(IV-1)*N )*T(K,KI-1)
 619                   WORK( K+(IV  )*N ) = -WORK( KI  +(IV  )*N )*T(K,KI  )
 620    80          CONTINUE
 621 *
 622 *              Solve upper quasi-triangular system:
 623 *              [ T(1:KI-2,1:KI-2) - (WR+i*WI) ]*X = SCALE*(WORK+i*WORK2)
 624 *
 625                JNXT = KI - 2
 626                DO 90 J = KI - 2, 1, -1
 627                   IF( J.GT.JNXT )
 628      $               GO TO 90
 629                   J1 = J
 630                   J2 = J
 631                   JNXT = J - 1
 632                   IF( J.GT.1 ) THEN
 633                      IF( T( J, J-1 ).NE.ZERO ) THEN
 634                         J1   = J - 1
 635                         JNXT = J - 2
 636                      END IF
 637                   END IF
 638 *
 639                   IF( J1.EQ.J2 ) THEN
 640 *
 641 *                    1-by-1 diagonal block
 642 *
 643                      CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
 644      $                            LDT, ONE, ONE, WORK( J+(IV-1)*N ), N,
 645      $                            WR, WI, X, 2, SCALE, XNORM, IERR )
 646 *
 647 *                    Scale X(1,1) and X(1,2) to avoid overflow when
 648 *                    updating the right-hand side.
 649 *
 650                      IF( XNORM.GT.ONE ) THEN
 651                         IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
 652                            X( 1, 1 ) = X( 1, 1 ) / XNORM
 653                            X( 1, 2 ) = X( 1, 2 ) / XNORM
 654                            SCALE = SCALE / XNORM
 655                         END IF
 656                      END IF
 657 *
 658 *                    Scale if necessary
 659 *
 660                      IF( SCALE.NE.ONE ) THEN
 661                         CALL DSCAL( KI, SCALE, WORK( 1+(IV-1)*N ), 1 )
 662                         CALL DSCAL( KI, SCALE, WORK( 1+(IV  )*N ), 1 )
 663                      END IF
 664                      WORK( J+(IV-1)*N ) = X( 1, 1 )
 665                      WORK( J+(IV  )*N ) = X( 1, 2 )
 666 *
 667 *                    Update the right-hand side
 668 *
 669                      CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
 670      $                           WORK( 1+(IV-1)*N ), 1 )
 671                      CALL DAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1,
 672      $                           WORK( 1+(IV  )*N ), 1 )
 673 *
 674                   ELSE
 675 *
 676 *                    2-by-2 diagonal block
 677 *
 678                      CALL DLALN2( .FALSE., 2, 2, SMIN, ONE,
 679      $                            T( J-1, J-1 ), LDT, ONE, ONE,
 680      $                            WORK( J-1+(IV-1)*N ), N, WR, WI, X, 2,
 681      $                            SCALE, XNORM, IERR )
 682 *
 683 *                    Scale X to avoid overflow when updating
 684 *                    the right-hand side.
 685 *
 686                      IF( XNORM.GT.ONE ) THEN
 687                         BETA = MAX( WORK( J-1 ), WORK( J ) )
 688                         IF( BETA.GT.BIGNUM / XNORM ) THEN
 689                            REC = ONE / XNORM
 690                            X( 1, 1 ) = X( 1, 1 )*REC
 691                            X( 1, 2 ) = X( 1, 2 )*REC
 692                            X( 2, 1 ) = X( 2, 1 )*REC
 693                            X( 2, 2 ) = X( 2, 2 )*REC
 694                            SCALE = SCALE*REC
 695                         END IF
 696                      END IF
 697 *
 698 *                    Scale if necessary
 699 *
 700                      IF( SCALE.NE.ONE ) THEN
 701                         CALL DSCAL( KI, SCALE, WORK( 1+(IV-1)*N ), 1 )
 702                         CALL DSCAL( KI, SCALE, WORK( 1+(IV  )*N ), 1 )
 703                      END IF
 704                      WORK( J-1+(IV-1)*N ) = X( 1, 1 )
 705                      WORK( J  +(IV-1)*N ) = X( 2, 1 )
 706                      WORK( J-1+(IV  )*N ) = X( 1, 2 )
 707                      WORK( J  +(IV  )*N ) = X( 2, 2 )
 708 *
 709 *                    Update the right-hand side
 710 *
 711                      CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
 712      $                           WORK( 1+(IV-1)*N   ), 1 )
 713                      CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
 714      $                           WORK( 1+(IV-1)*N   ), 1 )
 715                      CALL DAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
 716      $                           WORK( 1+(IV  )*N ), 1 )
 717                      CALL DAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1,
 718      $                           WORK( 1+(IV  )*N ), 1 )
 719                   END IF
 720    90          CONTINUE
 721 *
 722 *              Copy the vector x or Q*x to VR and normalize.
 723 *
 724                IF( .NOT.OVER ) THEN
 725 *                 ------------------------------
 726 *                 no back-transform: copy x to VR and normalize.
 727                   CALL DCOPY( KI, WORK( 1+(IV-1)*N ), 1, VR(1,IS-1), 1 )
 728                   CALL DCOPY( KI, WORK( 1+(IV  )*N ), 1, VR(1,IS  ), 1 )
 729 *
 730                   EMAX = ZERO
 731                   DO 100 K = 1, KI
 732                      EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+
 733      $                                 ABS( VR( K, IS   ) ) )
 734   100             CONTINUE
 735                   REMAX = ONE / EMAX
 736                   CALL DSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
 737                   CALL DSCAL( KI, REMAX, VR( 1, IS   ), 1 )
 738 *
 739                   DO 110 K = KI + 1, N
 740                      VR( K, IS-1 ) = ZERO
 741                      VR( K, IS   ) = ZERO
 742   110             CONTINUE
 743 *
 744                ELSE IF( NB.EQ.1 ) THEN
 745 *                 ------------------------------
 746 *                 version 1: back-transform each vector with GEMV, Q*x.
 747                   IF( KI.GT.2 ) THEN
 748                      CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
 749      $                           WORK( 1    + (IV-1)*N ), 1,
 750      $                           WORK( KI-1 + (IV-1)*N ), VR(1,KI-1), 1)
 751                      CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
 752      $                           WORK( 1  + (IV)*N ), 1,
 753      $                           WORK( KI + (IV)*N ), VR( 1, KI ), 1 )
 754                   ELSE
 755                      CALL DSCAL( N, WORK(KI-1+(IV-1)*N), VR(1,KI-1), 1)
 756                      CALL DSCAL( N, WORK(KI  +(IV  )*N), VR(1,KI  ), 1)
 757                   END IF
 758 *
 759                   EMAX = ZERO
 760                   DO 120 K = 1, N
 761                      EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+
 762      $                                 ABS( VR( K, KI   ) ) )
 763   120             CONTINUE
 764                   REMAX = ONE / EMAX
 765                   CALL DSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
 766                   CALL DSCAL( N, REMAX, VR( 1, KI   ), 1 )
 767 *
 768                ELSE
 769 *                 ------------------------------
 770 *                 version 2: back-transform block of vectors with GEMM
 771 *                 zero out below vector
 772                   DO K = KI + 1, N
 773                      WORK( K + (IV-1)*N ) = ZERO
 774                      WORK( K + (IV  )*N ) = ZERO
 775                   END DO
 776                   ISCOMPLEX( IV-1 ) = -IP
 777                   ISCOMPLEX( IV   ) =  IP
 778                   IV = IV - 1
 779 *                 back-transform and normalization is done below
 780                END IF
 781             END IF
 782
 783             IF( NB.GT.1 ) THEN
 784 *              --------------------------------------------------------
 785 *              Blocked version of back-transform
 786 *              For complex case, KI2 includes both vectors (KI-1 and KI)
 787                IF( IP.EQ.0 ) THEN
 788                   KI2 = KI
 789                ELSE
 790                   KI2 = KI - 1
 791                END IF
 792
 793 *              Columns IV:NB of work are valid vectors.
 794 *              When the number of vectors stored reaches NB-1 or NB,
 795 *              or if this was last vector, do the GEMM
 796                IF( (IV.LE.2) .OR. (KI2.EQ.1) ) THEN
 797                   CALL DGEMM( 'N', 'N', N, NB-IV+1, KI2+NB-IV, ONE,
 798      $                        VR, LDVR,
 799      $                        WORK( 1 + (IV)*N    ), N,
 800      $                        ZERO,
 801      $                        WORK( 1 + (NB+IV)*N ), N )
 802 *                 normalize vectors
 803                   DO K = IV, NB
 804                      IF( ISCOMPLEX(K).EQ.0 ) THEN
 805 *                       real eigenvector
 806                         II = IDAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
 807                         REMAX = ONE / ABS( WORK( II + (NB+K)*N ) )
 808                      ELSE IF( ISCOMPLEX(K).EQ.1 ) THEN
 809 *                       first eigenvector of conjugate pair
 810                         EMAX = ZERO
 811                         DO II = 1, N
 812                            EMAX = MAX( EMAX,
 813      $                                 ABS( WORK( II + (NB+K  )*N ) )+
 814      $                                 ABS( WORK( II + (NB+K+1)*N ) ) )
 815                         END DO
 816                         REMAX = ONE / EMAX
 817 *                    else if ISCOMPLEX(K).EQ.-1
 818 *                       second eigenvector of conjugate pair
 819 *                       reuse same REMAX as previous K
 820                      END IF
 821                      CALL DSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
 822                   END DO
 823                   CALL DLACPY( 'F', N, NB-IV+1,
 824      $                         WORK( 1 + (NB+IV)*N ), N,
 825      $                         VR( 1, KI2 ), LDVR )
 826                   IV = NB
 827                ELSE
 828                   IV = IV - 1
 829                END IF
 830             END IF ! blocked back-transform
 831 *
 832             IS = IS - 1
 833             IF( IP.NE.0 )
 834      $         IS = IS - 1
 835   140    CONTINUE
 836       END IF
 837
 838       IF( LEFTV ) THEN
 839 *
 840 *        ============================================================
 841 *        Compute left eigenvectors.
 842 *
 843 *        IV is index of column in current block.
 844 *        For complex left vector, uses IV for real part and IV+1 for complex part.
 845 *        Non-blocked version always uses IV=1;
 846 *        blocked     version starts with IV=1, goes up to NB-1 or NB.
 847 *        (Note the "0-th" column is used for 1-norms computed above.)
 848          IV = 1
 849          IP = 0
 850          IS = 1
 851          DO 260 KI = 1, N
 852             IF( IP.EQ.1 ) THEN
 853 *              previous iteration (ki-1) was first of conjugate pair,
 854 *              so this ki is second of conjugate pair; skip to end of loop
 855                IP = -1
 856                GO TO 260
 857             ELSE IF( KI.EQ.N ) THEN
 858 *              last column, so this ki must be real eigenvalue
 859                IP = 0
 860             ELSE IF( T( KI+1, KI ).EQ.ZERO ) THEN
 861 *              zero on sub-diagonal, so this ki is real eigenvalue
 862                IP = 0
 863             ELSE
 864 *              non-zero on sub-diagonal, so this ki is first of conjugate pair
 865                IP = 1
 866             END IF
 867 *
 868             IF( SOMEV ) THEN
 869                IF( .NOT.SELECT( KI ) )
 870      $            GO TO 260
 871             END IF
 872 *
 873 *           Compute the KI-th eigenvalue (WR,WI).
 874 *
 875             WR = T( KI, KI )
 876             WI = ZERO
 877             IF( IP.NE.0 )
 878      $         WI = SQRT( ABS( T( KI, KI+1 ) ) )*
 879      $              SQRT( ABS( T( KI+1, KI ) ) )
 880             SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
 881 *
 882             IF( IP.EQ.0 ) THEN
 883 *
 884 *              --------------------------------------------------------
 885 *              Real left eigenvector
 886 *
 887                WORK( KI + IV*N ) = ONE
 888 *
 889 *              Form right-hand side.
 890 *
 891                DO 160 K = KI + 1, N
 892                   WORK( K + IV*N ) = -T( KI, K )
 893   160          CONTINUE
 894 *
 895 *              Solve transposed quasi-triangular system:
 896 *              [ T(KI+1:N,KI+1:N) - WR ]**T * X = SCALE*WORK
 897 *
 898                VMAX = ONE
 899                VCRIT = BIGNUM
 900 *
 901                JNXT = KI + 1
 902                DO 170 J = KI + 1, N
 903                   IF( J.LT.JNXT )
 904      $               GO TO 170
 905                   J1 = J
 906                   J2 = J
 907                   JNXT = J + 1
 908                   IF( J.LT.N ) THEN
 909                      IF( T( J+1, J ).NE.ZERO ) THEN
 910                         J2 = J + 1
 911                         JNXT = J + 2
 912                      END IF
 913                   END IF
 914 *
 915                   IF( J1.EQ.J2 ) THEN
 916 *
 917 *                    1-by-1 diagonal block
 918 *
 919 *                    Scale if necessary to avoid overflow when forming
 920 *                    the right-hand side.
 921 *
 922                      IF( WORK( J ).GT.VCRIT ) THEN
 923                         REC = ONE / VMAX
 924                         CALL DSCAL( N-KI+1, REC, WORK( KI+IV*N ), 1 )
 925                         VMAX = ONE
 926                         VCRIT = BIGNUM
 927                      END IF
 928 *
 929                      WORK( J+IV*N ) = WORK( J+IV*N ) -
 930      $                                DDOT( J-KI-1, T( KI+1, J ), 1,
 931      $                                      WORK( KI+1+IV*N ), 1 )
 932 *
 933 *                    Solve [ T(J,J) - WR ]**T * X = WORK
 934 *
 935                      CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
 936      $                            LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
 937      $                            ZERO, X, 2, SCALE, XNORM, IERR )
 938 *
 939 *                    Scale if necessary
 940 *
 941                      IF( SCALE.NE.ONE )
 942      $                  CALL DSCAL( N-KI+1, SCALE, WORK( KI+IV*N ), 1 )
 943                      WORK( J+IV*N ) = X( 1, 1 )
 944                      VMAX = MAX( ABS( WORK( J+IV*N ) ), VMAX )
 945                      VCRIT = BIGNUM / VMAX
 946 *
 947                   ELSE
 948 *
 949 *                    2-by-2 diagonal block
 950 *
 951 *                    Scale if necessary to avoid overflow when forming
 952 *                    the right-hand side.
 953 *
 954                      BETA = MAX( WORK( J ), WORK( J+1 ) )
 955                      IF( BETA.GT.VCRIT ) THEN
 956                         REC = ONE / VMAX
 957                         CALL DSCAL( N-KI+1, REC, WORK( KI+IV*N ), 1 )
 958                         VMAX = ONE
 959                         VCRIT = BIGNUM
 960                      END IF
 961 *
 962                      WORK( J+IV*N ) = WORK( J+IV*N ) -
 963      $                                DDOT( J-KI-1, T( KI+1, J ), 1,
 964      $                                      WORK( KI+1+IV*N ), 1 )
 965 *
 966                      WORK( J+1+IV*N ) = WORK( J+1+IV*N ) -
 967      $                                  DDOT( J-KI-1, T( KI+1, J+1 ), 1,
 968      $                                        WORK( KI+1+IV*N ), 1 )
 969 *
 970 *                    Solve
 971 *                    [ T(J,J)-WR   T(J,J+1)      ]**T * X = SCALE*( WORK1 )
 972 *                    [ T(J+1,J)    T(J+1,J+1)-WR ]                ( WORK2 )
 973 *
 974                      CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ),
 975      $                            LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
 976      $                            ZERO, X, 2, SCALE, XNORM, IERR )
 977 *
 978 *                    Scale if necessary
 979 *
 980                      IF( SCALE.NE.ONE )
 981      $                  CALL DSCAL( N-KI+1, SCALE, WORK( KI+IV*N ), 1 )
 982                      WORK( J  +IV*N ) = X( 1, 1 )
 983                      WORK( J+1+IV*N ) = X( 2, 1 )
 984 *
 985                      VMAX = MAX( ABS( WORK( J  +IV*N ) ),
 986      $                           ABS( WORK( J+1+IV*N ) ), VMAX )
 987                      VCRIT = BIGNUM / VMAX
 988 *
 989                   END IF
 990   170          CONTINUE
 991 *
 992 *              Copy the vector x or Q*x to VL and normalize.
 993 *
 994                IF( .NOT.OVER ) THEN
 995 *                 ------------------------------
 996 *                 no back-transform: copy x to VL and normalize.
 997                   CALL DCOPY( N-KI+1, WORK( KI + IV*N ), 1,
 998      $                                VL( KI, IS ), 1 )
 999 *
1000                   II = IDAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
1001                   REMAX = ONE / ABS( VL( II, IS ) )
1002                   CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
1003 *
1004                   DO 180 K = 1, KI - 1
1005                      VL( K, IS ) = ZERO
1006   180             CONTINUE
1007 *
1008                ELSE IF( NB.EQ.1 ) THEN
1009 *                 ------------------------------
1010 *                 version 1: back-transform each vector with GEMV, Q*x.
1011                   IF( KI.LT.N )
1012      $               CALL DGEMV( 'N', N, N-KI, ONE,
1013      $                           VL( 1, KI+1 ), LDVL,
1014      $                           WORK( KI+1 + IV*N ), 1,
1015      $                           WORK( KI   + IV*N ), VL( 1, KI ), 1 )
1016 *
1017                   II = IDAMAX( N, VL( 1, KI ), 1 )
1018                   REMAX = ONE / ABS( VL( II, KI ) )
1019                   CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
1020 *
1021                ELSE
1022 *                 ------------------------------
1023 *                 version 2: back-transform block of vectors with GEMM
1024 *                 zero out above vector
1025 *                 could go from KI-NV+1 to KI-1
1026                   DO K = 1, KI - 1
1027                      WORK( K + IV*N ) = ZERO
1028                   END DO
1029                   ISCOMPLEX( IV ) = IP
1030 *                 back-transform and normalization is done below
1031                END IF
1032             ELSE
1033 *
1034 *              --------------------------------------------------------
1035 *              Complex left eigenvector.
1036 *
1037 *              Initial solve:
1038 *              [ ( T(KI,KI)    T(KI,KI+1)  )**T - (WR - I* WI) ]*X = 0.
1039 *              [ ( T(KI+1,KI) T(KI+1,KI+1) )                   ]
1040 *
1041                IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN
1042                   WORK( KI   + (IV  )*N ) = WI / T( KI, KI+1 )
1043                   WORK( KI+1 + (IV+1)*N ) = ONE
1044                ELSE
1045                   WORK( KI   + (IV  )*N ) = ONE
1046                   WORK( KI+1 + (IV+1)*N ) = -WI / T( KI+1, KI )
1047                END IF
1048                WORK( KI+1 + (IV  )*N ) = ZERO
1049                WORK( KI   + (IV+1)*N ) = ZERO
1050 *
1051 *              Form right-hand side.
1052 *
1053                DO 190 K = KI + 2, N
1054                   WORK( K+(IV  )*N ) = -WORK( KI  +(IV  )*N )*T(KI,  K)
1055                   WORK( K+(IV+1)*N ) = -WORK( KI+1+(IV+1)*N )*T(KI+1,K)
1056   190          CONTINUE
1057 *
1058 *              Solve transposed quasi-triangular system:
1059 *              [ T(KI+2:N,KI+2:N)**T - (WR-i*WI) ]*X = WORK1+i*WORK2
1060 *
1061                VMAX = ONE
1062                VCRIT = BIGNUM
1063 *
1064                JNXT = KI + 2
1065                DO 200 J = KI + 2, N
1066                   IF( J.LT.JNXT )
1067      $               GO TO 200
1068                   J1 = J
1069                   J2 = J
1070                   JNXT = J + 1
1071                   IF( J.LT.N ) THEN
1072                      IF( T( J+1, J ).NE.ZERO ) THEN
1073                         J2 = J + 1
1074                         JNXT = J + 2
1075                      END IF
1076                   END IF
1077 *
1078                   IF( J1.EQ.J2 ) THEN
1079 *
1080 *                    1-by-1 diagonal block
1081 *
1082 *                    Scale if necessary to avoid overflow when
1083 *                    forming the right-hand side elements.
1084 *
1085                      IF( WORK( J ).GT.VCRIT ) THEN
1086                         REC = ONE / VMAX
1087                         CALL DSCAL( N-KI+1, REC, WORK(KI+(IV  )*N), 1 )
1088                         CALL DSCAL( N-KI+1, REC, WORK(KI+(IV+1)*N), 1 )
1089                         VMAX = ONE
1090                         VCRIT = BIGNUM
1091                      END IF
1092 *
1093                      WORK( J+(IV  )*N ) = WORK( J+(IV)*N ) -
1094      $                                  DDOT( J-KI-2, T( KI+2, J ), 1,
1095      $                                        WORK( KI+2+(IV)*N ), 1 )
1096                      WORK( J+(IV+1)*N ) = WORK( J+(IV+1)*N ) -
1097      $                                  DDOT( J-KI-2, T( KI+2, J ), 1,
1098      $                                        WORK( KI+2+(IV+1)*N ), 1 )
1099 *
1100 *                    Solve [ T(J,J)-(WR-i*WI) ]*(X11+i*X12)= WK+I*WK2
1101 *
1102                      CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
1103      $                            LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
1104      $                            -WI, X, 2, SCALE, XNORM, IERR )
1105 *
1106 *                    Scale if necessary
1107 *
1108                      IF( SCALE.NE.ONE ) THEN
1109                         CALL DSCAL( N-KI+1, SCALE, WORK(KI+(IV  )*N), 1)
1110                         CALL DSCAL( N-KI+1, SCALE, WORK(KI+(IV+1)*N), 1)
1111                      END IF
1112                      WORK( J+(IV  )*N ) = X( 1, 1 )
1113                      WORK( J+(IV+1)*N ) = X( 1, 2 )
1114                      VMAX = MAX( ABS( WORK( J+(IV  )*N ) ),
1115      $                           ABS( WORK( J+(IV+1)*N ) ), VMAX )
1116                      VCRIT = BIGNUM / VMAX
1117 *
1118                   ELSE
1119 *
1120 *                    2-by-2 diagonal block
1121 *
1122 *                    Scale if necessary to avoid overflow when forming
1123 *                    the right-hand side elements.
1124 *
1125                      BETA = MAX( WORK( J ), WORK( J+1 ) )
1126                      IF( BETA.GT.VCRIT ) THEN
1127                         REC = ONE / VMAX
1128                         CALL DSCAL( N-KI+1, REC, WORK(KI+(IV  )*N), 1 )
1129                         CALL DSCAL( N-KI+1, REC, WORK(KI+(IV+1)*N), 1 )
1130                         VMAX = ONE
1131                         VCRIT = BIGNUM
1132                      END IF
1133 *
1134                      WORK( J  +(IV  )*N ) = WORK( J+(IV)*N ) -
1135      $                                DDOT( J-KI-2, T( KI+2, J ), 1,
1136      $                                      WORK( KI+2+(IV)*N ), 1 )
1137 *
1138                      WORK( J  +(IV+1)*N ) = WORK( J+(IV+1)*N ) -
1139      $                                DDOT( J-KI-2, T( KI+2, J ), 1,
1140      $                                      WORK( KI+2+(IV+1)*N ), 1 )
1141 *
1142                      WORK( J+1+(IV  )*N ) = WORK( J+1+(IV)*N ) -
1143      $                                DDOT( J-KI-2, T( KI+2, J+1 ), 1,
1144      $                                      WORK( KI+2+(IV)*N ), 1 )
1145 *
1146                      WORK( J+1+(IV+1)*N ) = WORK( J+1+(IV+1)*N ) -
1147      $                                DDOT( J-KI-2, T( KI+2, J+1 ), 1,
1148      $                                      WORK( KI+2+(IV+1)*N ), 1 )
1149 *
1150 *                    Solve 2-by-2 complex linear equation
1151 *                    [ (T(j,j)   T(j,j+1)  )**T - (wr-i*wi)*I ]*X = SCALE*B
1152 *                    [ (T(j+1,j) T(j+1,j+1))                  ]
1153 *
1154                      CALL DLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ),
1155      $                            LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
1156      $                            -WI, X, 2, SCALE, XNORM, IERR )
1157 *
1158 *                    Scale if necessary
1159 *
1160                      IF( SCALE.NE.ONE ) THEN
1161                         CALL DSCAL( N-KI+1, SCALE, WORK(KI+(IV  )*N), 1)
1162                         CALL DSCAL( N-KI+1, SCALE, WORK(KI+(IV+1)*N), 1)
1163                      END IF
1164                      WORK( J  +(IV  )*N ) = X( 1, 1 )
1165                      WORK( J  +(IV+1)*N ) = X( 1, 2 )
1166                      WORK( J+1+(IV  )*N ) = X( 2, 1 )
1167                      WORK( J+1+(IV+1)*N ) = X( 2, 2 )
1168                      VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ),
1169      $                           ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ),
1170      $                           VMAX )
1171                      VCRIT = BIGNUM / VMAX
1172 *
1173                   END IF
1174   200          CONTINUE
1175 *
1176 *              Copy the vector x or Q*x to VL and normalize.
1177 *
1178                IF( .NOT.OVER ) THEN
1179 *                 ------------------------------
1180 *                 no back-transform: copy x to VL and normalize.
1181                   CALL DCOPY( N-KI+1, WORK( KI + (IV  )*N ), 1,
1182      $                        VL( KI, IS   ), 1 )
1183                   CALL DCOPY( N-KI+1, WORK( KI + (IV+1)*N ), 1,
1184      $                        VL( KI, IS+1 ), 1 )
1185 *
1186                   EMAX = ZERO
1187                   DO 220 K = KI, N
1188                      EMAX = MAX( EMAX, ABS( VL( K, IS   ) )+
1189      $                                 ABS( VL( K, IS+1 ) ) )
1190   220             CONTINUE
1191                   REMAX = ONE / EMAX
1192                   CALL DSCAL( N-KI+1, REMAX, VL( KI, IS   ), 1 )
1193                   CALL DSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 )
1194 *
1195                   DO 230 K = 1, KI - 1
1196                      VL( K, IS   ) = ZERO
1197                      VL( K, IS+1 ) = ZERO
1198   230             CONTINUE
1199 *
1200                ELSE IF( NB.EQ.1 ) THEN
1201 *                 ------------------------------
1202 *                 version 1: back-transform each vector with GEMV, Q*x.
1203                   IF( KI.LT.N-1 ) THEN
1204                      CALL DGEMV( 'N', N, N-KI-1, ONE,
1205      $                           VL( 1, KI+2 ), LDVL,
1206      $                           WORK( KI+2 + (IV)*N ), 1,
1207      $                           WORK( KI   + (IV)*N ),
1208      $                           VL( 1, KI ), 1 )
1209                      CALL DGEMV( 'N', N, N-KI-1, ONE,
1210      $                           VL( 1, KI+2 ), LDVL,
1211      $                           WORK( KI+2 + (IV+1)*N ), 1,
1212      $                           WORK( KI+1 + (IV+1)*N ),
1213      $                           VL( 1, KI+1 ), 1 )
1214                   ELSE
1215                      CALL DSCAL( N, WORK(KI+  (IV  )*N), VL(1, KI  ), 1)
1216                      CALL DSCAL( N, WORK(KI+1+(IV+1)*N), VL(1, KI+1), 1)
1217                   END IF
1218 *
1219                   EMAX = ZERO
1220                   DO 240 K = 1, N
1221                      EMAX = MAX( EMAX, ABS( VL( K, KI   ) )+
1222      $                                 ABS( VL( K, KI+1 ) ) )
1223   240             CONTINUE
1224                   REMAX = ONE / EMAX
1225                   CALL DSCAL( N, REMAX, VL( 1, KI   ), 1 )
1226                   CALL DSCAL( N, REMAX, VL( 1, KI+1 ), 1 )
1227 *
1228                ELSE
1229 *                 ------------------------------
1230 *                 version 2: back-transform block of vectors with GEMM
1231 *                 zero out above vector
1232 *                 could go from KI-NV+1 to KI-1
1233                   DO K = 1, KI - 1
1234                      WORK( K + (IV  )*N ) = ZERO
1235                      WORK( K + (IV+1)*N ) = ZERO
1236                   END DO
1237                   ISCOMPLEX( IV   ) =  IP
1238                   ISCOMPLEX( IV+1 ) = -IP
1239                   IV = IV + 1
1240 *                 back-transform and normalization is done below
1241                END IF
1242             END IF
1243
1244             IF( NB.GT.1 ) THEN
1245 *              --------------------------------------------------------
1246 *              Blocked version of back-transform
1247 *              For complex case, KI2 includes both vectors (KI and KI+1)
1248                IF( IP.EQ.0 ) THEN
1249                   KI2 = KI
1250                ELSE
1251                   KI2 = KI + 1
1252                END IF
1253
1254 *              Columns 1:IV of work are valid vectors.
1255 *              When the number of vectors stored reaches NB-1 or NB,
1256 *              or if this was last vector, do the GEMM
1257                IF( (IV.GE.NB-1) .OR. (KI2.EQ.N) ) THEN
1258                   CALL DGEMM( 'N', 'N', N, IV, N-KI2+IV, ONE,
1259      $                        VL( 1, KI2-IV+1 ), LDVL,
1260      $                        WORK( KI2-IV+1 + (1)*N ), N,
1261      $                        ZERO,
1262      $                        WORK( 1 + (NB+1)*N ), N )
1263 *                 normalize vectors
1264                   DO K = 1, IV
1265                      IF( ISCOMPLEX(K).EQ.0) THEN
1266 *                       real eigenvector
1267                         II = IDAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
1268                         REMAX = ONE / ABS( WORK( II + (NB+K)*N ) )
1269                      ELSE IF( ISCOMPLEX(K).EQ.1) THEN
1270 *                       first eigenvector of conjugate pair
1271                         EMAX = ZERO
1272                         DO II = 1, N
1273                            EMAX = MAX( EMAX,
1274      $                                 ABS( WORK( II + (NB+K  )*N ) )+
1275      $                                 ABS( WORK( II + (NB+K+1)*N ) ) )
1276                         END DO
1277                         REMAX = ONE / EMAX
1278 *                    else if ISCOMPLEX(K).EQ.-1
1279 *                       second eigenvector of conjugate pair
1280 *                       reuse same REMAX as previous K
1281                      END IF
1282                      CALL DSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
1283                   END DO
1284                   CALL DLACPY( 'F', N, IV,
1285      $                         WORK( 1 + (NB+1)*N ), N,
1286      $                         VL( 1, KI2-IV+1 ), LDVL )
1287                   IV = 1
1288                ELSE
1289                   IV = IV + 1
1290                END IF
1291             END IF ! blocked back-transform
1292 *
1293             IS = IS + 1
1294             IF( IP.NE.0 )
1295      $         IS = IS + 1
1296   260    CONTINUE
1297       END IF
1298 *
1299       RETURN
1300 *
1301 *     End of DTREVC3
1302 *
1303       END