SRC/cggesx.f

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