SRC/cgges.f

   1 *> \brief <b> CGGES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices</b>
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CGGES + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgges.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgges.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgges.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CGGES( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB,
  22 *                         SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK,
  23 *                         LWORK, RWORK, BWORK, INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          JOBVSL, JOBVSR, SORT
  27 *       INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       LOGICAL            BWORK( * )
  31 *       REAL               RWORK( * )
  32 *       COMPLEX            A( LDA, * ), ALPHA( * ), B( LDB, * ),
  33 *      $                   BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
  34 *      $                   WORK( * )
  35 *       ..
  36 *       .. Function Arguments ..
  37 *       LOGICAL            SELCTG
  38 *       EXTERNAL           SELCTG
  39 *       ..
  40 *
  41 *
  42 *> \par Purpose:
  43 *  =============
  44 *>
  45 *> \verbatim
  46 *>
  47 *> CGGES computes for a pair of N-by-N complex nonsymmetric matrices
  48 *> (A,B), the generalized eigenvalues, the generalized complex Schur
  49 *> form (S, T), and optionally left and/or right Schur vectors (VSL
  50 *> and VSR). This gives the generalized Schur factorization
  51 *>
  52 *>         (A,B) = ( (VSL)*S*(VSR)**H, (VSL)*T*(VSR)**H )
  53 *>
  54 *> where (VSR)**H is the conjugate-transpose of VSR.
  55 *>
  56 *> Optionally, it also orders the eigenvalues so that a selected cluster
  57 *> of eigenvalues appears in the leading diagonal blocks of the upper
  58 *> triangular matrix S and the upper triangular matrix T. The leading
  59 *> columns of VSL and VSR then form an unitary basis for the
  60 *> corresponding left and right eigenspaces (deflating subspaces).
  61 *>
  62 *> (If only the generalized eigenvalues are needed, use the driver
  63 *> CGGEV instead, which is faster.)
  64 *>
  65 *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
  66 *> or a ratio alpha/beta = w, such that  A - w*B is singular.  It is
  67 *> usually represented as the pair (alpha,beta), as there is a
  68 *> reasonable interpretation for beta=0, and even for both being zero.
  69 *>
  70 *> A pair of matrices (S,T) is in generalized complex Schur form if S
  71 *> and T are upper triangular and, in addition, the diagonal elements
  72 *> of T are non-negative real numbers.
  73 *> \endverbatim
  74 *
  75 *  Arguments:
  76 *  ==========
  77 *
  78 *> \param[in] JOBVSL
  79 *> \verbatim
  80 *>          JOBVSL is CHARACTER*1
  81 *>          = 'N':  do not compute the left Schur vectors;
  82 *>          = 'V':  compute the left Schur vectors.
  83 *> \endverbatim
  84 *>
  85 *> \param[in] JOBVSR
  86 *> \verbatim
  87 *>          JOBVSR is CHARACTER*1
  88 *>          = 'N':  do not compute the right Schur vectors;
  89 *>          = 'V':  compute the right Schur vectors.
  90 *> \endverbatim
  91 *>
  92 *> \param[in] SORT
  93 *> \verbatim
  94 *>          SORT is CHARACTER*1
  95 *>          Specifies whether or not to order the eigenvalues on the
  96 *>          diagonal of the generalized Schur form.
  97 *>          = 'N':  Eigenvalues are not ordered;
  98 *>          = 'S':  Eigenvalues are ordered (see SELCTG).
  99 *> \endverbatim
 100 *>
 101 *> \param[in] SELCTG
 102 *> \verbatim
 103 *>          SELCTG is a LOGICAL FUNCTION of two COMPLEX arguments
 104 *>          SELCTG must be declared EXTERNAL in the calling subroutine.
 105 *>          If SORT = 'N', SELCTG is not referenced.
 106 *>          If SORT = 'S', SELCTG is used to select eigenvalues to sort
 107 *>          to the top left of the Schur form.
 108 *>          An eigenvalue ALPHA(j)/BETA(j) is selected if
 109 *>          SELCTG(ALPHA(j),BETA(j)) is true.
 110 *>
 111 *>          Note that a selected complex eigenvalue may no longer satisfy
 112 *>          SELCTG(ALPHA(j),BETA(j)) = .TRUE. after ordering, since
 113 *>          ordering may change the value of complex eigenvalues
 114 *>          (especially if the eigenvalue is ill-conditioned), in this
 115 *>          case INFO is set to N+2 (See INFO below).
 116 *> \endverbatim
 117 *>
 118 *> \param[in] N
 119 *> \verbatim
 120 *>          N is INTEGER
 121 *>          The order of the matrices A, B, VSL, and VSR.  N >= 0.
 122 *> \endverbatim
 123 *>
 124 *> \param[in,out] A
 125 *> \verbatim
 126 *>          A is COMPLEX array, dimension (LDA, N)
 127 *>          On entry, the first of the pair of matrices.
 128 *>          On exit, A has been overwritten by its generalized Schur
 129 *>          form S.
 130 *> \endverbatim
 131 *>
 132 *> \param[in] LDA
 133 *> \verbatim
 134 *>          LDA is INTEGER
 135 *>          The leading dimension of A.  LDA >= max(1,N).
 136 *> \endverbatim
 137 *>
 138 *> \param[in,out] B
 139 *> \verbatim
 140 *>          B is COMPLEX array, dimension (LDB, N)
 141 *>          On entry, the second of the pair of matrices.
 142 *>          On exit, B has been overwritten by its generalized Schur
 143 *>          form T.
 144 *> \endverbatim
 145 *>
 146 *> \param[in] LDB
 147 *> \verbatim
 148 *>          LDB is INTEGER
 149 *>          The leading dimension of B.  LDB >= max(1,N).
 150 *> \endverbatim
 151 *>
 152 *> \param[out] SDIM
 153 *> \verbatim
 154 *>          SDIM is INTEGER
 155 *>          If SORT = 'N', SDIM = 0.
 156 *>          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
 157 *>          for which SELCTG is true.
 158 *> \endverbatim
 159 *>
 160 *> \param[out] ALPHA
 161 *> \verbatim
 162 *>          ALPHA is COMPLEX array, dimension (N)
 163 *> \endverbatim
 164 *>
 165 *> \param[out] BETA
 166 *> \verbatim
 167 *>          BETA is COMPLEX array, dimension (N)
 168 *>          On exit,  ALPHA(j)/BETA(j), j=1,...,N, will be the
 169 *>          generalized eigenvalues.  ALPHA(j), j=1,...,N  and  BETA(j),
 170 *>          j=1,...,N  are the diagonals of the complex Schur form (A,B)
 171 *>          output by CGGES. The  BETA(j) will be non-negative real.
 172 *>
 173 *>          Note: the quotients ALPHA(j)/BETA(j) may easily over- or
 174 *>          underflow, and BETA(j) may even be zero.  Thus, the user
 175 *>          should avoid naively computing the ratio alpha/beta.
 176 *>          However, ALPHA will be always less than and usually
 177 *>          comparable with norm(A) in magnitude, and BETA always less
 178 *>          than and usually comparable with norm(B).
 179 *> \endverbatim
 180 *>
 181 *> \param[out] VSL
 182 *> \verbatim
 183 *>          VSL is COMPLEX array, dimension (LDVSL,N)
 184 *>          If JOBVSL = 'V', VSL will contain the left Schur vectors.
 185 *>          Not referenced if JOBVSL = 'N'.
 186 *> \endverbatim
 187 *>
 188 *> \param[in] LDVSL
 189 *> \verbatim
 190 *>          LDVSL is INTEGER
 191 *>          The leading dimension of the matrix VSL. LDVSL >= 1, and
 192 *>          if JOBVSL = 'V', LDVSL >= N.
 193 *> \endverbatim
 194 *>
 195 *> \param[out] VSR
 196 *> \verbatim
 197 *>          VSR is COMPLEX array, dimension (LDVSR,N)
 198 *>          If JOBVSR = 'V', VSR will contain the right Schur vectors.
 199 *>          Not referenced if JOBVSR = 'N'.
 200 *> \endverbatim
 201 *>
 202 *> \param[in] LDVSR
 203 *> \verbatim
 204 *>          LDVSR is INTEGER
 205 *>          The leading dimension of the matrix VSR. LDVSR >= 1, and
 206 *>          if JOBVSR = 'V', LDVSR >= N.
 207 *> \endverbatim
 208 *>
 209 *> \param[out] WORK
 210 *> \verbatim
 211 *>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
 212 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 213 *> \endverbatim
 214 *>
 215 *> \param[in] LWORK
 216 *> \verbatim
 217 *>          LWORK is INTEGER
 218 *>          The dimension of the array WORK.  LWORK >= max(1,2*N).
 219 *>          For good performance, LWORK must generally be larger.
 220 *>
 221 *>          If LWORK = -1, then a workspace query is assumed; the routine
 222 *>          only calculates the optimal size of the WORK array, returns
 223 *>          this value as the first entry of the WORK array, and no error
 224 *>          message related to LWORK is issued by XERBLA.
 225 *> \endverbatim
 226 *>
 227 *> \param[out] RWORK
 228 *> \verbatim
 229 *>          RWORK is REAL array, dimension (8*N)
 230 *> \endverbatim
 231 *>
 232 *> \param[out] BWORK
 233 *> \verbatim
 234 *>          BWORK is LOGICAL array, dimension (N)
 235 *>          Not referenced if SORT = 'N'.
 236 *> \endverbatim
 237 *>
 238 *> \param[out] INFO
 239 *> \verbatim
 240 *>          INFO is INTEGER
 241 *>          = 0:  successful exit
 242 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 243 *>          =1,...,N:
 244 *>                The QZ iteration failed.  (A,B) are not in Schur
 245 *>                form, but ALPHA(j) and BETA(j) should be correct for
 246 *>                j=INFO+1,...,N.
 247 *>          > N:  =N+1: other than QZ iteration failed in CHGEQZ
 248 *>                =N+2: after reordering, roundoff changed values of
 249 *>                      some complex eigenvalues so that leading
 250 *>                      eigenvalues in the Generalized Schur form no
 251 *>                      longer satisfy SELCTG=.TRUE.  This could also
 252 *>                      be caused due to scaling.
 253 *>                =N+3: reordering failed in CTGSEN.
 254 *> \endverbatim
 255 *
 256 *  Authors:
 257 *  ========
 258 *
 259 *> \author Univ. of Tennessee
 260 *> \author Univ. of California Berkeley
 261 *> \author Univ. of Colorado Denver
 262 *> \author NAG Ltd.
 263 *
 264 *> \date November 2015
 265 *
 266 *> \ingroup complexGEeigen
 267 *
 268 *  =====================================================================
 269       SUBROUTINE CGGES( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB,
 270      $                  SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK,
 271      $                  LWORK, RWORK, BWORK, INFO )
 272 *
 273 *  -- LAPACK driver routine (version 3.6.0) --
 274 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 275 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 276 *     November 2015
 277 *
 278 *     .. Scalar Arguments ..
 279       CHARACTER          JOBVSL, JOBVSR, SORT
 280       INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
 281 *     ..
 282 *     .. Array Arguments ..
 283       LOGICAL            BWORK( * )
 284       REAL               RWORK( * )
 285       COMPLEX            A( LDA, * ), ALPHA( * ), B( LDB, * ),
 286      $                   BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
 287      $                   WORK( * )
 288 *     ..
 289 *     .. Function Arguments ..
 290       LOGICAL            SELCTG
 291       EXTERNAL           SELCTG
 292 *     ..
 293 *
 294 *  =====================================================================
 295 *
 296 *     .. Parameters ..
 297       REAL               ZERO, ONE
 298       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
 299       COMPLEX            CZERO, CONE
 300       PARAMETER          ( CZERO = ( 0.0E0, 0.0E0 ),
 301      $                   CONE = ( 1.0E0, 0.0E0 ) )
 302 *     ..
 303 *     .. Local Scalars ..
 304       LOGICAL            CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
 305      $                   LQUERY, WANTST
 306       INTEGER            I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
 307      $                   ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK, LWKMIN,
 308      $                   LWKOPT
 309       REAL               ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
 310      $                   PVSR, SMLNUM
 311 *     ..
 312 *     .. Local Arrays ..
 313       INTEGER            IDUM( 1 )
 314       REAL               DIF( 2 )
 315 *     ..
 316 *     .. External Subroutines ..
 317       EXTERNAL           CGEQRF, CGGBAK, CGGBAL, CGGHRD, CHGEQZ, CLACPY,
 318      $                   CLASCL, CLASET, CTGSEN, CUNGQR, CUNMQR, SLABAD,
 319      $                   XERBLA
 320 *     ..
 321 *     .. External Functions ..
 322       LOGICAL            LSAME
 323       INTEGER            ILAENV
 324       REAL               CLANGE, SLAMCH
 325       EXTERNAL           LSAME, ILAENV, CLANGE, SLAMCH
 326 *     ..
 327 *     .. Intrinsic Functions ..
 328       INTRINSIC          MAX, SQRT
 329 *     ..
 330 *     .. Executable Statements ..
 331 *
 332 *     Decode the input arguments
 333 *
 334       IF( LSAME( JOBVSL, 'N' ) ) THEN
 335          IJOBVL = 1
 336          ILVSL = .FALSE.
 337       ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
 338          IJOBVL = 2
 339          ILVSL = .TRUE.
 340       ELSE
 341          IJOBVL = -1
 342          ILVSL = .FALSE.
 343       END IF
 344 *
 345       IF( LSAME( JOBVSR, 'N' ) ) THEN
 346          IJOBVR = 1
 347          ILVSR = .FALSE.
 348       ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
 349          IJOBVR = 2
 350          ILVSR = .TRUE.
 351       ELSE
 352          IJOBVR = -1
 353          ILVSR = .FALSE.
 354       END IF
 355 *
 356       WANTST = LSAME( SORT, 'S' )
 357 *
 358 *     Test the input arguments
 359 *
 360       INFO = 0
 361       LQUERY = ( LWORK.EQ.-1 )
 362       IF( IJOBVL.LE.0 ) THEN
 363          INFO = -1
 364       ELSE IF( IJOBVR.LE.0 ) THEN
 365          INFO = -2
 366       ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
 367          INFO = -3
 368       ELSE IF( N.LT.0 ) THEN
 369          INFO = -5
 370       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 371          INFO = -7
 372       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 373          INFO = -9
 374       ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
 375          INFO = -14
 376       ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
 377          INFO = -16
 378       END IF
 379 *
 380 *     Compute workspace
 381 *      (Note: Comments in the code beginning "Workspace:" describe the
 382 *       minimal amount of workspace needed at that point in the code,
 383 *       as well as the preferred amount for good performance.
 384 *       NB refers to the optimal block size for the immediately
 385 *       following subroutine, as returned by ILAENV.)
 386 *
 387       IF( INFO.EQ.0 ) THEN
 388          LWKMIN = MAX( 1, 2*N )
 389          LWKOPT = MAX( 1, N + N*ILAENV( 1, 'CGEQRF', ' ', N, 1, N, 0 ) )
 390          LWKOPT = MAX( LWKOPT, N +
 391      $                 N*ILAENV( 1, 'CUNMQR', ' ', N, 1, N, -1 ) )
 392          IF( ILVSL ) THEN
 393             LWKOPT = MAX( LWKOPT, N +
 394      $                    N*ILAENV( 1, 'CUNGQR', ' ', N, 1, N, -1 ) )
 395          END IF
 396          WORK( 1 ) = LWKOPT
 397 *
 398          IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY )
 399      $      INFO = -18
 400       END IF
 401 *
 402       IF( INFO.NE.0 ) THEN
 403          CALL XERBLA( 'CGGES ', -INFO )
 404          RETURN
 405       ELSE IF( LQUERY ) THEN
 406          RETURN
 407       END IF
 408 *
 409 *     Quick return if possible
 410 *
 411       IF( N.EQ.0 ) THEN
 412          SDIM = 0
 413          RETURN
 414       END IF
 415 *
 416 *     Get machine constants
 417 *
 418       EPS = SLAMCH( 'P' )
 419       SMLNUM = SLAMCH( 'S' )
 420       BIGNUM = ONE / SMLNUM
 421       CALL SLABAD( SMLNUM, BIGNUM )
 422       SMLNUM = SQRT( SMLNUM ) / EPS
 423       BIGNUM = ONE / SMLNUM
 424 *
 425 *     Scale A if max element outside range [SMLNUM,BIGNUM]
 426 *
 427       ANRM = CLANGE( 'M', N, N, A, LDA, RWORK )
 428       ILASCL = .FALSE.
 429       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 430          ANRMTO = SMLNUM
 431          ILASCL = .TRUE.
 432       ELSE IF( ANRM.GT.BIGNUM ) THEN
 433          ANRMTO = BIGNUM
 434          ILASCL = .TRUE.
 435       END IF
 436 *
 437       IF( ILASCL )
 438      $   CALL CLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
 439 *
 440 *     Scale B if max element outside range [SMLNUM,BIGNUM]
 441 *
 442       BNRM = CLANGE( 'M', N, N, B, LDB, RWORK )
 443       ILBSCL = .FALSE.
 444       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
 445          BNRMTO = SMLNUM
 446          ILBSCL = .TRUE.
 447       ELSE IF( BNRM.GT.BIGNUM ) THEN
 448          BNRMTO = BIGNUM
 449          ILBSCL = .TRUE.
 450       END IF
 451 *
 452       IF( ILBSCL )
 453      $   CALL CLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
 454 *
 455 *     Permute the matrix to make it more nearly triangular
 456 *     (Real Workspace: need 6*N)
 457 *
 458       ILEFT = 1
 459       IRIGHT = N + 1
 460       IRWRK = IRIGHT + N
 461       CALL CGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
 462      $             RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
 463 *
 464 *     Reduce B to triangular form (QR decomposition of B)
 465 *     (Complex Workspace: need N, prefer N*NB)
 466 *
 467       IROWS = IHI + 1 - ILO
 468       ICOLS = N + 1 - ILO
 469       ITAU = 1
 470       IWRK = ITAU + IROWS
 471       CALL CGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
 472      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
 473 *
 474 *     Apply the orthogonal transformation to matrix A
 475 *     (Complex Workspace: need N, prefer N*NB)
 476 *
 477       CALL CUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
 478      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
 479      $             LWORK+1-IWRK, IERR )
 480 *
 481 *     Initialize VSL
 482 *     (Complex Workspace: need N, prefer N*NB)
 483 *
 484       IF( ILVSL ) THEN
 485          CALL CLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL )
 486          IF( IROWS.GT.1 ) THEN
 487             CALL CLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
 488      $                   VSL( ILO+1, ILO ), LDVSL )
 489          END IF
 490          CALL CUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
 491      $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
 492       END IF
 493 *
 494 *     Initialize VSR
 495 *
 496       IF( ILVSR )
 497      $   CALL CLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR )
 498 *
 499 *     Reduce to generalized Hessenberg form
 500 *     (Workspace: none needed)
 501 *
 502       CALL CGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
 503      $             LDVSL, VSR, LDVSR, IERR )
 504 *
 505       SDIM = 0
 506 *
 507 *     Perform QZ algorithm, computing Schur vectors if desired
 508 *     (Complex Workspace: need N)
 509 *     (Real Workspace: need N)
 510 *
 511       IWRK = ITAU
 512       CALL CHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
 513      $             ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWRK ),
 514      $             LWORK+1-IWRK, RWORK( IRWRK ), IERR )
 515       IF( IERR.NE.0 ) THEN
 516          IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
 517             INFO = IERR
 518          ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
 519             INFO = IERR - N
 520          ELSE
 521             INFO = N + 1
 522          END IF
 523          GO TO 30
 524       END IF
 525 *
 526 *     Sort eigenvalues ALPHA/BETA if desired
 527 *     (Workspace: none needed)
 528 *
 529       IF( WANTST ) THEN
 530 *
 531 *        Undo scaling on eigenvalues before selecting
 532 *
 533          IF( ILASCL )
 534      $      CALL CLASCL( 'G', 0, 0, ANRM, ANRMTO, N, 1, ALPHA, N, IERR )
 535          IF( ILBSCL )
 536      $      CALL CLASCL( 'G', 0, 0, BNRM, BNRMTO, N, 1, BETA, N, IERR )
 537 *
 538 *        Select eigenvalues
 539 *
 540          DO 10 I = 1, N
 541             BWORK( I ) = SELCTG( ALPHA( I ), BETA( I ) )
 542    10    CONTINUE
 543 *
 544          CALL CTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, ALPHA,
 545      $                BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PVSL, PVSR,
 546      $                DIF, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1, IERR )
 547          IF( IERR.EQ.1 )
 548      $      INFO = N + 3
 549 *
 550       END IF
 551 *
 552 *     Apply back-permutation to VSL and VSR
 553 *     (Workspace: none needed)
 554 *
 555       IF( ILVSL )
 556      $   CALL CGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
 557      $                RWORK( IRIGHT ), N, VSL, LDVSL, IERR )
 558       IF( ILVSR )
 559      $   CALL CGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
 560      $                RWORK( IRIGHT ), N, VSR, LDVSR, IERR )
 561 *
 562 *     Undo scaling
 563 *
 564       IF( ILASCL ) THEN
 565          CALL CLASCL( 'U', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
 566          CALL CLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
 567       END IF
 568 *
 569       IF( ILBSCL ) THEN
 570          CALL CLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
 571          CALL CLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
 572       END IF
 573 *
 574       IF( WANTST ) THEN
 575 *
 576 *        Check if reordering is correct
 577 *
 578          LASTSL = .TRUE.
 579          SDIM = 0
 580          DO 20 I = 1, N
 581             CURSL = SELCTG( ALPHA( I ), BETA( I ) )
 582             IF( CURSL )
 583      $         SDIM = SDIM + 1
 584             IF( CURSL .AND. .NOT.LASTSL )
 585      $         INFO = N + 2
 586             LASTSL = CURSL
 587    20    CONTINUE
 588 *
 589       END IF
 590 *
 591    30 CONTINUE
 592 *
 593       WORK( 1 ) = LWKOPT
 594 *
 595       RETURN
 596 *
 597 *     End of CGGES
 598 *
 599       END