SRC/ctgevc.f

   1 *> \brief \b CTGEVC
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CTGEVC + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctgevc.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctgevc.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctgevc.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
  22 *                          LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          HOWMNY, SIDE
  26 *       INTEGER            INFO, LDP, LDS, LDVL, LDVR, M, MM, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       LOGICAL            SELECT( * )
  30 *       REAL               RWORK( * )
  31 *       COMPLEX            P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
  32 *      $                   VR( LDVR, * ), WORK( * )
  33 *       ..
  34 *
  35 *
  36 *
  37 *> \par Purpose:
  38 *  =============
  39 *>
  40 *> \verbatim
  41 *>
  42 *> CTGEVC computes some or all of the right and/or left eigenvectors of
  43 *> a pair of complex matrices (S,P), where S and P are upper triangular.
  44 *> Matrix pairs of this type are produced by the generalized Schur
  45 *> factorization of a complex matrix pair (A,B):
  46 *>
  47 *>    A = Q*S*Z**H,  B = Q*P*Z**H
  48 *>
  49 *> as computed by CGGHRD + CHGEQZ.
  50 *>
  51 *> The right eigenvector x and the left eigenvector y of (S,P)
  52 *> corresponding to an eigenvalue w are defined by:
  53 *>
  54 *>    S*x = w*P*x,  (y**H)*S = w*(y**H)*P,
  55 *>
  56 *> where y**H denotes the conjugate tranpose of y.
  57 *> The eigenvalues are not input to this routine, but are computed
  58 *> directly from the diagonal elements of S and P.
  59 *>
  60 *> This routine returns the matrices X and/or Y of right and left
  61 *> eigenvectors of (S,P), or the products Z*X and/or Q*Y,
  62 *> where Z and Q are input matrices.
  63 *> If Q and Z are the unitary factors from the generalized Schur
  64 *> factorization of a matrix pair (A,B), then Z*X and Q*Y
  65 *> are the matrices of right and left eigenvectors of (A,B).
  66 *> \endverbatim
  67 *
  68 *  Arguments:
  69 *  ==========
  70 *
  71 *> \param[in] SIDE
  72 *> \verbatim
  73 *>          SIDE is CHARACTER*1
  74 *>          = 'R': compute right eigenvectors only;
  75 *>          = 'L': compute left eigenvectors only;
  76 *>          = 'B': compute both right and left eigenvectors.
  77 *> \endverbatim
  78 *>
  79 *> \param[in] HOWMNY
  80 *> \verbatim
  81 *>          HOWMNY is CHARACTER*1
  82 *>          = 'A': compute all right and/or left eigenvectors;
  83 *>          = 'B': compute all right and/or left eigenvectors,
  84 *>                 backtransformed by the matrices in VR and/or VL;
  85 *>          = 'S': compute selected right and/or left eigenvectors,
  86 *>                 specified by the logical array SELECT.
  87 *> \endverbatim
  88 *>
  89 *> \param[in] SELECT
  90 *> \verbatim
  91 *>          SELECT is LOGICAL array, dimension (N)
  92 *>          If HOWMNY='S', SELECT specifies the eigenvectors to be
  93 *>          computed.  The eigenvector corresponding to the j-th
  94 *>          eigenvalue is computed if SELECT(j) = .TRUE..
  95 *>          Not referenced if HOWMNY = 'A' or 'B'.
  96 *> \endverbatim
  97 *>
  98 *> \param[in] N
  99 *> \verbatim
 100 *>          N is INTEGER
 101 *>          The order of the matrices S and P.  N >= 0.
 102 *> \endverbatim
 103 *>
 104 *> \param[in] S
 105 *> \verbatim
 106 *>          S is COMPLEX array, dimension (LDS,N)
 107 *>          The upper triangular matrix S from a generalized Schur
 108 *>          factorization, as computed by CHGEQZ.
 109 *> \endverbatim
 110 *>
 111 *> \param[in] LDS
 112 *> \verbatim
 113 *>          LDS is INTEGER
 114 *>          The leading dimension of array S.  LDS >= max(1,N).
 115 *> \endverbatim
 116 *>
 117 *> \param[in] P
 118 *> \verbatim
 119 *>          P is COMPLEX array, dimension (LDP,N)
 120 *>          The upper triangular matrix P from a generalized Schur
 121 *>          factorization, as computed by CHGEQZ.  P must have real
 122 *>          diagonal elements.
 123 *> \endverbatim
 124 *>
 125 *> \param[in] LDP
 126 *> \verbatim
 127 *>          LDP is INTEGER
 128 *>          The leading dimension of array P.  LDP >= max(1,N).
 129 *> \endverbatim
 130 *>
 131 *> \param[in,out] VL
 132 *> \verbatim
 133 *>          VL is COMPLEX array, dimension (LDVL,MM)
 134 *>          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
 135 *>          contain an N-by-N matrix Q (usually the unitary matrix Q
 136 *>          of left Schur vectors returned by CHGEQZ).
 137 *>          On exit, if SIDE = 'L' or 'B', VL contains:
 138 *>          if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
 139 *>          if HOWMNY = 'B', the matrix Q*Y;
 140 *>          if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
 141 *>                      SELECT, stored consecutively in the columns of
 142 *>                      VL, in the same order as their eigenvalues.
 143 *>          Not referenced if SIDE = 'R'.
 144 *> \endverbatim
 145 *>
 146 *> \param[in] LDVL
 147 *> \verbatim
 148 *>          LDVL is INTEGER
 149 *>          The leading dimension of array VL.  LDVL >= 1, and if
 150 *>          SIDE = 'L' or 'l' or 'B' or 'b', LDVL >= N.
 151 *> \endverbatim
 152 *>
 153 *> \param[in,out] VR
 154 *> \verbatim
 155 *>          VR is COMPLEX array, dimension (LDVR,MM)
 156 *>          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
 157 *>          contain an N-by-N matrix Q (usually the unitary matrix Z
 158 *>          of right Schur vectors returned by CHGEQZ).
 159 *>          On exit, if SIDE = 'R' or 'B', VR contains:
 160 *>          if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
 161 *>          if HOWMNY = 'B', the matrix Z*X;
 162 *>          if HOWMNY = 'S', the right eigenvectors of (S,P) specified by
 163 *>                      SELECT, stored consecutively in the columns of
 164 *>                      VR, in the same order as their eigenvalues.
 165 *>          Not referenced if SIDE = 'L'.
 166 *> \endverbatim
 167 *>
 168 *> \param[in] LDVR
 169 *> \verbatim
 170 *>          LDVR is INTEGER
 171 *>          The leading dimension of the array VR.  LDVR >= 1, and if
 172 *>          SIDE = 'R' or 'B', LDVR >= N.
 173 *> \endverbatim
 174 *>
 175 *> \param[in] MM
 176 *> \verbatim
 177 *>          MM is INTEGER
 178 *>          The number of columns in the arrays VL and/or VR. MM >= M.
 179 *> \endverbatim
 180 *>
 181 *> \param[out] M
 182 *> \verbatim
 183 *>          M is INTEGER
 184 *>          The number of columns in the arrays VL and/or VR actually
 185 *>          used to store the eigenvectors.  If HOWMNY = 'A' or 'B', M
 186 *>          is set to N.  Each selected eigenvector occupies one column.
 187 *> \endverbatim
 188 *>
 189 *> \param[out] WORK
 190 *> \verbatim
 191 *>          WORK is COMPLEX array, dimension (2*N)
 192 *> \endverbatim
 193 *>
 194 *> \param[out] RWORK
 195 *> \verbatim
 196 *>          RWORK is REAL array, dimension (2*N)
 197 *> \endverbatim
 198 *>
 199 *> \param[out] INFO
 200 *> \verbatim
 201 *>          INFO is INTEGER
 202 *>          = 0:  successful exit.
 203 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 204 *> \endverbatim
 205 *
 206 *  Authors:
 207 *  ========
 208 *
 209 *> \author Univ. of Tennessee
 210 *> \author Univ. of California Berkeley
 211 *> \author Univ. of Colorado Denver
 212 *> \author NAG Ltd.
 213 *
 214 *> \date November 2011
 215 *
 216 *> \ingroup complexGEcomputational
 217 *
 218 *  =====================================================================
 219       SUBROUTINE CTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
 220      $                   LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
 221 *
 222 *  -- LAPACK computational routine (version 3.4.0) --
 223 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 224 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 225 *     November 2011
 226 *
 227 *     .. Scalar Arguments ..
 228       CHARACTER          HOWMNY, SIDE
 229       INTEGER            INFO, LDP, LDS, LDVL, LDVR, M, MM, N
 230 *     ..
 231 *     .. Array Arguments ..
 232       LOGICAL            SELECT( * )
 233       REAL               RWORK( * )
 234       COMPLEX            P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
 235      $                   VR( LDVR, * ), WORK( * )
 236 *     ..
 237 *
 238 *
 239 *  =====================================================================
 240 *
 241 *     .. Parameters ..
 242       REAL               ZERO, ONE
 243       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 244       COMPLEX            CZERO, CONE
 245       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
 246      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
 247 *     ..
 248 *     .. Local Scalars ..
 249       LOGICAL            COMPL, COMPR, ILALL, ILBACK, ILBBAD, ILCOMP,
 250      $                   LSA, LSB
 251       INTEGER            I, IBEG, IEIG, IEND, IHWMNY, IM, ISIDE, ISRC,
 252      $                   J, JE, JR
 253       REAL               ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG,
 254      $                   BIGNUM, BNORM, BSCALE, DMIN, SAFMIN, SBETA,
 255      $                   SCALE, SMALL, TEMP, ULP, XMAX
 256       COMPLEX            BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X
 257 *     ..
 258 *     .. External Functions ..
 259       LOGICAL            LSAME
 260       REAL               SLAMCH
 261       COMPLEX            CLADIV
 262       EXTERNAL           LSAME, SLAMCH, CLADIV
 263 *     ..
 264 *     .. External Subroutines ..
 265       EXTERNAL           CGEMV, SLABAD, XERBLA
 266 *     ..
 267 *     .. Intrinsic Functions ..
 268       INTRINSIC          ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL
 269 *     ..
 270 *     .. Statement Functions ..
 271       REAL               ABS1
 272 *     ..
 273 *     .. Statement Function definitions ..
 274       ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
 275 *     ..
 276 *     .. Executable Statements ..
 277 *
 278 *     Decode and Test the input parameters
 279 *
 280       IF( LSAME( HOWMNY, 'A' ) ) THEN
 281          IHWMNY = 1
 282          ILALL = .TRUE.
 283          ILBACK = .FALSE.
 284       ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN
 285          IHWMNY = 2
 286          ILALL = .FALSE.
 287          ILBACK = .FALSE.
 288       ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN
 289          IHWMNY = 3
 290          ILALL = .TRUE.
 291          ILBACK = .TRUE.
 292       ELSE
 293          IHWMNY = -1
 294       END IF
 295 *
 296       IF( LSAME( SIDE, 'R' ) ) THEN
 297          ISIDE = 1
 298          COMPL = .FALSE.
 299          COMPR = .TRUE.
 300       ELSE IF( LSAME( SIDE, 'L' ) ) THEN
 301          ISIDE = 2
 302          COMPL = .TRUE.
 303          COMPR = .FALSE.
 304       ELSE IF( LSAME( SIDE, 'B' ) ) THEN
 305          ISIDE = 3
 306          COMPL = .TRUE.
 307          COMPR = .TRUE.
 308       ELSE
 309          ISIDE = -1
 310       END IF
 311 *
 312       INFO = 0
 313       IF( ISIDE.LT.0 ) THEN
 314          INFO = -1
 315       ELSE IF( IHWMNY.LT.0 ) THEN
 316          INFO = -2
 317       ELSE IF( N.LT.0 ) THEN
 318          INFO = -4
 319       ELSE IF( LDS.LT.MAX( 1, N ) ) THEN
 320          INFO = -6
 321       ELSE IF( LDP.LT.MAX( 1, N ) ) THEN
 322          INFO = -8
 323       END IF
 324       IF( INFO.NE.0 ) THEN
 325          CALL XERBLA( 'CTGEVC', -INFO )
 326          RETURN
 327       END IF
 328 *
 329 *     Count the number of eigenvectors
 330 *
 331       IF( .NOT.ILALL ) THEN
 332          IM = 0
 333          DO 10 J = 1, N
 334             IF( SELECT( J ) )
 335      $         IM = IM + 1
 336    10    CONTINUE
 337       ELSE
 338          IM = N
 339       END IF
 340 *
 341 *     Check diagonal of B
 342 *
 343       ILBBAD = .FALSE.
 344       DO 20 J = 1, N
 345          IF( AIMAG( P( J, J ) ).NE.ZERO )
 346      $      ILBBAD = .TRUE.
 347    20 CONTINUE
 348 *
 349       IF( ILBBAD ) THEN
 350          INFO = -7
 351       ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN
 352          INFO = -10
 353       ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN
 354          INFO = -12
 355       ELSE IF( MM.LT.IM ) THEN
 356          INFO = -13
 357       END IF
 358       IF( INFO.NE.0 ) THEN
 359          CALL XERBLA( 'CTGEVC', -INFO )
 360          RETURN
 361       END IF
 362 *
 363 *     Quick return if possible
 364 *
 365       M = IM
 366       IF( N.EQ.0 )
 367      $   RETURN
 368 *
 369 *     Machine Constants
 370 *
 371       SAFMIN = SLAMCH( 'Safe minimum' )
 372       BIG = ONE / SAFMIN
 373       CALL SLABAD( SAFMIN, BIG )
 374       ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
 375       SMALL = SAFMIN*N / ULP
 376       BIG = ONE / SMALL
 377       BIGNUM = ONE / ( SAFMIN*N )
 378 *
 379 *     Compute the 1-norm of each column of the strictly upper triangular
 380 *     part of A and B to check for possible overflow in the triangular
 381 *     solver.
 382 *
 383       ANORM = ABS1( S( 1, 1 ) )
 384       BNORM = ABS1( P( 1, 1 ) )
 385       RWORK( 1 ) = ZERO
 386       RWORK( N+1 ) = ZERO
 387       DO 40 J = 2, N
 388          RWORK( J ) = ZERO
 389          RWORK( N+J ) = ZERO
 390          DO 30 I = 1, J - 1
 391             RWORK( J ) = RWORK( J ) + ABS1( S( I, J ) )
 392             RWORK( N+J ) = RWORK( N+J ) + ABS1( P( I, J ) )
 393    30    CONTINUE
 394          ANORM = MAX( ANORM, RWORK( J )+ABS1( S( J, J ) ) )
 395          BNORM = MAX( BNORM, RWORK( N+J )+ABS1( P( J, J ) ) )
 396    40 CONTINUE
 397 *
 398       ASCALE = ONE / MAX( ANORM, SAFMIN )
 399       BSCALE = ONE / MAX( BNORM, SAFMIN )
 400 *
 401 *     Left eigenvectors
 402 *
 403       IF( COMPL ) THEN
 404          IEIG = 0
 405 *
 406 *        Main loop over eigenvalues
 407 *
 408          DO 140 JE = 1, N
 409             IF( ILALL ) THEN
 410                ILCOMP = .TRUE.
 411             ELSE
 412                ILCOMP = SELECT( JE )
 413             END IF
 414             IF( ILCOMP ) THEN
 415                IEIG = IEIG + 1
 416 *
 417                IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
 418      $             ABS( REAL( P( JE, JE ) ) ).LE.SAFMIN ) THEN
 419 *
 420 *                 Singular matrix pencil -- return unit eigenvector
 421 *
 422                   DO 50 JR = 1, N
 423                      VL( JR, IEIG ) = CZERO
 424    50             CONTINUE
 425                   VL( IEIG, IEIG ) = CONE
 426                   GO TO 140
 427                END IF
 428 *
 429 *              Non-singular eigenvalue:
 430 *              Compute coefficients  a  and  b  in
 431 *                   H
 432 *                 y  ( a A - b B ) = 0
 433 *
 434                TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE,
 435      $                ABS( REAL( P( JE, JE ) ) )*BSCALE, SAFMIN )
 436                SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
 437                SBETA = ( TEMP*REAL( P( JE, JE ) ) )*BSCALE
 438                ACOEFF = SBETA*ASCALE
 439                BCOEFF = SALPHA*BSCALE
 440 *
 441 *              Scale to avoid underflow
 442 *
 443                LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL
 444                LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT.
 445      $               SMALL
 446 *
 447                SCALE = ONE
 448                IF( LSA )
 449      $            SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
 450                IF( LSB )
 451      $            SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )*
 452      $                    MIN( BNORM, BIG ) )
 453                IF( LSA .OR. LSB ) THEN
 454                   SCALE = MIN( SCALE, ONE /
 455      $                    ( SAFMIN*MAX( ONE, ABS( ACOEFF ),
 456      $                    ABS1( BCOEFF ) ) ) )
 457                   IF( LSA ) THEN
 458                      ACOEFF = ASCALE*( SCALE*SBETA )
 459                   ELSE
 460                      ACOEFF = SCALE*ACOEFF
 461                   END IF
 462                   IF( LSB ) THEN
 463                      BCOEFF = BSCALE*( SCALE*SALPHA )
 464                   ELSE
 465                      BCOEFF = SCALE*BCOEFF
 466                   END IF
 467                END IF
 468 *
 469                ACOEFA = ABS( ACOEFF )
 470                BCOEFA = ABS1( BCOEFF )
 471                XMAX = ONE
 472                DO 60 JR = 1, N
 473                   WORK( JR ) = CZERO
 474    60          CONTINUE
 475                WORK( JE ) = CONE
 476                DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
 477 *
 478 *                                              H
 479 *              Triangular solve of  (a A - b B)  y = 0
 480 *
 481 *                                      H
 482 *              (rowwise in  (a A - b B) , or columnwise in a A - b B)
 483 *
 484                DO 100 J = JE + 1, N
 485 *
 486 *                 Compute
 487 *                       j-1
 488 *                 SUM = sum  conjg( a*S(k,j) - b*P(k,j) )*x(k)
 489 *                       k=je
 490 *                 (Scale if necessary)
 491 *
 492                   TEMP = ONE / XMAX
 493                   IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GT.BIGNUM*
 494      $                TEMP ) THEN
 495                      DO 70 JR = JE, J - 1
 496                         WORK( JR ) = TEMP*WORK( JR )
 497    70                CONTINUE
 498                      XMAX = ONE
 499                   END IF
 500                   SUMA = CZERO
 501                   SUMB = CZERO
 502 *
 503                   DO 80 JR = JE, J - 1
 504                      SUMA = SUMA + CONJG( S( JR, J ) )*WORK( JR )
 505                      SUMB = SUMB + CONJG( P( JR, J ) )*WORK( JR )
 506    80             CONTINUE
 507                   SUM = ACOEFF*SUMA - CONJG( BCOEFF )*SUMB
 508 *
 509 *                 Form x(j) = - SUM / conjg( a*S(j,j) - b*P(j,j) )
 510 *
 511 *                 with scaling and perturbation of the denominator
 512 *
 513                   D = CONJG( ACOEFF*S( J, J )-BCOEFF*P( J, J ) )
 514                   IF( ABS1( D ).LE.DMIN )
 515      $               D = CMPLX( DMIN )
 516 *
 517                   IF( ABS1( D ).LT.ONE ) THEN
 518                      IF( ABS1( SUM ).GE.BIGNUM*ABS1( D ) ) THEN
 519                         TEMP = ONE / ABS1( SUM )
 520                         DO 90 JR = JE, J - 1
 521                            WORK( JR ) = TEMP*WORK( JR )
 522    90                   CONTINUE
 523                         XMAX = TEMP*XMAX
 524                         SUM = TEMP*SUM
 525                      END IF
 526                   END IF
 527                   WORK( J ) = CLADIV( -SUM, D )
 528                   XMAX = MAX( XMAX, ABS1( WORK( J ) ) )
 529   100          CONTINUE
 530 *
 531 *              Back transform eigenvector if HOWMNY='B'.
 532 *
 533                IF( ILBACK ) THEN
 534                   CALL CGEMV( 'N', N, N+1-JE, CONE, VL( 1, JE ), LDVL,
 535      $                        WORK( JE ), 1, CZERO, WORK( N+1 ), 1 )
 536                   ISRC = 2
 537                   IBEG = 1
 538                ELSE
 539                   ISRC = 1
 540                   IBEG = JE
 541                END IF
 542 *
 543 *              Copy and scale eigenvector into column of VL
 544 *
 545                XMAX = ZERO
 546                DO 110 JR = IBEG, N
 547                   XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) )
 548   110          CONTINUE
 549 *
 550                IF( XMAX.GT.SAFMIN ) THEN
 551                   TEMP = ONE / XMAX
 552                   DO 120 JR = IBEG, N
 553                      VL( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR )
 554   120             CONTINUE
 555                ELSE
 556                   IBEG = N + 1
 557                END IF
 558 *
 559                DO 130 JR = 1, IBEG - 1
 560                   VL( JR, IEIG ) = CZERO
 561   130          CONTINUE
 562 *
 563             END IF
 564   140    CONTINUE
 565       END IF
 566 *
 567 *     Right eigenvectors
 568 *
 569       IF( COMPR ) THEN
 570          IEIG = IM + 1
 571 *
 572 *        Main loop over eigenvalues
 573 *
 574          DO 250 JE = N, 1, -1
 575             IF( ILALL ) THEN
 576                ILCOMP = .TRUE.
 577             ELSE
 578                ILCOMP = SELECT( JE )
 579             END IF
 580             IF( ILCOMP ) THEN
 581                IEIG = IEIG - 1
 582 *
 583                IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
 584      $             ABS( REAL( P( JE, JE ) ) ).LE.SAFMIN ) THEN
 585 *
 586 *                 Singular matrix pencil -- return unit eigenvector
 587 *
 588                   DO 150 JR = 1, N
 589                      VR( JR, IEIG ) = CZERO
 590   150             CONTINUE
 591                   VR( IEIG, IEIG ) = CONE
 592                   GO TO 250
 593                END IF
 594 *
 595 *              Non-singular eigenvalue:
 596 *              Compute coefficients  a  and  b  in
 597 *
 598 *              ( a A - b B ) x  = 0
 599 *
 600                TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE,
 601      $                ABS( REAL( P( JE, JE ) ) )*BSCALE, SAFMIN )
 602                SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
 603                SBETA = ( TEMP*REAL( P( JE, JE ) ) )*BSCALE
 604                ACOEFF = SBETA*ASCALE
 605                BCOEFF = SALPHA*BSCALE
 606 *
 607 *              Scale to avoid underflow
 608 *
 609                LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL
 610                LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT.
 611      $               SMALL
 612 *
 613                SCALE = ONE
 614                IF( LSA )
 615      $            SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
 616                IF( LSB )
 617      $            SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )*
 618      $                    MIN( BNORM, BIG ) )
 619                IF( LSA .OR. LSB ) THEN
 620                   SCALE = MIN( SCALE, ONE /
 621      $                    ( SAFMIN*MAX( ONE, ABS( ACOEFF ),
 622      $                    ABS1( BCOEFF ) ) ) )
 623                   IF( LSA ) THEN
 624                      ACOEFF = ASCALE*( SCALE*SBETA )
 625                   ELSE
 626                      ACOEFF = SCALE*ACOEFF
 627                   END IF
 628                   IF( LSB ) THEN
 629                      BCOEFF = BSCALE*( SCALE*SALPHA )
 630                   ELSE
 631                      BCOEFF = SCALE*BCOEFF
 632                   END IF
 633                END IF
 634 *
 635                ACOEFA = ABS( ACOEFF )
 636                BCOEFA = ABS1( BCOEFF )
 637                XMAX = ONE
 638                DO 160 JR = 1, N
 639                   WORK( JR ) = CZERO
 640   160          CONTINUE
 641                WORK( JE ) = CONE
 642                DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
 643 *
 644 *              Triangular solve of  (a A - b B) x = 0  (columnwise)
 645 *
 646 *              WORK(1:j-1) contains sums w,
 647 *              WORK(j+1:JE) contains x
 648 *
 649                DO 170 JR = 1, JE - 1
 650                   WORK( JR ) = ACOEFF*S( JR, JE ) - BCOEFF*P( JR, JE )
 651   170          CONTINUE
 652                WORK( JE ) = CONE
 653 *
 654                DO 210 J = JE - 1, 1, -1
 655 *
 656 *                 Form x(j) := - w(j) / d
 657 *                 with scaling and perturbation of the denominator
 658 *
 659                   D = ACOEFF*S( J, J ) - BCOEFF*P( J, J )
 660                   IF( ABS1( D ).LE.DMIN )
 661      $               D = CMPLX( DMIN )
 662 *
 663                   IF( ABS1( D ).LT.ONE ) THEN
 664                      IF( ABS1( WORK( J ) ).GE.BIGNUM*ABS1( D ) ) THEN
 665                         TEMP = ONE / ABS1( WORK( J ) )
 666                         DO 180 JR = 1, JE
 667                            WORK( JR ) = TEMP*WORK( JR )
 668   180                   CONTINUE
 669                      END IF
 670                   END IF
 671 *
 672                   WORK( J ) = CLADIV( -WORK( J ), D )
 673 *
 674                   IF( J.GT.1 ) THEN
 675 *
 676 *                    w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
 677 *
 678                      IF( ABS1( WORK( J ) ).GT.ONE ) THEN
 679                         TEMP = ONE / ABS1( WORK( J ) )
 680                         IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GE.
 681      $                      BIGNUM*TEMP ) THEN
 682                            DO 190 JR = 1, JE
 683                               WORK( JR ) = TEMP*WORK( JR )
 684   190                      CONTINUE
 685                         END IF
 686                      END IF
 687 *
 688                      CA = ACOEFF*WORK( J )
 689                      CB = BCOEFF*WORK( J )
 690                      DO 200 JR = 1, J - 1
 691                         WORK( JR ) = WORK( JR ) + CA*S( JR, J ) -
 692      $                               CB*P( JR, J )
 693   200                CONTINUE
 694                   END IF
 695   210          CONTINUE
 696 *
 697 *              Back transform eigenvector if HOWMNY='B'.
 698 *
 699                IF( ILBACK ) THEN
 700                   CALL CGEMV( 'N', N, JE, CONE, VR, LDVR, WORK, 1,
 701      $                        CZERO, WORK( N+1 ), 1 )
 702                   ISRC = 2
 703                   IEND = N
 704                ELSE
 705                   ISRC = 1
 706                   IEND = JE
 707                END IF
 708 *
 709 *              Copy and scale eigenvector into column of VR
 710 *
 711                XMAX = ZERO
 712                DO 220 JR = 1, IEND
 713                   XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) )
 714   220          CONTINUE
 715 *
 716                IF( XMAX.GT.SAFMIN ) THEN
 717                   TEMP = ONE / XMAX
 718                   DO 230 JR = 1, IEND
 719                      VR( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR )
 720   230             CONTINUE
 721                ELSE
 722                   IEND = 0
 723                END IF
 724 *
 725                DO 240 JR = IEND + 1, N
 726                   VR( JR, IEIG ) = CZERO
 727   240          CONTINUE
 728 *
 729             END IF
 730   250    CONTINUE
 731       END IF
 732 *
 733       RETURN
 734 *
 735 *     End of CTGEVC
 736 *
 737       END