SRC/ztrevc3.f

   1 *> \brief \b ZTREVC3
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZTREVC3 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrevc3.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrevc3.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrevc3.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
  22 *      $                    LDVR, MM, M, WORK, LWORK, RWORK, LRWORK, 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   RWORK( * )
  31 *       COMPLEX*16         T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
  32 *      $                   WORK( * )
  33 *       ..
  34 *
  35 *
  36 *> \par Purpose:
  37 *  =============
  38 *>
  39 *> \verbatim
  40 *>
  41 *> ZTREVC3 computes some or all of the right and/or left eigenvectors of
  42 *> a complex upper triangular matrix T.
  43 *> Matrices of this type are produced by the Schur factorization of
  44 *> a complex general matrix:  A = Q*T*Q**H, as computed by ZHSEQR.
  45 *>
  46 *> The right eigenvector x and the left eigenvector y of T corresponding
  47 *> to an eigenvalue w are defined by:
  48 *>
  49 *>              T*x = w*x,     (y**H)*T = w*(y**H)
  50 *>
  51 *> where y**H denotes the conjugate transpose of the vector y.
  52 *> The eigenvalues are not input to this routine, but are read directly
  53 *> from the diagonal of T.
  54 *>
  55 *> This routine returns the matrices X and/or Y of right and left
  56 *> eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
  57 *> input matrix. If Q is the unitary factor that reduces a matrix A to
  58 *> Schur form T, then Q*X and Q*Y are the matrices of right and left
  59 *> eigenvectors of A.
  60 *>
  61 *> This uses a Level 3 BLAS version of the back transformation.
  62 *> \endverbatim
  63 *
  64 *  Arguments:
  65 *  ==========
  66 *
  67 *> \param[in] SIDE
  68 *> \verbatim
  69 *>          SIDE is CHARACTER*1
  70 *>          = 'R':  compute right eigenvectors only;
  71 *>          = 'L':  compute left eigenvectors only;
  72 *>          = 'B':  compute both right and left eigenvectors.
  73 *> \endverbatim
  74 *>
  75 *> \param[in] HOWMNY
  76 *> \verbatim
  77 *>          HOWMNY is CHARACTER*1
  78 *>          = 'A':  compute all right and/or left eigenvectors;
  79 *>          = 'B':  compute all right and/or left eigenvectors,
  80 *>                  backtransformed using the matrices supplied in
  81 *>                  VR and/or VL;
  82 *>          = 'S':  compute selected right and/or left eigenvectors,
  83 *>                  as indicated by the logical array SELECT.
  84 *> \endverbatim
  85 *>
  86 *> \param[in] SELECT
  87 *> \verbatim
  88 *>          SELECT is LOGICAL array, dimension (N)
  89 *>          If HOWMNY = 'S', SELECT specifies the eigenvectors to be
  90 *>          computed.
  91 *>          The eigenvector corresponding to the j-th eigenvalue is
  92 *>          computed if SELECT(j) = .TRUE..
  93 *>          Not referenced if HOWMNY = 'A' or 'B'.
  94 *> \endverbatim
  95 *>
  96 *> \param[in] N
  97 *> \verbatim
  98 *>          N is INTEGER
  99 *>          The order of the matrix T. N >= 0.
 100 *> \endverbatim
 101 *>
 102 *> \param[in,out] T
 103 *> \verbatim
 104 *>          T is COMPLEX*16 array, dimension (LDT,N)
 105 *>          The upper triangular matrix T.  T is modified, but restored
 106 *>          on exit.
 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 COMPLEX*16 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 unitary matrix Q of
 120 *>          Schur vectors returned by ZHSEQR).
 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 *>          Not referenced if SIDE = 'R'.
 129 *> \endverbatim
 130 *>
 131 *> \param[in] LDVL
 132 *> \verbatim
 133 *>          LDVL is INTEGER
 134 *>          The leading dimension of the array VL.
 135 *>          LDVL >= 1, and if SIDE = 'L' or 'B', LDVL >= N.
 136 *> \endverbatim
 137 *>
 138 *> \param[in,out] VR
 139 *> \verbatim
 140 *>          VR is COMPLEX*16 array, dimension (LDVR,MM)
 141 *>          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
 142 *>          contain an N-by-N matrix Q (usually the unitary matrix Q of
 143 *>          Schur vectors returned by ZHSEQR).
 144 *>          On exit, if SIDE = 'R' or 'B', VR contains:
 145 *>          if HOWMNY = 'A', the matrix X of right eigenvectors of T;
 146 *>          if HOWMNY = 'B', the matrix Q*X;
 147 *>          if HOWMNY = 'S', the right eigenvectors of T specified by
 148 *>                           SELECT, stored consecutively in the columns
 149 *>                           of VR, in the same order as their
 150 *>                           eigenvalues.
 151 *>          Not referenced if SIDE = 'L'.
 152 *> \endverbatim
 153 *>
 154 *> \param[in] LDVR
 155 *> \verbatim
 156 *>          LDVR is INTEGER
 157 *>          The leading dimension of the array VR.
 158 *>          LDVR >= 1, and if SIDE = 'R' or 'B', LDVR >= N.
 159 *> \endverbatim
 160 *>
 161 *> \param[in] MM
 162 *> \verbatim
 163 *>          MM is INTEGER
 164 *>          The number of columns in the arrays VL and/or VR. MM >= M.
 165 *> \endverbatim
 166 *>
 167 *> \param[out] M
 168 *> \verbatim
 169 *>          M is INTEGER
 170 *>          The number of columns in the arrays VL and/or VR actually
 171 *>          used to store the eigenvectors.
 172 *>          If HOWMNY = 'A' or 'B', M is set to N.
 173 *>          Each selected eigenvector occupies one column.
 174 *> \endverbatim
 175 *>
 176 *> \param[out] WORK
 177 *> \verbatim
 178 *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
 179 *> \endverbatim
 180 *>
 181 *> \param[in] LWORK
 182 *> \verbatim
 183 *>          LWORK is INTEGER
 184 *>          The dimension of array WORK. LWORK >= max(1,2*N).
 185 *>          For optimum performance, LWORK >= N + 2*N*NB, where NB is
 186 *>          the optimal blocksize.
 187 *>
 188 *>          If LWORK = -1, then a workspace query is assumed; the routine
 189 *>          only calculates the optimal size of the WORK array, returns
 190 *>          this value as the first entry of the WORK array, and no error
 191 *>          message related to LWORK is issued by XERBLA.
 192 *> \endverbatim
 193 *>
 194 *> \param[out] RWORK
 195 *> \verbatim
 196 *>          RWORK is DOUBLE PRECISION array, dimension (LRWORK)
 197 *> \endverbatim
 198 *>
 199 *> \param[in] LRWORK
 200 *> \verbatim
 201 *>          LRWORK is INTEGER
 202 *>          The dimension of array RWORK. LRWORK >= max(1,N).
 203 *>
 204 *>          If LRWORK = -1, then a workspace query is assumed; the routine
 205 *>          only calculates the optimal size of the RWORK array, returns
 206 *>          this value as the first entry of the RWORK array, and no error
 207 *>          message related to LRWORK is issued by XERBLA.
 208 *> \endverbatim
 209 *>
 210 *> \param[out] INFO
 211 *> \verbatim
 212 *>          INFO is INTEGER
 213 *>          = 0:  successful exit
 214 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 215 *> \endverbatim
 216 *
 217 *  Authors:
 218 *  ========
 219 *
 220 *> \author Univ. of Tennessee
 221 *> \author Univ. of California Berkeley
 222 *> \author Univ. of Colorado Denver
 223 *> \author NAG Ltd.
 224 *
 225 *> \date November 2011
 226 *
 227 *  @precisions fortran z -> c
 228 *
 229 *> \ingroup complex16OTHERcomputational
 230 *
 231 *> \par Further Details:
 232 *  =====================
 233 *>
 234 *> \verbatim
 235 *>
 236 *>  The algorithm used in this program is basically backward (forward)
 237 *>  substitution, with scaling to make the the code robust against
 238 *>  possible overflow.
 239 *>
 240 *>  Each eigenvector is normalized so that the element of largest
 241 *>  magnitude has magnitude 1; here the magnitude of a complex number
 242 *>  (x,y) is taken to be |x| + |y|.
 243 *> \endverbatim
 244 *>
 245 *  =====================================================================
 246       SUBROUTINE ZTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
 247      $                    LDVR, MM, M, WORK, LWORK, RWORK, LRWORK, INFO)
 248       IMPLICIT NONE
 249 *
 250 *  -- LAPACK computational routine (version 3.4.0) --
 251 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 252 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 253 *     November 2011
 254 *
 255 *     .. Scalar Arguments ..
 256       CHARACTER          HOWMNY, SIDE
 257       INTEGER            INFO, LDT, LDVL, LDVR, LWORK, LRWORK, M, MM, N
 258 *     ..
 259 *     .. Array Arguments ..
 260       LOGICAL            SELECT( * )
 261       DOUBLE PRECISION   RWORK( * )
 262       COMPLEX*16         T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
 263      $                   WORK( * )
 264 *     ..
 265 *
 266 *  =====================================================================
 267 *
 268 *     .. Parameters ..
 269       DOUBLE PRECISION   ZERO, ONE
 270       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 271       COMPLEX*16         CZERO, CONE
 272       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
 273      $                     CONE  = ( 1.0D+0, 0.0D+0 ) )
 274       INTEGER            NBMIN, NBMAX
 275       PARAMETER          ( NBMIN = 8, NBMAX = 128 )
 276 *     ..
 277 *     .. Local Scalars ..
 278       LOGICAL            ALLV, BOTHV, LEFTV, LQUERY, OVER, RIGHTV, SOMEV
 279       INTEGER            I, II, IS, J, K, KI, IV, MAXWRK, NB
 280       DOUBLE PRECISION   OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL
 281       COMPLEX*16         CDUM
 282 *     ..
 283 *     .. External Functions ..
 284       LOGICAL            LSAME
 285       INTEGER            ILAENV, IZAMAX
 286       DOUBLE PRECISION   DLAMCH, DZASUM
 287       EXTERNAL           LSAME, ILAENV, IZAMAX, DLAMCH, DZASUM
 288 *     ..
 289 *     .. External Subroutines ..
 290       EXTERNAL           XERBLA, ZCOPY, ZDSCAL, ZGEMV, ZLATRS
 291 *     ..
 292 *     .. Intrinsic Functions ..
 293       INTRINSIC          ABS, DBLE, DCMPLX, CONJG, AIMAG, MAX
 294 *     ..
 295 *     .. Statement Functions ..
 296       DOUBLE PRECISION   CABS1
 297 *     ..
 298 *     .. Statement Function definitions ..
 299       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( AIMAG( CDUM ) )
 300 *     ..
 301 *     .. Executable Statements ..
 302 *
 303 *     Decode and test the input parameters
 304 *
 305       BOTHV  = LSAME( SIDE, 'B' )
 306       RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
 307       LEFTV  = LSAME( SIDE, 'L' ) .OR. BOTHV
 308 *
 309       ALLV  = LSAME( HOWMNY, 'A' )
 310       OVER  = LSAME( HOWMNY, 'B' )
 311       SOMEV = LSAME( HOWMNY, 'S' )
 312 *
 313 *     Set M to the number of columns required to store the selected
 314 *     eigenvectors.
 315 *
 316       IF( SOMEV ) THEN
 317          M = 0
 318          DO 10 J = 1, N
 319             IF( SELECT( J ) )
 320      $         M = M + 1
 321    10    CONTINUE
 322       ELSE
 323          M = N
 324       END IF
 325 *
 326       INFO = 0
 327       NB = ILAENV( 1, 'ZTREVC', SIDE // HOWMNY, N, -1, -1, -1 )
 328       MAXWRK = N + 2*N*NB
 329       WORK(1) = MAXWRK
 330       RWORK(1) = N
 331       LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 )
 332       IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
 333          INFO = -1
 334       ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
 335          INFO = -2
 336       ELSE IF( N.LT.0 ) THEN
 337          INFO = -4
 338       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
 339          INFO = -6
 340       ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
 341          INFO = -8
 342       ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
 343          INFO = -10
 344       ELSE IF( MM.LT.M ) THEN
 345          INFO = -11
 346       ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN
 347          INFO = -14
 348       ELSE IF ( LRWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
 349          INFO = -16
 350       END IF
 351       IF( INFO.NE.0 ) THEN
 352          CALL XERBLA( 'ZTREVC3', -INFO )
 353          RETURN
 354       ELSE IF( LQUERY ) THEN
 355          RETURN
 356       END IF
 357 *
 358 *     Quick return if possible.
 359 *
 360       IF( N.EQ.0 )
 361      $   RETURN
 362 *
 363 *     Use blocked version of back-transformation if sufficient workspace.
 364 *     Zero-out the workspace to avoid potential NaN propagation.
 365 *
 366       IF( OVER .AND. LWORK .GE. N + 2*N*NBMIN ) THEN
 367          NB = (LWORK - N) / (2*N)
 368          NB = MIN( NB, NBMAX )
 369          CALL ZLASET( 'F', N, 1+2*NB, CZERO, CZERO, WORK, N )
 370       ELSE
 371          NB = 1
 372       END IF
 373 *
 374 *     Set the constants to control overflow.
 375 *
 376       UNFL = DLAMCH( 'Safe minimum' )
 377       OVFL = ONE / UNFL
 378       CALL DLABAD( UNFL, OVFL )
 379       ULP = DLAMCH( 'Precision' )
 380       SMLNUM = UNFL*( N / ULP )
 381 *
 382 *     Store the diagonal elements of T in working array WORK.
 383 *
 384       DO 20 I = 1, N
 385          WORK( I ) = T( I, I )
 386    20 CONTINUE
 387 *
 388 *     Compute 1-norm of each column of strictly upper triangular
 389 *     part of T to control overflow in triangular solver.
 390 *
 391       RWORK( 1 ) = ZERO
 392       DO 30 J = 2, N
 393          RWORK( J ) = DZASUM( J-1, T( 1, J ), 1 )
 394    30 CONTINUE
 395 *
 396       IF( RIGHTV ) THEN
 397 *
 398 *        ============================================================
 399 *        Compute right eigenvectors.
 400 *
 401 *        IV is index of column in current block.
 402 *        Non-blocked version always uses IV=NB=1;
 403 *        blocked     version starts with IV=NB, goes down to 1.
 404 *        (Note the "0-th" column is used to store the original diagonal.)
 405          IV = NB
 406          IS = M
 407          DO 80 KI = N, 1, -1
 408             IF( SOMEV ) THEN
 409                IF( .NOT.SELECT( KI ) )
 410      $            GO TO 80
 411             END IF
 412             SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
 413 *
 414 *           --------------------------------------------------------
 415 *           Complex right eigenvector
 416 *
 417             WORK( KI + IV*N ) = CONE
 418 *
 419 *           Form right-hand side.
 420 *
 421             DO 40 K = 1, KI - 1
 422                WORK( K + IV*N ) = -T( K, KI )
 423    40       CONTINUE
 424 *
 425 *           Solve upper triangular system:
 426 *           [ T(1:KI-1,1:KI-1) - T(KI,KI) ]*X = SCALE*WORK.
 427 *
 428             DO 50 K = 1, KI - 1
 429                T( K, K ) = T( K, K ) - T( KI, KI )
 430                IF( CABS1( T( K, K ) ).LT.SMIN )
 431      $            T( K, K ) = SMIN
 432    50       CONTINUE
 433 *
 434             IF( KI.GT.1 ) THEN
 435                CALL ZLATRS( 'Upper', 'No transpose', 'Non-unit', 'Y',
 436      $                      KI-1, T, LDT, WORK( 1 + IV*N ), SCALE,
 437      $                      RWORK, INFO )
 438                WORK( KI + IV*N ) = SCALE
 439             END IF
 440 *
 441 *           Copy the vector x or Q*x to VR and normalize.
 442 *
 443             IF( .NOT.OVER ) THEN
 444 *              ------------------------------
 445 *              no back-transform: copy x to VR and normalize.
 446                CALL ZCOPY( KI, WORK( 1 + IV*N ), 1, VR( 1, IS ), 1 )
 447 *
 448                II = IZAMAX( KI, VR( 1, IS ), 1 )
 449                REMAX = ONE / CABS1( VR( II, IS ) )
 450                CALL ZDSCAL( KI, REMAX, VR( 1, IS ), 1 )
 451 *
 452                DO 60 K = KI + 1, N
 453                   VR( K, IS ) = CZERO
 454    60          CONTINUE
 455 *
 456             ELSE IF( NB.EQ.1 ) THEN
 457 *              ------------------------------
 458 *              version 1: back-transform each vector with GEMV, Q*x.
 459                IF( KI.GT.1 )
 460      $            CALL ZGEMV( 'N', N, KI-1, CONE, VR, LDVR,
 461      $                        WORK( 1 + IV*N ), 1, DCMPLX( SCALE ),
 462      $                        VR( 1, KI ), 1 )
 463 *
 464                II = IZAMAX( N, VR( 1, KI ), 1 )
 465                REMAX = ONE / CABS1( VR( II, KI ) )
 466                CALL ZDSCAL( N, REMAX, VR( 1, KI ), 1 )
 467 *
 468             ELSE
 469 *              ------------------------------
 470 *              version 2: back-transform block of vectors with GEMM
 471 *              zero out below vector
 472                DO K = KI + 1, N
 473                   WORK( K + IV*N ) = CZERO
 474                END DO
 475 *
 476 *              Columns IV:NB of work are valid vectors.
 477 *              When the number of vectors stored reaches NB,
 478 *              or if this was last vector, do the GEMM
 479                IF( (IV.EQ.1) .OR. (KI.EQ.1) ) THEN
 480                   CALL ZGEMM( 'N', 'N', N, NB-IV+1, KI+NB-IV, CONE,
 481      $                        VR, LDVR,
 482      $                        WORK( 1 + (IV)*N    ), N,
 483      $                        CZERO,
 484      $                        WORK( 1 + (NB+IV)*N ), N )
 485 *                 normalize vectors
 486                   DO K = IV, NB
 487                      II = IZAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
 488                      REMAX = ONE / CABS1( WORK( II + (NB+K)*N ) )
 489                      CALL ZDSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
 490                   END DO
 491                   CALL ZLACPY( 'F', N, NB-IV+1,
 492      $                         WORK( 1 + (NB+IV)*N ), N,
 493      $                         VR( 1, KI ), LDVR )
 494                   IV = NB
 495                ELSE
 496                   IV = IV - 1
 497                END IF
 498             END IF
 499 *
 500 *           Restore the original diagonal elements of T.
 501 *
 502             DO 70 K = 1, KI - 1
 503                T( K, K ) = WORK( K )
 504    70       CONTINUE
 505 *
 506             IS = IS - 1
 507    80    CONTINUE
 508       END IF
 509 *
 510       IF( LEFTV ) THEN
 511 *
 512 *        ============================================================
 513 *        Compute left eigenvectors.
 514 *
 515 *        IV is index of column in current block.
 516 *        Non-blocked version always uses IV=1;
 517 *        blocked     version starts with IV=1, goes up to NB.
 518 *        (Note the "0-th" column is used to store the original diagonal.)
 519          IV = 1
 520          IS = 1
 521          DO 130 KI = 1, N
 522 *
 523             IF( SOMEV ) THEN
 524                IF( .NOT.SELECT( KI ) )
 525      $            GO TO 130
 526             END IF
 527             SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
 528 *
 529 *           --------------------------------------------------------
 530 *           Complex left eigenvector
 531 *
 532             WORK( KI + IV*N ) = CONE
 533 *
 534 *           Form right-hand side.
 535 *
 536             DO 90 K = KI + 1, N
 537                WORK( K + IV*N ) = -CONJG( T( KI, K ) )
 538    90       CONTINUE
 539 *
 540 *           Solve conjugate-transposed triangular system:
 541 *           [ T(KI+1:N,KI+1:N) - T(KI,KI) ]**H * X = SCALE*WORK.
 542 *
 543             DO 100 K = KI + 1, N
 544                T( K, K ) = T( K, K ) - T( KI, KI )
 545                IF( CABS1( T( K, K ) ).LT.SMIN )
 546      $            T( K, K ) = SMIN
 547   100       CONTINUE
 548 *
 549             IF( KI.LT.N ) THEN
 550                CALL ZLATRS( 'Upper', 'Conjugate transpose', 'Non-unit',
 551      $                      'Y', N-KI, T( KI+1, KI+1 ), LDT,
 552      $                      WORK( KI+1 + IV*N ), SCALE, RWORK, INFO )
 553                WORK( KI + IV*N ) = SCALE
 554             END IF
 555 *
 556 *           Copy the vector x or Q*x to VL and normalize.
 557 *
 558             IF( .NOT.OVER ) THEN
 559 *              ------------------------------
 560 *              no back-transform: copy x to VL and normalize.
 561                CALL ZCOPY( N-KI+1, WORK( KI + IV*N ), 1, VL(KI,IS), 1 )
 562 *
 563                II = IZAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
 564                REMAX = ONE / CABS1( VL( II, IS ) )
 565                CALL ZDSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
 566 *
 567                DO 110 K = 1, KI - 1
 568                   VL( K, IS ) = CZERO
 569   110          CONTINUE
 570 *
 571             ELSE IF( NB.EQ.1 ) THEN
 572 *              ------------------------------
 573 *              version 1: back-transform each vector with GEMV, Q*x.
 574                IF( KI.LT.N )
 575      $            CALL ZGEMV( 'N', N, N-KI, CONE, VL( 1, KI+1 ), LDVL,
 576      $                        WORK( KI+1 + IV*N ), 1, DCMPLX( SCALE ),
 577      $                        VL( 1, KI ), 1 )
 578 *
 579                II = IZAMAX( N, VL( 1, KI ), 1 )
 580                REMAX = ONE / CABS1( VL( II, KI ) )
 581                CALL ZDSCAL( N, REMAX, VL( 1, KI ), 1 )
 582 *
 583             ELSE
 584 *              ------------------------------
 585 *              version 2: back-transform block of vectors with GEMM
 586 *              zero out above vector
 587 *              could go from KI-NV+1 to KI-1
 588                DO K = 1, KI - 1
 589                   WORK( K + IV*N ) = CZERO
 590                END DO
 591 *
 592 *              Columns 1:IV of work are valid vectors.
 593 *              When the number of vectors stored reaches NB,
 594 *              or if this was last vector, do the GEMM
 595                IF( (IV.EQ.NB) .OR. (KI.EQ.N) ) THEN
 596                   CALL ZGEMM( 'N', 'N', N, IV, N-KI+IV, CONE,
 597      $                        VL( 1, KI-IV+1 ), LDVL,
 598      $                        WORK( KI-IV+1 + (1)*N ), N,
 599      $                        CZERO,
 600      $                        WORK( 1 + (NB+1)*N ), N )
 601 *                 normalize vectors
 602                   DO K = 1, IV
 603                      II = IZAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
 604                      REMAX = ONE / CABS1( WORK( II + (NB+K)*N ) )
 605                      CALL ZDSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
 606                   END DO
 607                   CALL ZLACPY( 'F', N, IV,
 608      $                         WORK( 1 + (NB+1)*N ), N,
 609      $                         VL( 1, KI-IV+1 ), LDVL )
 610                   IV = 1
 611                ELSE
 612                   IV = IV + 1
 613                END IF
 614             END IF
 615 *
 616 *           Restore the original diagonal elements of T.
 617 *
 618             DO 120 K = KI + 1, N
 619                T( K, K ) = WORK( K )
 620   120       CONTINUE
 621 *
 622             IS = IS + 1
 623   130    CONTINUE
 624       END IF
 625 *
 626       RETURN
 627 *
 628 *     End of ZTREVC3
 629 *
 630       END