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