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