SRC/cgges3.f

   1 *> \brief <b> CGGES3 computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices (blocked algorithm)</b>
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CGGES3 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgges3.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgges3.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgges3.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CGGES3( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B,
  22 *      $                   LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR,
  23 *      $                   WORK, 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 *> CGGES3 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 CGGES3. 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.
 219 *>
 220 *>          If LWORK = -1, then a workspace query is assumed; the routine
 221 *>          only calculates the optimal size of the WORK array, returns
 222 *>          this value as the first entry of the WORK array, and no error
 223 *>          message related to LWORK is issued by XERBLA.
 224 *> \endverbatim
 225 *>
 226 *> \param[out] RWORK
 227 *> \verbatim
 228 *>          RWORK is REAL array, dimension (8*N)
 229 *> \endverbatim
 230 *>
 231 *> \param[out] BWORK
 232 *> \verbatim
 233 *>          BWORK is LOGICAL array, dimension (N)
 234 *>          Not referenced if SORT = 'N'.
 235 *> \endverbatim
 236 *>
 237 *> \param[out] INFO
 238 *> \verbatim
 239 *>          INFO is INTEGER
 240 *>          = 0:  successful exit
 241 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 242 *>          =1,...,N:
 243 *>                The QZ iteration failed.  (A,B) are not in Schur
 244 *>                form, but ALPHA(j) and BETA(j) should be correct for
 245 *>                j=INFO+1,...,N.
 246 *>          > N:  =N+1: other than QZ iteration failed in CHGEQZ
 247 *>                =N+2: after reordering, roundoff changed values of
 248 *>                      some complex eigenvalues so that leading
 249 *>                      eigenvalues in the Generalized Schur form no
 250 *>                      longer satisfy SELCTG=.TRUE.  This could also
 251 *>                      be caused due to scaling.
 252 *>                =N+3: reordering failed in CTGSEN.
 253 *> \endverbatim
 254 *
 255 *  Authors:
 256 *  ========
 257 *
 258 *> \author Univ. of Tennessee
 259 *> \author Univ. of California Berkeley
 260 *> \author Univ. of Colorado Denver
 261 *> \author NAG Ltd.
 262 *
 263 *> \date January 2015
 264 *
 265 *> \ingroup complexGEeigen
 266 *
 267 *  =====================================================================
 268       SUBROUTINE CGGES3( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B,
 269      $                   LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR,
 270      $                   WORK, LWORK, RWORK, BWORK, INFO )
 271 *
 272 *  -- LAPACK driver routine (version 3.6.1) --
 273 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 274 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 275 *     January 2015
 276 *
 277 *     .. Scalar Arguments ..
 278       CHARACTER          JOBVSL, JOBVSR, SORT
 279       INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
 280 *     ..
 281 *     .. Array Arguments ..
 282       LOGICAL            BWORK( * )
 283       REAL               RWORK( * )
 284       COMPLEX            A( LDA, * ), ALPHA( * ), B( LDB, * ),
 285      $                   BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
 286      $                   WORK( * )
 287 *     ..
 288 *     .. Function Arguments ..
 289       LOGICAL            SELCTG
 290       EXTERNAL           SELCTG
 291 *     ..
 292 *
 293 *  =====================================================================
 294 *
 295 *     .. Parameters ..
 296       REAL               ZERO, ONE
 297       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
 298       COMPLEX            CZERO, CONE
 299       PARAMETER          ( CZERO = ( 0.0E0, 0.0E0 ),
 300      $                   CONE = ( 1.0E0, 0.0E0 ) )
 301 *     ..
 302 *     .. Local Scalars ..
 303       LOGICAL            CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
 304      $                   LQUERY, WANTST
 305       INTEGER            I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
 306      $                   ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK, LWKOPT
 307       REAL               ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
 308      $                   PVSR, SMLNUM
 309 *     ..
 310 *     .. Local Arrays ..
 311       INTEGER            IDUM( 1 )
 312       REAL               DIF( 2 )
 313 *     ..
 314 *     .. External Subroutines ..
 315       EXTERNAL           CGEQRF, CGGBAK, CGGBAL, CGGHD3, CHGEQZ, CLACPY,
 316      $                   CLASCL, CLASET, CTGSEN, CUNGQR, CUNMQR, SLABAD,
 317      $                   XERBLA
 318 *     ..
 319 *     .. External Functions ..
 320       LOGICAL            LSAME
 321       REAL               CLANGE, SLAMCH
 322       EXTERNAL           LSAME, CLANGE, SLAMCH
 323 *     ..
 324 *     .. Intrinsic Functions ..
 325       INTRINSIC          MAX, SQRT
 326 *     ..
 327 *     .. Executable Statements ..
 328 *
 329 *     Decode the input arguments
 330 *
 331       IF( LSAME( JOBVSL, 'N' ) ) THEN
 332          IJOBVL = 1
 333          ILVSL = .FALSE.
 334       ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
 335          IJOBVL = 2
 336          ILVSL = .TRUE.
 337       ELSE
 338          IJOBVL = -1
 339          ILVSL = .FALSE.
 340       END IF
 341 *
 342       IF( LSAME( JOBVSR, 'N' ) ) THEN
 343          IJOBVR = 1
 344          ILVSR = .FALSE.
 345       ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
 346          IJOBVR = 2
 347          ILVSR = .TRUE.
 348       ELSE
 349          IJOBVR = -1
 350          ILVSR = .FALSE.
 351       END IF
 352 *
 353       WANTST = LSAME( SORT, 'S' )
 354 *
 355 *     Test the input arguments
 356 *
 357       INFO = 0
 358       LQUERY = ( LWORK.EQ.-1 )
 359       IF( IJOBVL.LE.0 ) THEN
 360          INFO = -1
 361       ELSE IF( IJOBVR.LE.0 ) THEN
 362          INFO = -2
 363       ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
 364          INFO = -3
 365       ELSE IF( N.LT.0 ) THEN
 366          INFO = -5
 367       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 368          INFO = -7
 369       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 370          INFO = -9
 371       ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
 372          INFO = -14
 373       ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
 374          INFO = -16
 375       ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN
 376          INFO = -18
 377       END IF
 378 *
 379 *     Compute workspace
 380 *
 381       IF( INFO.EQ.0 ) THEN
 382          CALL CGEQRF( N, N, B, LDB, WORK, WORK, -1, IERR )
 383          LWKOPT = MAX( 1,  N + INT ( WORK( 1 ) ) )
 384          CALL CUNMQR( 'L', 'C', N, N, N, B, LDB, WORK, A, LDA, WORK,
 385      $                -1, IERR )
 386          LWKOPT = MAX( LWKOPT, N + INT ( WORK( 1 ) ) )
 387          IF( ILVSL ) THEN
 388             CALL CUNGQR( N, N, N, VSL, LDVSL, WORK, WORK, -1,
 389      $                   IERR )
 390             LWKOPT = MAX( LWKOPT, N + INT ( WORK( 1 ) ) )
 391          END IF
 392          CALL CGGHD3( JOBVSL, JOBVSR, N, 1, N, A, LDA, B, LDB, VSL,
 393      $                LDVSL, VSR, LDVSR, WORK, -1, IERR )
 394          LWKOPT = MAX( LWKOPT, N + INT ( WORK( 1 ) ) )
 395          CALL CHGEQZ( 'S', JOBVSL, JOBVSR, N, 1, N, A, LDA, B, LDB,
 396      $                ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK, -1,
 397      $                RWORK, IERR )
 398          LWKOPT = MAX( LWKOPT, INT ( WORK( 1 ) ) )
 399          IF( WANTST ) THEN
 400             CALL CTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB,
 401      $                   ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, SDIM,
 402      $                   PVSL, PVSR, DIF, WORK, -1, IDUM, 1, IERR )
 403             LWKOPT = MAX( LWKOPT, INT ( WORK( 1 ) ) )
 404          END IF
 405          WORK( 1 ) = CMPLX( LWKOPT )
 406       END IF
 407
 408 *
 409       IF( INFO.NE.0 ) THEN
 410          CALL XERBLA( 'CGGES3 ', -INFO )
 411          RETURN
 412       ELSE IF( LQUERY ) THEN
 413          RETURN
 414       END IF
 415 *
 416 *     Quick return if possible
 417 *
 418       IF( N.EQ.0 ) THEN
 419          SDIM = 0
 420          RETURN
 421       END IF
 422 *
 423 *     Get machine constants
 424 *
 425       EPS = SLAMCH( 'P' )
 426       SMLNUM = SLAMCH( 'S' )
 427       BIGNUM = ONE / SMLNUM
 428       CALL SLABAD( SMLNUM, BIGNUM )
 429       SMLNUM = SQRT( SMLNUM ) / EPS
 430       BIGNUM = ONE / SMLNUM
 431 *
 432 *     Scale A if max element outside range [SMLNUM,BIGNUM]
 433 *
 434       ANRM = CLANGE( 'M', N, N, A, LDA, RWORK )
 435       ILASCL = .FALSE.
 436       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 437          ANRMTO = SMLNUM
 438          ILASCL = .TRUE.
 439       ELSE IF( ANRM.GT.BIGNUM ) THEN
 440          ANRMTO = BIGNUM
 441          ILASCL = .TRUE.
 442       END IF
 443 *
 444       IF( ILASCL )
 445      $   CALL CLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
 446 *
 447 *     Scale B if max element outside range [SMLNUM,BIGNUM]
 448 *
 449       BNRM = CLANGE( 'M', N, N, B, LDB, RWORK )
 450       ILBSCL = .FALSE.
 451       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
 452          BNRMTO = SMLNUM
 453          ILBSCL = .TRUE.
 454       ELSE IF( BNRM.GT.BIGNUM ) THEN
 455          BNRMTO = BIGNUM
 456          ILBSCL = .TRUE.
 457       END IF
 458 *
 459       IF( ILBSCL )
 460      $   CALL CLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
 461 *
 462 *     Permute the matrix to make it more nearly triangular
 463 *
 464       ILEFT = 1
 465       IRIGHT = N + 1
 466       IRWRK = IRIGHT + N
 467       CALL CGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
 468      $             RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
 469 *
 470 *     Reduce B to triangular form (QR decomposition of B)
 471 *
 472       IROWS = IHI + 1 - ILO
 473       ICOLS = N + 1 - ILO
 474       ITAU = 1
 475       IWRK = ITAU + IROWS
 476       CALL CGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
 477      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
 478 *
 479 *     Apply the orthogonal transformation to matrix A
 480 *
 481       CALL CUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
 482      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
 483      $             LWORK+1-IWRK, IERR )
 484 *
 485 *     Initialize VSL
 486 *
 487       IF( ILVSL ) THEN
 488          CALL CLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL )
 489          IF( IROWS.GT.1 ) THEN
 490             CALL CLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
 491      $                   VSL( ILO+1, ILO ), LDVSL )
 492          END IF
 493          CALL CUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
 494      $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
 495       END IF
 496 *
 497 *     Initialize VSR
 498 *
 499       IF( ILVSR )
 500      $   CALL CLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR )
 501 *
 502 *     Reduce to generalized Hessenberg form
 503 *
 504       CALL CGGHD3( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
 505      $             LDVSL, VSR, LDVSR, WORK( IWRK ), LWORK+1-IWRK, IERR )
 506 *
 507       SDIM = 0
 508 *
 509 *     Perform QZ algorithm, computing Schur vectors if desired
 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 *
 528       IF( WANTST ) THEN
 529 *
 530 *        Undo scaling on eigenvalues before selecting
 531 *
 532          IF( ILASCL )
 533      $      CALL CLASCL( 'G', 0, 0, ANRM, ANRMTO, N, 1, ALPHA, N, IERR )
 534          IF( ILBSCL )
 535      $      CALL CLASCL( 'G', 0, 0, BNRM, BNRMTO, N, 1, BETA, N, IERR )
 536 *
 537 *        Select eigenvalues
 538 *
 539          DO 10 I = 1, N
 540             BWORK( I ) = SELCTG( ALPHA( I ), BETA( I ) )
 541    10    CONTINUE
 542 *
 543          CALL CTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, ALPHA,
 544      $                BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PVSL, PVSR,
 545      $                DIF, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1, IERR )
 546          IF( IERR.EQ.1 )
 547      $      INFO = N + 3
 548 *
 549       END IF
 550 *
 551 *     Apply back-permutation to VSL and VSR
 552 *
 553       IF( ILVSL )
 554      $   CALL CGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
 555      $                RWORK( IRIGHT ), N, VSL, LDVSL, IERR )
 556       IF( ILVSR )
 557      $   CALL CGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
 558      $                RWORK( IRIGHT ), N, VSR, LDVSR, IERR )
 559 *
 560 *     Undo scaling
 561 *
 562       IF( ILASCL ) THEN
 563          CALL CLASCL( 'U', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
 564          CALL CLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
 565       END IF
 566 *
 567       IF( ILBSCL ) THEN
 568          CALL CLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
 569          CALL CLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
 570       END IF
 571 *
 572       IF( WANTST ) THEN
 573 *
 574 *        Check if reordering is correct
 575 *
 576          LASTSL = .TRUE.
 577          SDIM = 0
 578          DO 20 I = 1, N
 579             CURSL = SELCTG( ALPHA( I ), BETA( I ) )
 580             IF( CURSL )
 581      $         SDIM = SDIM + 1
 582             IF( CURSL .AND. .NOT.LASTSL )
 583      $         INFO = N + 2
 584             LASTSL = CURSL
 585    20    CONTINUE
 586 *
 587       END IF
 588 *
 589    30 CONTINUE
 590 *
 591       WORK( 1 ) = CMPLX( LWKOPT )
 592 *
 593       RETURN
 594 *
 595 *     End of CGGES3
 596 *
 597       END