LAPACK: SRC/ctgevc.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
00002      $                   LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          HOWMNY, SIDE
00011       INTEGER            INFO, LDP, LDS, LDVL, LDVR, M, MM, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       LOGICAL            SELECT( * )
00015       REAL               RWORK( * )
00016       COMPLEX            P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
00017      $                   VR( LDVR, * ), WORK( * )
00018 *     ..
00019 *
00020 *
00021 *  Purpose
00022 *  =======
00023 *
00024 *  CTGEVC computes some or all of the right and/or left eigenvectors of
00025 *  a pair of complex matrices (S,P), where S and P are upper triangular.
00026 *  Matrix pairs of this type are produced by the generalized Schur
00027 *  factorization of a complex matrix pair (A,B):
00028 *  
00029 *     A = Q*S*Z**H,  B = Q*P*Z**H
00030 *  
00031 *  as computed by CGGHRD + CHGEQZ.
00032 *  
00033 *  The right eigenvector x and the left eigenvector y of (S,P)
00034 *  corresponding to an eigenvalue w are defined by:
00035 *  
00036 *     S*x = w*P*x,  (y**H)*S = w*(y**H)*P,
00037 *  
00038 *  where y**H denotes the conjugate tranpose of y.
00039 *  The eigenvalues are not input to this routine, but are computed
00040 *  directly from the diagonal elements of S and P.
00041 *  
00042 *  This routine returns the matrices X and/or Y of right and left
00043 *  eigenvectors of (S,P), or the products Z*X and/or Q*Y,
00044 *  where Z and Q are input matrices.
00045 *  If Q and Z are the unitary factors from the generalized Schur
00046 *  factorization of a matrix pair (A,B), then Z*X and Q*Y
00047 *  are the matrices of right and left eigenvectors of (A,B).
00048 *
00049 *  Arguments
00050 *  =========
00051 *
00052 *  SIDE    (input) CHARACTER*1
00053 *          = 'R': compute right eigenvectors only;
00054 *          = 'L': compute left eigenvectors only;
00055 *          = 'B': compute both right and left eigenvectors.
00056 *
00057 *  HOWMNY  (input) CHARACTER*1
00058 *          = 'A': compute all right and/or left eigenvectors;
00059 *          = 'B': compute all right and/or left eigenvectors,
00060 *                 backtransformed by the matrices in VR and/or VL;
00061 *          = 'S': compute selected right and/or left eigenvectors,
00062 *                 specified by the logical array SELECT.
00063 *
00064 *  SELECT  (input) LOGICAL array, dimension (N)
00065 *          If HOWMNY='S', SELECT specifies the eigenvectors to be
00066 *          computed.  The eigenvector corresponding to the j-th
00067 *          eigenvalue is computed if SELECT(j) = .TRUE..
00068 *          Not referenced if HOWMNY = 'A' or 'B'.
00069 *
00070 *  N       (input) INTEGER
00071 *          The order of the matrices S and P.  N >= 0.
00072 *
00073 *  S       (input) COMPLEX array, dimension (LDS,N)
00074 *          The upper triangular matrix S from a generalized Schur
00075 *          factorization, as computed by CHGEQZ.
00076 *
00077 *  LDS     (input) INTEGER
00078 *          The leading dimension of array S.  LDS >= max(1,N).
00079 *
00080 *  P       (input) COMPLEX array, dimension (LDP,N)
00081 *          The upper triangular matrix P from a generalized Schur
00082 *          factorization, as computed by CHGEQZ.  P must have real
00083 *          diagonal elements.
00084 *
00085 *  LDP     (input) INTEGER
00086 *          The leading dimension of array P.  LDP >= max(1,N).
00087 *
00088 *  VL      (input/output) COMPLEX array, dimension (LDVL,MM)
00089 *          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
00090 *          contain an N-by-N matrix Q (usually the unitary matrix Q
00091 *          of left Schur vectors returned by CHGEQZ).
00092 *          On exit, if SIDE = 'L' or 'B', VL contains:
00093 *          if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
00094 *          if HOWMNY = 'B', the matrix Q*Y;
00095 *          if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
00096 *                      SELECT, stored consecutively in the columns of
00097 *                      VL, in the same order as their eigenvalues.
00098 *          Not referenced if SIDE = 'R'.
00099 *
00100 *  LDVL    (input) INTEGER
00101 *          The leading dimension of array VL.  LDVL >= 1, and if
00102 *          SIDE = 'L' or 'l' or 'B' or 'b', LDVL >= N.
00103 *
00104 *  VR      (input/output) COMPLEX array, dimension (LDVR,MM)
00105 *          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
00106 *          contain an N-by-N matrix Q (usually the unitary matrix Z
00107 *          of right Schur vectors returned by CHGEQZ).
00108 *          On exit, if SIDE = 'R' or 'B', VR contains:
00109 *          if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
00110 *          if HOWMNY = 'B', the matrix Z*X;
00111 *          if HOWMNY = 'S', the right eigenvectors of (S,P) specified by
00112 *                      SELECT, stored consecutively in the columns of
00113 *                      VR, in the same order as their eigenvalues.
00114 *          Not referenced if SIDE = 'L'.
00115 *
00116 *  LDVR    (input) INTEGER
00117 *          The leading dimension of the array VR.  LDVR >= 1, and if
00118 *          SIDE = 'R' or 'B', LDVR >= N.
00119 *
00120 *  MM      (input) INTEGER
00121 *          The number of columns in the arrays VL and/or VR. MM >= M.
00122 *
00123 *  M       (output) INTEGER
00124 *          The number of columns in the arrays VL and/or VR actually
00125 *          used to store the eigenvectors.  If HOWMNY = 'A' or 'B', M
00126 *          is set to N.  Each selected eigenvector occupies one column.
00127 *
00128 *  WORK    (workspace) COMPLEX array, dimension (2*N)
00129 *
00130 *  RWORK   (workspace) REAL array, dimension (2*N)
00131 *
00132 *  INFO    (output) INTEGER
00133 *          = 0:  successful exit.
00134 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00135 *
00136 *  =====================================================================
00137 *
00138 *     .. Parameters ..
00139       REAL               ZERO, ONE
00140       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00141       COMPLEX            CZERO, CONE
00142       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
00143      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
00144 *     ..
00145 *     .. Local Scalars ..
00146       LOGICAL            COMPL, COMPR, ILALL, ILBACK, ILBBAD, ILCOMP,
00147      $                   LSA, LSB
00148       INTEGER            I, IBEG, IEIG, IEND, IHWMNY, IM, ISIDE, ISRC,
00149      $                   J, JE, JR
00150       REAL               ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG,
00151      $                   BIGNUM, BNORM, BSCALE, DMIN, SAFMIN, SBETA,
00152      $                   SCALE, SMALL, TEMP, ULP, XMAX
00153       COMPLEX            BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X
00154 *     ..
00155 *     .. External Functions ..
00156       LOGICAL            LSAME
00157       REAL               SLAMCH
00158       COMPLEX            CLADIV
00159       EXTERNAL           LSAME, SLAMCH, CLADIV
00160 *     ..
00161 *     .. External Subroutines ..
00162       EXTERNAL           CGEMV, SLABAD, XERBLA
00163 *     ..
00164 *     .. Intrinsic Functions ..
00165       INTRINSIC          ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL
00166 *     ..
00167 *     .. Statement Functions ..
00168       REAL               ABS1
00169 *     ..
00170 *     .. Statement Function definitions ..
00171       ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
00172 *     ..
00173 *     .. Executable Statements ..
00174 *
00175 *     Decode and Test the input parameters
00176 *
00177       IF( LSAME( HOWMNY, 'A' ) ) THEN
00178          IHWMNY = 1
00179          ILALL = .TRUE.
00180          ILBACK = .FALSE.
00181       ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN
00182          IHWMNY = 2
00183          ILALL = .FALSE.
00184          ILBACK = .FALSE.
00185       ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN
00186          IHWMNY = 3
00187          ILALL = .TRUE.
00188          ILBACK = .TRUE.
00189       ELSE
00190          IHWMNY = -1
00191       END IF
00192 *
00193       IF( LSAME( SIDE, 'R' ) ) THEN
00194          ISIDE = 1
00195          COMPL = .FALSE.
00196          COMPR = .TRUE.
00197       ELSE IF( LSAME( SIDE, 'L' ) ) THEN
00198          ISIDE = 2
00199          COMPL = .TRUE.
00200          COMPR = .FALSE.
00201       ELSE IF( LSAME( SIDE, 'B' ) ) THEN
00202          ISIDE = 3
00203          COMPL = .TRUE.
00204          COMPR = .TRUE.
00205       ELSE
00206          ISIDE = -1
00207       END IF
00208 *
00209       INFO = 0
00210       IF( ISIDE.LT.0 ) THEN
00211          INFO = -1
00212       ELSE IF( IHWMNY.LT.0 ) THEN
00213          INFO = -2
00214       ELSE IF( N.LT.0 ) THEN
00215          INFO = -4
00216       ELSE IF( LDS.LT.MAX( 1, N ) ) THEN
00217          INFO = -6
00218       ELSE IF( LDP.LT.MAX( 1, N ) ) THEN
00219          INFO = -8
00220       END IF
00221       IF( INFO.NE.0 ) THEN
00222          CALL XERBLA( 'CTGEVC', -INFO )
00223          RETURN
00224       END IF
00225 *
00226 *     Count the number of eigenvectors
00227 *
00228       IF( .NOT.ILALL ) THEN
00229          IM = 0
00230          DO 10 J = 1, N
00231             IF( SELECT( J ) )
00232      $         IM = IM + 1
00233    10    CONTINUE
00234       ELSE
00235          IM = N
00236       END IF
00237 *
00238 *     Check diagonal of B
00239 *
00240       ILBBAD = .FALSE.
00241       DO 20 J = 1, N
00242          IF( AIMAG( P( J, J ) ).NE.ZERO )
00243      $      ILBBAD = .TRUE.
00244    20 CONTINUE
00245 *
00246       IF( ILBBAD ) THEN
00247          INFO = -7
00248       ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN
00249          INFO = -10
00250       ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN
00251          INFO = -12
00252       ELSE IF( MM.LT.IM ) THEN
00253          INFO = -13
00254       END IF
00255       IF( INFO.NE.0 ) THEN
00256          CALL XERBLA( 'CTGEVC', -INFO )
00257          RETURN
00258       END IF
00259 *
00260 *     Quick return if possible
00261 *
00262       M = IM
00263       IF( N.EQ.0 )
00264      $   RETURN
00265 *
00266 *     Machine Constants
00267 *
00268       SAFMIN = SLAMCH( 'Safe minimum' )
00269       BIG = ONE / SAFMIN
00270       CALL SLABAD( SAFMIN, BIG )
00271       ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
00272       SMALL = SAFMIN*N / ULP
00273       BIG = ONE / SMALL
00274       BIGNUM = ONE / ( SAFMIN*N )
00275 *
00276 *     Compute the 1-norm of each column of the strictly upper triangular
00277 *     part of A and B to check for possible overflow in the triangular
00278 *     solver.
00279 *
00280       ANORM = ABS1( S( 1, 1 ) )
00281       BNORM = ABS1( P( 1, 1 ) )
00282       RWORK( 1 ) = ZERO
00283       RWORK( N+1 ) = ZERO
00284       DO 40 J = 2, N
00285          RWORK( J ) = ZERO
00286          RWORK( N+J ) = ZERO
00287          DO 30 I = 1, J - 1
00288             RWORK( J ) = RWORK( J ) + ABS1( S( I, J ) )
00289             RWORK( N+J ) = RWORK( N+J ) + ABS1( P( I, J ) )
00290    30    CONTINUE
00291          ANORM = MAX( ANORM, RWORK( J )+ABS1( S( J, J ) ) )
00292          BNORM = MAX( BNORM, RWORK( N+J )+ABS1( P( J, J ) ) )
00293    40 CONTINUE
00294 *
00295       ASCALE = ONE / MAX( ANORM, SAFMIN )
00296       BSCALE = ONE / MAX( BNORM, SAFMIN )
00297 *
00298 *     Left eigenvectors
00299 *
00300       IF( COMPL ) THEN
00301          IEIG = 0
00302 *
00303 *        Main loop over eigenvalues
00304 *
00305          DO 140 JE = 1, N
00306             IF( ILALL ) THEN
00307                ILCOMP = .TRUE.
00308             ELSE
00309                ILCOMP = SELECT( JE )
00310             END IF
00311             IF( ILCOMP ) THEN
00312                IEIG = IEIG + 1
00313 *
00314                IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
00315      $             ABS( REAL( P( JE, JE ) ) ).LE.SAFMIN ) THEN
00316 *
00317 *                 Singular matrix pencil -- return unit eigenvector
00318 *
00319                   DO 50 JR = 1, N
00320                      VL( JR, IEIG ) = CZERO
00321    50             CONTINUE
00322                   VL( IEIG, IEIG ) = CONE
00323                   GO TO 140
00324                END IF
00325 *
00326 *              Non-singular eigenvalue:
00327 *              Compute coefficients  a  and  b  in
00328 *                   H
00329 *                 y  ( a A - b B ) = 0
00330 *
00331                TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE,
00332      $                ABS( REAL( P( JE, JE ) ) )*BSCALE, SAFMIN )
00333                SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
00334                SBETA = ( TEMP*REAL( P( JE, JE ) ) )*BSCALE
00335                ACOEFF = SBETA*ASCALE
00336                BCOEFF = SALPHA*BSCALE
00337 *
00338 *              Scale to avoid underflow
00339 *
00340                LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL
00341                LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT.
00342      $               SMALL
00343 *
00344                SCALE = ONE
00345                IF( LSA )
00346      $            SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
00347                IF( LSB )
00348      $            SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )*
00349      $                    MIN( BNORM, BIG ) )
00350                IF( LSA .OR. LSB ) THEN
00351                   SCALE = MIN( SCALE, ONE /
00352      $                    ( SAFMIN*MAX( ONE, ABS( ACOEFF ),
00353      $                    ABS1( BCOEFF ) ) ) )
00354                   IF( LSA ) THEN
00355                      ACOEFF = ASCALE*( SCALE*SBETA )
00356                   ELSE
00357                      ACOEFF = SCALE*ACOEFF
00358                   END IF
00359                   IF( LSB ) THEN
00360                      BCOEFF = BSCALE*( SCALE*SALPHA )
00361                   ELSE
00362                      BCOEFF = SCALE*BCOEFF
00363                   END IF
00364                END IF
00365 *
00366                ACOEFA = ABS( ACOEFF )
00367                BCOEFA = ABS1( BCOEFF )
00368                XMAX = ONE
00369                DO 60 JR = 1, N
00370                   WORK( JR ) = CZERO
00371    60          CONTINUE
00372                WORK( JE ) = CONE
00373                DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
00374 *
00375 *                                              H
00376 *              Triangular solve of  (a A - b B)  y = 0
00377 *
00378 *                                      H
00379 *              (rowwise in  (a A - b B) , or columnwise in a A - b B)
00380 *
00381                DO 100 J = JE + 1, N
00382 *
00383 *                 Compute
00384 *                       j-1
00385 *                 SUM = sum  conjg( a*S(k,j) - b*P(k,j) )*x(k)
00386 *                       k=je
00387 *                 (Scale if necessary)
00388 *
00389                   TEMP = ONE / XMAX
00390                   IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GT.BIGNUM*
00391      $                TEMP ) THEN
00392                      DO 70 JR = JE, J - 1
00393                         WORK( JR ) = TEMP*WORK( JR )
00394    70                CONTINUE
00395                      XMAX = ONE
00396                   END IF
00397                   SUMA = CZERO
00398                   SUMB = CZERO
00399 *
00400                   DO 80 JR = JE, J - 1
00401                      SUMA = SUMA + CONJG( S( JR, J ) )*WORK( JR )
00402                      SUMB = SUMB + CONJG( P( JR, J ) )*WORK( JR )
00403    80             CONTINUE
00404                   SUM = ACOEFF*SUMA - CONJG( BCOEFF )*SUMB
00405 *
00406 *                 Form x(j) = - SUM / conjg( a*S(j,j) - b*P(j,j) )
00407 *
00408 *                 with scaling and perturbation of the denominator
00409 *
00410                   D = CONJG( ACOEFF*S( J, J )-BCOEFF*P( J, J ) )
00411                   IF( ABS1( D ).LE.DMIN )
00412      $               D = CMPLX( DMIN )
00413 *
00414                   IF( ABS1( D ).LT.ONE ) THEN
00415                      IF( ABS1( SUM ).GE.BIGNUM*ABS1( D ) ) THEN
00416                         TEMP = ONE / ABS1( SUM )
00417                         DO 90 JR = JE, J - 1
00418                            WORK( JR ) = TEMP*WORK( JR )
00419    90                   CONTINUE
00420                         XMAX = TEMP*XMAX
00421                         SUM = TEMP*SUM
00422                      END IF
00423                   END IF
00424                   WORK( J ) = CLADIV( -SUM, D )
00425                   XMAX = MAX( XMAX, ABS1( WORK( J ) ) )
00426   100          CONTINUE
00427 *
00428 *              Back transform eigenvector if HOWMNY='B'.
00429 *
00430                IF( ILBACK ) THEN
00431                   CALL CGEMV( 'N', N, N+1-JE, CONE, VL( 1, JE ), LDVL,
00432      $                        WORK( JE ), 1, CZERO, WORK( N+1 ), 1 )
00433                   ISRC = 2
00434                   IBEG = 1
00435                ELSE
00436                   ISRC = 1
00437                   IBEG = JE
00438                END IF
00439 *
00440 *              Copy and scale eigenvector into column of VL
00441 *
00442                XMAX = ZERO
00443                DO 110 JR = IBEG, N
00444                   XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) )
00445   110          CONTINUE
00446 *
00447                IF( XMAX.GT.SAFMIN ) THEN
00448                   TEMP = ONE / XMAX
00449                   DO 120 JR = IBEG, N
00450                      VL( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR )
00451   120             CONTINUE
00452                ELSE
00453                   IBEG = N + 1
00454                END IF
00455 *
00456                DO 130 JR = 1, IBEG - 1
00457                   VL( JR, IEIG ) = CZERO
00458   130          CONTINUE
00459 *
00460             END IF
00461   140    CONTINUE
00462       END IF
00463 *
00464 *     Right eigenvectors
00465 *
00466       IF( COMPR ) THEN
00467          IEIG = IM + 1
00468 *
00469 *        Main loop over eigenvalues
00470 *
00471          DO 250 JE = N, 1, -1
00472             IF( ILALL ) THEN
00473                ILCOMP = .TRUE.
00474             ELSE
00475                ILCOMP = SELECT( JE )
00476             END IF
00477             IF( ILCOMP ) THEN
00478                IEIG = IEIG - 1
00479 *
00480                IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
00481      $             ABS( REAL( P( JE, JE ) ) ).LE.SAFMIN ) THEN
00482 *
00483 *                 Singular matrix pencil -- return unit eigenvector
00484 *
00485                   DO 150 JR = 1, N
00486                      VR( JR, IEIG ) = CZERO
00487   150             CONTINUE
00488                   VR( IEIG, IEIG ) = CONE
00489                   GO TO 250
00490                END IF
00491 *
00492 *              Non-singular eigenvalue:
00493 *              Compute coefficients  a  and  b  in
00494 *
00495 *              ( a A - b B ) x  = 0
00496 *
00497                TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE,
00498      $                ABS( REAL( P( JE, JE ) ) )*BSCALE, SAFMIN )
00499                SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
00500                SBETA = ( TEMP*REAL( P( JE, JE ) ) )*BSCALE
00501                ACOEFF = SBETA*ASCALE
00502                BCOEFF = SALPHA*BSCALE
00503 *
00504 *              Scale to avoid underflow
00505 *
00506                LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL
00507                LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT.
00508      $               SMALL
00509 *
00510                SCALE = ONE
00511                IF( LSA )
00512      $            SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
00513                IF( LSB )
00514      $            SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )*
00515      $                    MIN( BNORM, BIG ) )
00516                IF( LSA .OR. LSB ) THEN
00517                   SCALE = MIN( SCALE, ONE /
00518      $                    ( SAFMIN*MAX( ONE, ABS( ACOEFF ),
00519      $                    ABS1( BCOEFF ) ) ) )
00520                   IF( LSA ) THEN
00521                      ACOEFF = ASCALE*( SCALE*SBETA )
00522                   ELSE
00523                      ACOEFF = SCALE*ACOEFF
00524                   END IF
00525                   IF( LSB ) THEN
00526                      BCOEFF = BSCALE*( SCALE*SALPHA )
00527                   ELSE
00528                      BCOEFF = SCALE*BCOEFF
00529                   END IF
00530                END IF
00531 *
00532                ACOEFA = ABS( ACOEFF )
00533                BCOEFA = ABS1( BCOEFF )
00534                XMAX = ONE
00535                DO 160 JR = 1, N
00536                   WORK( JR ) = CZERO
00537   160          CONTINUE
00538                WORK( JE ) = CONE
00539                DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
00540 *
00541 *              Triangular solve of  (a A - b B) x = 0  (columnwise)
00542 *
00543 *              WORK(1:j-1) contains sums w,
00544 *              WORK(j+1:JE) contains x
00545 *
00546                DO 170 JR = 1, JE - 1
00547                   WORK( JR ) = ACOEFF*S( JR, JE ) - BCOEFF*P( JR, JE )
00548   170          CONTINUE
00549                WORK( JE ) = CONE
00550 *
00551                DO 210 J = JE - 1, 1, -1
00552 *
00553 *                 Form x(j) := - w(j) / d
00554 *                 with scaling and perturbation of the denominator
00555 *
00556                   D = ACOEFF*S( J, J ) - BCOEFF*P( J, J )
00557                   IF( ABS1( D ).LE.DMIN )
00558      $               D = CMPLX( DMIN )
00559 *
00560                   IF( ABS1( D ).LT.ONE ) THEN
00561                      IF( ABS1( WORK( J ) ).GE.BIGNUM*ABS1( D ) ) THEN
00562                         TEMP = ONE / ABS1( WORK( J ) )
00563                         DO 180 JR = 1, JE
00564                            WORK( JR ) = TEMP*WORK( JR )
00565   180                   CONTINUE
00566                      END IF
00567                   END IF
00568 *
00569                   WORK( J ) = CLADIV( -WORK( J ), D )
00570 *
00571                   IF( J.GT.1 ) THEN
00572 *
00573 *                    w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
00574 *
00575                      IF( ABS1( WORK( J ) ).GT.ONE ) THEN
00576                         TEMP = ONE / ABS1( WORK( J ) )
00577                         IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GE.
00578      $                      BIGNUM*TEMP ) THEN
00579                            DO 190 JR = 1, JE
00580                               WORK( JR ) = TEMP*WORK( JR )
00581   190                      CONTINUE
00582                         END IF
00583                      END IF
00584 *
00585                      CA = ACOEFF*WORK( J )
00586                      CB = BCOEFF*WORK( J )
00587                      DO 200 JR = 1, J - 1
00588                         WORK( JR ) = WORK( JR ) + CA*S( JR, J ) -
00589      $                               CB*P( JR, J )
00590   200                CONTINUE
00591                   END IF
00592   210          CONTINUE
00593 *
00594 *              Back transform eigenvector if HOWMNY='B'.
00595 *
00596                IF( ILBACK ) THEN
00597                   CALL CGEMV( 'N', N, JE, CONE, VR, LDVR, WORK, 1,
00598      $                        CZERO, WORK( N+1 ), 1 )
00599                   ISRC = 2
00600                   IEND = N
00601                ELSE
00602                   ISRC = 1
00603                   IEND = JE
00604                END IF
00605 *
00606 *              Copy and scale eigenvector into column of VR
00607 *
00608                XMAX = ZERO
00609                DO 220 JR = 1, IEND
00610                   XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) )
00611   220          CONTINUE
00612 *
00613                IF( XMAX.GT.SAFMIN ) THEN
00614                   TEMP = ONE / XMAX
00615                   DO 230 JR = 1, IEND
00616                      VR( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR )
00617   230             CONTINUE
00618                ELSE
00619                   IEND = 0
00620                END IF
00621 *
00622                DO 240 JR = IEND + 1, N
00623                   VR( JR, IEIG ) = CZERO
00624   240          CONTINUE
00625 *
00626             END IF
00627   250    CONTINUE
00628       END IF
00629 *
00630       RETURN
00631 *
00632 *     End of CTGEVC
00633 *
00634       END