LAPACK: SRC/cggesx.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CGGESX( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA,
00002      $                   B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR,
00003      $                   LDVSR, RCONDE, RCONDV, WORK, LWORK, RWORK,
00004      $                   IWORK, LIWORK, BWORK, INFO )
00005 *
00006 *  -- LAPACK driver routine (version 3.2) --
00007 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00008 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00009 *     November 2006
00010 *
00011 *     .. Scalar Arguments ..
00012       CHARACTER          JOBVSL, JOBVSR, SENSE, SORT
00013       INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N,
00014      $                   SDIM
00015 *     ..
00016 *     .. Array Arguments ..
00017       LOGICAL            BWORK( * )
00018       INTEGER            IWORK( * )
00019       REAL               RCONDE( 2 ), RCONDV( 2 ), RWORK( * )
00020       COMPLEX            A( LDA, * ), ALPHA( * ), B( LDB, * ),
00021      $                   BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
00022      $                   WORK( * )
00023 *     ..
00024 *     .. Function Arguments ..
00025       LOGICAL            SELCTG
00026       EXTERNAL           SELCTG
00027 *     ..
00028 *
00029 *  Purpose
00030 *  =======
00031 *
00032 *  CGGESX computes for a pair of N-by-N complex nonsymmetric matrices
00033 *  (A,B), the generalized eigenvalues, the complex Schur form (S,T),
00034 *  and, optionally, the left and/or right matrices of Schur vectors (VSL
00035 *  and VSR).  This gives the generalized Schur factorization
00036 *
00037 *       (A,B) = ( (VSL) S (VSR)**H, (VSL) T (VSR)**H )
00038 *
00039 *  where (VSR)**H is the conjugate-transpose of VSR.
00040 *
00041 *  Optionally, it also orders the eigenvalues so that a selected cluster
00042 *  of eigenvalues appears in the leading diagonal blocks of the upper
00043 *  triangular matrix S and the upper triangular matrix T; computes
00044 *  a reciprocal condition number for the average of the selected
00045 *  eigenvalues (RCONDE); and computes a reciprocal condition number for
00046 *  the right and left deflating subspaces corresponding to the selected
00047 *  eigenvalues (RCONDV). The leading columns of VSL and VSR then form
00048 *  an orthonormal basis for the corresponding left and right eigenspaces
00049 *  (deflating subspaces).
00050 *
00051 *  A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
00052 *  or a ratio alpha/beta = w, such that  A - w*B is singular.  It is
00053 *  usually represented as the pair (alpha,beta), as there is a
00054 *  reasonable interpretation for beta=0 or for both being zero.
00055 *
00056 *  A pair of matrices (S,T) is in generalized complex Schur form if T is
00057 *  upper triangular with non-negative diagonal and S is upper
00058 *  triangular.
00059 *
00060 *  Arguments
00061 *  =========
00062 *
00063 *  JOBVSL  (input) CHARACTER*1
00064 *          = 'N':  do not compute the left Schur vectors;
00065 *          = 'V':  compute the left Schur vectors.
00066 *
00067 *  JOBVSR  (input) CHARACTER*1
00068 *          = 'N':  do not compute the right Schur vectors;
00069 *          = 'V':  compute the right Schur vectors.
00070 *
00071 *  SORT    (input) CHARACTER*1
00072 *          Specifies whether or not to order the eigenvalues on the
00073 *          diagonal of the generalized Schur form.
00074 *          = 'N':  Eigenvalues are not ordered;
00075 *          = 'S':  Eigenvalues are ordered (see SELCTG).
00076 *
00077 *  SELCTG  (external procedure) LOGICAL FUNCTION of two COMPLEX arguments
00078 *          SELCTG must be declared EXTERNAL in the calling subroutine.
00079 *          If SORT = 'N', SELCTG is not referenced.
00080 *          If SORT = 'S', SELCTG is used to select eigenvalues to sort
00081 *          to the top left of the Schur form.
00082 *          Note that a selected complex eigenvalue may no longer satisfy
00083 *          SELCTG(ALPHA(j),BETA(j)) = .TRUE. after ordering, since
00084 *          ordering may change the value of complex eigenvalues
00085 *          (especially if the eigenvalue is ill-conditioned), in this
00086 *          case INFO is set to N+3 see INFO below).
00087 *
00088 *  SENSE   (input) CHARACTER*1
00089 *          Determines which reciprocal condition numbers are computed.
00090 *          = 'N' : None are computed;
00091 *          = 'E' : Computed for average of selected eigenvalues only;
00092 *          = 'V' : Computed for selected deflating subspaces only;
00093 *          = 'B' : Computed for both.
00094 *          If SENSE = 'E', 'V', or 'B', SORT must equal 'S'.
00095 *
00096 *  N       (input) INTEGER
00097 *          The order of the matrices A, B, VSL, and VSR.  N >= 0.
00098 *
00099 *  A       (input/output) COMPLEX array, dimension (LDA, N)
00100 *          On entry, the first of the pair of matrices.
00101 *          On exit, A has been overwritten by its generalized Schur
00102 *          form S.
00103 *
00104 *  LDA     (input) INTEGER
00105 *          The leading dimension of A.  LDA >= max(1,N).
00106 *
00107 *  B       (input/output) COMPLEX array, dimension (LDB, N)
00108 *          On entry, the second of the pair of matrices.
00109 *          On exit, B has been overwritten by its generalized Schur
00110 *          form T.
00111 *
00112 *  LDB     (input) INTEGER
00113 *          The leading dimension of B.  LDB >= max(1,N).
00114 *
00115 *  SDIM    (output) INTEGER
00116 *          If SORT = 'N', SDIM = 0.
00117 *          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
00118 *          for which SELCTG is true.
00119 *
00120 *  ALPHA   (output) COMPLEX array, dimension (N)
00121 *  BETA    (output) COMPLEX array, dimension (N)
00122 *          On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
00123 *          generalized eigenvalues.  ALPHA(j) and BETA(j),j=1,...,N  are
00124 *          the diagonals of the complex Schur form (S,T).  BETA(j) will
00125 *          be non-negative real.
00126 *
00127 *          Note: the quotients ALPHA(j)/BETA(j) may easily over- or
00128 *          underflow, and BETA(j) may even be zero.  Thus, the user
00129 *          should avoid naively computing the ratio alpha/beta.
00130 *          However, ALPHA will be always less than and usually
00131 *          comparable with norm(A) in magnitude, and BETA always less
00132 *          than and usually comparable with norm(B).
00133 *
00134 *  VSL     (output) COMPLEX array, dimension (LDVSL,N)
00135 *          If JOBVSL = 'V', VSL will contain the left Schur vectors.
00136 *          Not referenced if JOBVSL = 'N'.
00137 *
00138 *  LDVSL   (input) INTEGER
00139 *          The leading dimension of the matrix VSL. LDVSL >=1, and
00140 *          if JOBVSL = 'V', LDVSL >= N.
00141 *
00142 *  VSR     (output) COMPLEX array, dimension (LDVSR,N)
00143 *          If JOBVSR = 'V', VSR will contain the right Schur vectors.
00144 *          Not referenced if JOBVSR = 'N'.
00145 *
00146 *  LDVSR   (input) INTEGER
00147 *          The leading dimension of the matrix VSR. LDVSR >= 1, and
00148 *          if JOBVSR = 'V', LDVSR >= N.
00149 *
00150 *  RCONDE  (output) REAL array, dimension ( 2 )
00151 *          If SENSE = 'E' or 'B', RCONDE(1) and RCONDE(2) contain the
00152 *          reciprocal condition numbers for the average of the selected
00153 *          eigenvalues.
00154 *          Not referenced if SENSE = 'N' or 'V'.
00155 *
00156 *  RCONDV  (output) REAL array, dimension ( 2 )
00157 *          If SENSE = 'V' or 'B', RCONDV(1) and RCONDV(2) contain the
00158 *          reciprocal condition number for the selected deflating
00159 *          subspaces.
00160 *          Not referenced if SENSE = 'N' or 'E'.
00161 *
00162 *  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
00163 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00164 *
00165 *  LWORK   (input) INTEGER
00166 *          The dimension of the array WORK.
00167 *          If N = 0, LWORK >= 1, else if SENSE = 'E', 'V', or 'B',
00168 *          LWORK >= MAX(1,2*N,2*SDIM*(N-SDIM)), else
00169 *          LWORK >= MAX(1,2*N).  Note that 2*SDIM*(N-SDIM) <= N*N/2.
00170 *          Note also that an error is only returned if
00171 *          LWORK < MAX(1,2*N), but if SENSE = 'E' or 'V' or 'B' this may
00172 *          not be large enough.
00173 *
00174 *          If LWORK = -1, then a workspace query is assumed; the routine
00175 *          only calculates the bound on the optimal size of the WORK
00176 *          array and the minimum size of the IWORK array, returns these
00177 *          values as the first entries of the WORK and IWORK arrays, and
00178 *          no error message related to LWORK or LIWORK is issued by
00179 *          XERBLA.
00180 *
00181 *  RWORK   (workspace) REAL array, dimension ( 8*N )
00182 *          Real workspace.
00183 *
00184 *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
00185 *          On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
00186 *
00187 *  LIWORK  (input) INTEGER
00188 *          The dimension of the array WORK.
00189 *          If SENSE = 'N' or N = 0, LIWORK >= 1, otherwise
00190 *          LIWORK >= N+2.
00191 *
00192 *          If LIWORK = -1, then a workspace query is assumed; the
00193 *          routine only calculates the bound on the optimal size of the
00194 *          WORK array and the minimum size of the IWORK array, returns
00195 *          these values as the first entries of the WORK and IWORK
00196 *          arrays, and no error message related to LWORK or LIWORK is
00197 *          issued by XERBLA.
00198 *
00199 *  BWORK   (workspace) LOGICAL array, dimension (N)
00200 *          Not referenced if SORT = 'N'.
00201 *
00202 *  INFO    (output) INTEGER
00203 *          = 0:  successful exit
00204 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00205 *          = 1,...,N:
00206 *                The QZ iteration failed.  (A,B) are not in Schur
00207 *                form, but ALPHA(j) and BETA(j) should be correct for
00208 *                j=INFO+1,...,N.
00209 *          > N:  =N+1: other than QZ iteration failed in CHGEQZ
00210 *                =N+2: after reordering, roundoff changed values of
00211 *                      some complex eigenvalues so that leading
00212 *                      eigenvalues in the Generalized Schur form no
00213 *                      longer satisfy SELCTG=.TRUE.  This could also
00214 *                      be caused due to scaling.
00215 *                =N+3: reordering failed in CTGSEN.
00216 *
00217 *  =====================================================================
00218 *
00219 *     .. Parameters ..
00220       REAL               ZERO, ONE
00221       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00222       COMPLEX            CZERO, CONE
00223       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
00224      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
00225 *     ..
00226 *     .. Local Scalars ..
00227       LOGICAL            CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
00228      $                   LQUERY, WANTSB, WANTSE, WANTSN, WANTST, WANTSV
00229       INTEGER            I, ICOLS, IERR, IHI, IJOB, IJOBVL, IJOBVR,
00230      $                   ILEFT, ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK,
00231      $                   LIWMIN, LWRK, MAXWRK, MINWRK
00232       REAL               ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PL,
00233      $                   PR, SMLNUM
00234 *     ..
00235 *     .. Local Arrays ..
00236       REAL               DIF( 2 )
00237 *     ..
00238 *     .. External Subroutines ..
00239       EXTERNAL           CGEQRF, CGGBAK, CGGBAL, CGGHRD, CHGEQZ, CLACPY,
00240      $                   CLASCL, CLASET, CTGSEN, CUNGQR, CUNMQR, SLABAD,
00241      $                   XERBLA
00242 *     ..
00243 *     .. External Functions ..
00244       LOGICAL            LSAME
00245       INTEGER            ILAENV
00246       REAL               CLANGE, SLAMCH
00247       EXTERNAL           LSAME, ILAENV, CLANGE, SLAMCH
00248 *     ..
00249 *     .. Intrinsic Functions ..
00250       INTRINSIC          MAX, SQRT
00251 *     ..
00252 *     .. Executable Statements ..
00253 *
00254 *     Decode the input arguments
00255 *
00256       IF( LSAME( JOBVSL, 'N' ) ) THEN
00257          IJOBVL = 1
00258          ILVSL = .FALSE.
00259       ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
00260          IJOBVL = 2
00261          ILVSL = .TRUE.
00262       ELSE
00263          IJOBVL = -1
00264          ILVSL = .FALSE.
00265       END IF
00266 *
00267       IF( LSAME( JOBVSR, 'N' ) ) THEN
00268          IJOBVR = 1
00269          ILVSR = .FALSE.
00270       ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
00271          IJOBVR = 2
00272          ILVSR = .TRUE.
00273       ELSE
00274          IJOBVR = -1
00275          ILVSR = .FALSE.
00276       END IF
00277 *
00278       WANTST = LSAME( SORT, 'S' )
00279       WANTSN = LSAME( SENSE, 'N' )
00280       WANTSE = LSAME( SENSE, 'E' )
00281       WANTSV = LSAME( SENSE, 'V' )
00282       WANTSB = LSAME( SENSE, 'B' )
00283       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
00284       IF( WANTSN ) THEN
00285          IJOB = 0
00286       ELSE IF( WANTSE ) THEN
00287          IJOB = 1
00288       ELSE IF( WANTSV ) THEN
00289          IJOB = 2
00290       ELSE IF( WANTSB ) THEN
00291          IJOB = 4
00292       END IF
00293 *
00294 *     Test the input arguments
00295 *
00296       INFO = 0
00297       IF( IJOBVL.LE.0 ) THEN
00298          INFO = -1
00299       ELSE IF( IJOBVR.LE.0 ) THEN
00300          INFO = -2
00301       ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
00302          INFO = -3
00303       ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR.
00304      $         ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN
00305          INFO = -5
00306       ELSE IF( N.LT.0 ) THEN
00307          INFO = -6
00308       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00309          INFO = -8
00310       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00311          INFO = -10
00312       ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
00313          INFO = -15
00314       ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
00315          INFO = -17
00316       END IF
00317 *
00318 *     Compute workspace
00319 *      (Note: Comments in the code beginning "Workspace:" describe the
00320 *       minimal amount of workspace needed at that point in the code,
00321 *       as well as the preferred amount for good performance.
00322 *       NB refers to the optimal block size for the immediately
00323 *       following subroutine, as returned by ILAENV.)
00324 *
00325       IF( INFO.EQ.0 ) THEN
00326          IF( N.GT.0) THEN
00327             MINWRK = 2*N
00328             MAXWRK = N*(1 + ILAENV( 1, 'CGEQRF', ' ', N, 1, N, 0 ) )
00329             MAXWRK = MAX( MAXWRK, N*( 1 +
00330      $                    ILAENV( 1, 'CUNMQR', ' ', N, 1, N, -1 ) ) )
00331             IF( ILVSL ) THEN
00332                MAXWRK = MAX( MAXWRK, N*( 1 +
00333      $                       ILAENV( 1, 'CUNGQR', ' ', N, 1, N, -1 ) ) )
00334             END IF
00335             LWRK = MAXWRK
00336             IF( IJOB.GE.1 )
00337      $         LWRK = MAX( LWRK, N*N/2 )
00338          ELSE
00339             MINWRK = 1
00340             MAXWRK = 1
00341             LWRK   = 1
00342          END IF
00343          WORK( 1 ) = LWRK
00344          IF( WANTSN .OR. N.EQ.0 ) THEN
00345             LIWMIN = 1
00346          ELSE
00347             LIWMIN = N + 2
00348          END IF
00349          IWORK( 1 ) = LIWMIN
00350 *
00351          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
00352             INFO = -21
00353          ELSE IF( LIWORK.LT.LIWMIN  .AND. .NOT.LQUERY) THEN
00354             INFO = -24
00355          END IF
00356       END IF
00357 *
00358       IF( INFO.NE.0 ) THEN
00359          CALL XERBLA( 'CGGESX', -INFO )
00360          RETURN
00361       ELSE IF (LQUERY) THEN
00362          RETURN
00363       END IF
00364 *
00365 *     Quick return if possible
00366 *
00367       IF( N.EQ.0 ) THEN
00368          SDIM = 0
00369          RETURN
00370       END IF
00371 *
00372 *     Get machine constants
00373 *
00374       EPS = SLAMCH( 'P' )
00375       SMLNUM = SLAMCH( 'S' )
00376       BIGNUM = ONE / SMLNUM
00377       CALL SLABAD( SMLNUM, BIGNUM )
00378       SMLNUM = SQRT( SMLNUM ) / EPS
00379       BIGNUM = ONE / SMLNUM
00380 *
00381 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00382 *
00383       ANRM = CLANGE( 'M', N, N, A, LDA, RWORK )
00384       ILASCL = .FALSE.
00385       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00386          ANRMTO = SMLNUM
00387          ILASCL = .TRUE.
00388       ELSE IF( ANRM.GT.BIGNUM ) THEN
00389          ANRMTO = BIGNUM
00390          ILASCL = .TRUE.
00391       END IF
00392       IF( ILASCL )
00393      $   CALL CLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
00394 *
00395 *     Scale B if max element outside range [SMLNUM,BIGNUM]
00396 *
00397       BNRM = CLANGE( 'M', N, N, B, LDB, RWORK )
00398       ILBSCL = .FALSE.
00399       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
00400          BNRMTO = SMLNUM
00401          ILBSCL = .TRUE.
00402       ELSE IF( BNRM.GT.BIGNUM ) THEN
00403          BNRMTO = BIGNUM
00404          ILBSCL = .TRUE.
00405       END IF
00406       IF( ILBSCL )
00407      $   CALL CLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
00408 *
00409 *     Permute the matrix to make it more nearly triangular
00410 *     (Real Workspace: need 6*N)
00411 *
00412       ILEFT = 1
00413       IRIGHT = N + 1
00414       IRWRK = IRIGHT + N
00415       CALL CGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
00416      $             RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
00417 *
00418 *     Reduce B to triangular form (QR decomposition of B)
00419 *     (Complex Workspace: need N, prefer N*NB)
00420 *
00421       IROWS = IHI + 1 - ILO
00422       ICOLS = N + 1 - ILO
00423       ITAU = 1
00424       IWRK = ITAU + IROWS
00425       CALL CGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
00426      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
00427 *
00428 *     Apply the unitary transformation to matrix A
00429 *     (Complex Workspace: need N, prefer N*NB)
00430 *
00431       CALL CUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
00432      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
00433      $             LWORK+1-IWRK, IERR )
00434 *
00435 *     Initialize VSL
00436 *     (Complex Workspace: need N, prefer N*NB)
00437 *
00438       IF( ILVSL ) THEN
00439          CALL CLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL )
00440          IF( IROWS.GT.1 ) THEN
00441             CALL CLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
00442      $                   VSL( ILO+1, ILO ), LDVSL )
00443          END IF
00444          CALL CUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
00445      $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
00446       END IF
00447 *
00448 *     Initialize VSR
00449 *
00450       IF( ILVSR )
00451      $   CALL CLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR )
00452 *
00453 *     Reduce to generalized Hessenberg form
00454 *     (Workspace: none needed)
00455 *
00456       CALL CGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
00457      $             LDVSL, VSR, LDVSR, IERR )
00458 *
00459       SDIM = 0
00460 *
00461 *     Perform QZ algorithm, computing Schur vectors if desired
00462 *     (Complex Workspace: need N)
00463 *     (Real Workspace:    need N)
00464 *
00465       IWRK = ITAU
00466       CALL CHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
00467      $             ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWRK ),
00468      $             LWORK+1-IWRK, RWORK( IRWRK ), IERR )
00469       IF( IERR.NE.0 ) THEN
00470          IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
00471             INFO = IERR
00472          ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
00473             INFO = IERR - N
00474          ELSE
00475             INFO = N + 1
00476          END IF
00477          GO TO 40
00478       END IF
00479 *
00480 *     Sort eigenvalues ALPHA/BETA and compute the reciprocal of
00481 *     condition number(s)
00482 *
00483       IF( WANTST ) THEN
00484 *
00485 *        Undo scaling on eigenvalues before SELCTGing
00486 *
00487          IF( ILASCL )
00488      $      CALL CLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
00489          IF( ILBSCL )
00490      $      CALL CLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
00491 *
00492 *        Select eigenvalues
00493 *
00494          DO 10 I = 1, N
00495             BWORK( I ) = SELCTG( ALPHA( I ), BETA( I ) )
00496    10    CONTINUE
00497 *
00498 *        Reorder eigenvalues, transform Generalized Schur vectors, and
00499 *        compute reciprocal condition numbers
00500 *        (Complex Workspace: If IJOB >= 1, need MAX(1, 2*SDIM*(N-SDIM))
00501 *                            otherwise, need 1 )
00502 *
00503          CALL CTGSEN( IJOB, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB,
00504      $                ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PL, PR,
00505      $                DIF, WORK( IWRK ), LWORK-IWRK+1, IWORK, LIWORK,
00506      $                IERR )
00507 *
00508          IF( IJOB.GE.1 )
00509      $      MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) )
00510          IF( IERR.EQ.-21 ) THEN
00511 *
00512 *            not enough complex workspace
00513 *
00514             INFO = -21
00515          ELSE
00516             IF( IJOB.EQ.1 .OR. IJOB.EQ.4 ) THEN
00517                RCONDE( 1 ) = PL
00518                RCONDE( 2 ) = PR
00519             END IF
00520             IF( IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
00521                RCONDV( 1 ) = DIF( 1 )
00522                RCONDV( 2 ) = DIF( 2 )
00523             END IF
00524             IF( IERR.EQ.1 )
00525      $         INFO = N + 3
00526          END IF
00527 *
00528       END IF
00529 *
00530 *     Apply permutation to VSL and VSR
00531 *     (Workspace: none needed)
00532 *
00533       IF( ILVSL )
00534      $   CALL CGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
00535      $                RWORK( IRIGHT ), N, VSL, LDVSL, IERR )
00536 *
00537       IF( ILVSR )
00538      $   CALL CGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
00539      $                RWORK( IRIGHT ), N, VSR, LDVSR, IERR )
00540 *
00541 *     Undo scaling
00542 *
00543       IF( ILASCL ) THEN
00544          CALL CLASCL( 'U', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
00545          CALL CLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
00546       END IF
00547 *
00548       IF( ILBSCL ) THEN
00549          CALL CLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
00550          CALL CLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
00551       END IF
00552 *
00553       IF( WANTST ) THEN
00554 *
00555 *        Check if reordering is correct
00556 *
00557          LASTSL = .TRUE.
00558          SDIM = 0
00559          DO 30 I = 1, N
00560             CURSL = SELCTG( ALPHA( I ), BETA( I ) )
00561             IF( CURSL )
00562      $         SDIM = SDIM + 1
00563             IF( CURSL .AND. .NOT.LASTSL )
00564      $         INFO = N + 2
00565             LASTSL = CURSL
00566    30    CONTINUE
00567 *
00568       END IF
00569 *
00570    40 CONTINUE
00571 *
00572       WORK( 1 ) = MAXWRK
00573       IWORK( 1 ) = LIWMIN
00574 *
00575       RETURN
00576 *
00577 *     End of CGGESX
00578 *
00579       END