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