LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sggevx.f
Go to the documentation of this file.
00001 *> \brief <b> SGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors 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 SGGEVX + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sggevx.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sggevx.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sggevx.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
00022 *                          ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
00023 *                          IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
00024 *                          RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
00025 * 
00026 *       .. Scalar Arguments ..
00027 *       CHARACTER          BALANC, JOBVL, JOBVR, SENSE
00028 *       INTEGER            IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
00029 *       REAL               ABNRM, BBNRM
00030 *       ..
00031 *       .. Array Arguments ..
00032 *       LOGICAL            BWORK( * )
00033 *       INTEGER            IWORK( * )
00034 *       REAL               A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
00035 *      $                   B( LDB, * ), BETA( * ), LSCALE( * ),
00036 *      $                   RCONDE( * ), RCONDV( * ), RSCALE( * ),
00037 *      $                   VL( LDVL, * ), VR( LDVR, * ), WORK( * )
00038 *       ..
00039 *  
00040 *
00041 *> \par Purpose:
00042 *  =============
00043 *>
00044 *> \verbatim
00045 *>
00046 *> SGGEVX computes for a pair of N-by-N real nonsymmetric matrices (A,B)
00047 *> the generalized eigenvalues, and optionally, the left and/or right
00048 *> generalized eigenvectors.
00049 *>
00050 *> Optionally also, it computes a balancing transformation to improve
00051 *> the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
00052 *> LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for
00053 *> the eigenvalues (RCONDE), and reciprocal condition numbers for the
00054 *> right eigenvectors (RCONDV).
00055 *>
00056 *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar
00057 *> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
00058 *> singular. It is usually represented as the pair (alpha,beta), as
00059 *> there is a reasonable interpretation for beta=0, and even for both
00060 *> being zero.
00061 *>
00062 *> The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
00063 *> of (A,B) satisfies
00064 *>
00065 *>                  A * v(j) = lambda(j) * B * v(j) .
00066 *>
00067 *> The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
00068 *> of (A,B) satisfies
00069 *>
00070 *>                  u(j)**H * A  = lambda(j) * u(j)**H * B.
00071 *>
00072 *> where u(j)**H is the conjugate-transpose of u(j).
00073 *>
00074 *> \endverbatim
00075 *
00076 *  Arguments:
00077 *  ==========
00078 *
00079 *> \param[in] BALANC
00080 *> \verbatim
00081 *>          BALANC is CHARACTER*1
00082 *>          Specifies the balance option to be performed.
00083 *>          = 'N':  do not diagonally scale or permute;
00084 *>          = 'P':  permute only;
00085 *>          = 'S':  scale only;
00086 *>          = 'B':  both permute and scale.
00087 *>          Computed reciprocal condition numbers will be for the
00088 *>          matrices after permuting and/or balancing. Permuting does
00089 *>          not change condition numbers (in exact arithmetic), but
00090 *>          balancing does.
00091 *> \endverbatim
00092 *>
00093 *> \param[in] JOBVL
00094 *> \verbatim
00095 *>          JOBVL is CHARACTER*1
00096 *>          = 'N':  do not compute the left generalized eigenvectors;
00097 *>          = 'V':  compute the left generalized eigenvectors.
00098 *> \endverbatim
00099 *>
00100 *> \param[in] JOBVR
00101 *> \verbatim
00102 *>          JOBVR is CHARACTER*1
00103 *>          = 'N':  do not compute the right generalized eigenvectors;
00104 *>          = 'V':  compute the right generalized eigenvectors.
00105 *> \endverbatim
00106 *>
00107 *> \param[in] SENSE
00108 *> \verbatim
00109 *>          SENSE is CHARACTER*1
00110 *>          Determines which reciprocal condition numbers are computed.
00111 *>          = 'N': none are computed;
00112 *>          = 'E': computed for eigenvalues only;
00113 *>          = 'V': computed for eigenvectors only;
00114 *>          = 'B': computed for eigenvalues and eigenvectors.
00115 *> \endverbatim
00116 *>
00117 *> \param[in] N
00118 *> \verbatim
00119 *>          N is INTEGER
00120 *>          The order of the matrices A, B, VL, and VR.  N >= 0.
00121 *> \endverbatim
00122 *>
00123 *> \param[in,out] A
00124 *> \verbatim
00125 *>          A is REAL array, dimension (LDA, N)
00126 *>          On entry, the matrix A in the pair (A,B).
00127 *>          On exit, A has been overwritten. If JOBVL='V' or JOBVR='V'
00128 *>          or both, then A contains the first part of the real Schur
00129 *>          form of the "balanced" versions of the input A and B.
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 REAL array, dimension (LDB, N)
00141 *>          On entry, the matrix B in the pair (A,B).
00142 *>          On exit, B has been overwritten. If JOBVL='V' or JOBVR='V'
00143 *>          or both, then B contains the second part of the real Schur
00144 *>          form of the "balanced" versions of the input A and B.
00145 *> \endverbatim
00146 *>
00147 *> \param[in] LDB
00148 *> \verbatim
00149 *>          LDB is INTEGER
00150 *>          The leading dimension of B.  LDB >= max(1,N).
00151 *> \endverbatim
00152 *>
00153 *> \param[out] ALPHAR
00154 *> \verbatim
00155 *>          ALPHAR is REAL array, dimension (N)
00156 *> \endverbatim
00157 *>
00158 *> \param[out] ALPHAI
00159 *> \verbatim
00160 *>          ALPHAI is REAL array, dimension (N)
00161 *> \endverbatim
00162 *>
00163 *> \param[out] BETA
00164 *> \verbatim
00165 *>          BETA is REAL array, dimension (N)
00166 *>          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
00167 *>          be the generalized eigenvalues.  If ALPHAI(j) is zero, then
00168 *>          the j-th eigenvalue is real; if positive, then the j-th and
00169 *>          (j+1)-st eigenvalues are a complex conjugate pair, with
00170 *>          ALPHAI(j+1) negative.
00171 *>
00172 *>          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
00173 *>          may easily over- or underflow, and BETA(j) may even be zero.
00174 *>          Thus, the user should avoid naively computing the ratio
00175 *>          ALPHA/BETA. However, ALPHAR and ALPHAI will be always less
00176 *>          than and usually comparable with norm(A) in magnitude, and
00177 *>          BETA always less than and usually comparable with norm(B).
00178 *> \endverbatim
00179 *>
00180 *> \param[out] VL
00181 *> \verbatim
00182 *>          VL is REAL array, dimension (LDVL,N)
00183 *>          If JOBVL = 'V', the left eigenvectors u(j) are stored one
00184 *>          after another in the columns of VL, in the same order as
00185 *>          their eigenvalues. If the j-th eigenvalue is real, then
00186 *>          u(j) = VL(:,j), the j-th column of VL. If the j-th and
00187 *>          (j+1)-th eigenvalues form a complex conjugate pair, then
00188 *>          u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1).
00189 *>          Each eigenvector will be scaled so the largest component have
00190 *>          abs(real part) + abs(imag. part) = 1.
00191 *>          Not referenced if JOBVL = 'N'.
00192 *> \endverbatim
00193 *>
00194 *> \param[in] LDVL
00195 *> \verbatim
00196 *>          LDVL is INTEGER
00197 *>          The leading dimension of the matrix VL. LDVL >= 1, and
00198 *>          if JOBVL = 'V', LDVL >= N.
00199 *> \endverbatim
00200 *>
00201 *> \param[out] VR
00202 *> \verbatim
00203 *>          VR is REAL array, dimension (LDVR,N)
00204 *>          If JOBVR = 'V', the right eigenvectors v(j) are stored one
00205 *>          after another in the columns of VR, in the same order as
00206 *>          their eigenvalues. If the j-th eigenvalue is real, then
00207 *>          v(j) = VR(:,j), the j-th column of VR. If the j-th and
00208 *>          (j+1)-th eigenvalues form a complex conjugate pair, then
00209 *>          v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
00210 *>          Each eigenvector will be scaled so the largest component have
00211 *>          abs(real part) + abs(imag. part) = 1.
00212 *>          Not referenced if JOBVR = 'N'.
00213 *> \endverbatim
00214 *>
00215 *> \param[in] LDVR
00216 *> \verbatim
00217 *>          LDVR is INTEGER
00218 *>          The leading dimension of the matrix VR. LDVR >= 1, and
00219 *>          if JOBVR = 'V', LDVR >= N.
00220 *> \endverbatim
00221 *>
00222 *> \param[out] ILO
00223 *> \verbatim
00224 *>          ILO is INTEGER
00225 *> \endverbatim
00226 *>
00227 *> \param[out] IHI
00228 *> \verbatim
00229 *>          IHI is INTEGER
00230 *>          ILO and IHI are integer values such that on exit
00231 *>          A(i,j) = 0 and B(i,j) = 0 if i > j and
00232 *>          j = 1,...,ILO-1 or i = IHI+1,...,N.
00233 *>          If BALANC = 'N' or 'S', ILO = 1 and IHI = N.
00234 *> \endverbatim
00235 *>
00236 *> \param[out] LSCALE
00237 *> \verbatim
00238 *>          LSCALE is REAL array, dimension (N)
00239 *>          Details of the permutations and scaling factors applied
00240 *>          to the left side of A and B.  If PL(j) is the index of the
00241 *>          row interchanged with row j, and DL(j) is the scaling
00242 *>          factor applied to row j, then
00243 *>            LSCALE(j) = PL(j)  for j = 1,...,ILO-1
00244 *>                      = DL(j)  for j = ILO,...,IHI
00245 *>                      = PL(j)  for j = IHI+1,...,N.
00246 *>          The order in which the interchanges are made is N to IHI+1,
00247 *>          then 1 to ILO-1.
00248 *> \endverbatim
00249 *>
00250 *> \param[out] RSCALE
00251 *> \verbatim
00252 *>          RSCALE is REAL array, dimension (N)
00253 *>          Details of the permutations and scaling factors applied
00254 *>          to the right side of A and B.  If PR(j) is the index of the
00255 *>          column interchanged with column j, and DR(j) is the scaling
00256 *>          factor applied to column j, then
00257 *>            RSCALE(j) = PR(j)  for j = 1,...,ILO-1
00258 *>                      = DR(j)  for j = ILO,...,IHI
00259 *>                      = PR(j)  for j = IHI+1,...,N
00260 *>          The order in which the interchanges are made is N to IHI+1,
00261 *>          then 1 to ILO-1.
00262 *> \endverbatim
00263 *>
00264 *> \param[out] ABNRM
00265 *> \verbatim
00266 *>          ABNRM is REAL
00267 *>          The one-norm of the balanced matrix A.
00268 *> \endverbatim
00269 *>
00270 *> \param[out] BBNRM
00271 *> \verbatim
00272 *>          BBNRM is REAL
00273 *>          The one-norm of the balanced matrix B.
00274 *> \endverbatim
00275 *>
00276 *> \param[out] RCONDE
00277 *> \verbatim
00278 *>          RCONDE is REAL array, dimension (N)
00279 *>          If SENSE = 'E' or 'B', the reciprocal condition numbers of
00280 *>          the eigenvalues, stored in consecutive elements of the array.
00281 *>          For a complex conjugate pair of eigenvalues two consecutive
00282 *>          elements of RCONDE are set to the same value. Thus RCONDE(j),
00283 *>          RCONDV(j), and the j-th columns of VL and VR all correspond
00284 *>          to the j-th eigenpair.
00285 *>          If SENSE = 'N' or 'V', RCONDE is not referenced.
00286 *> \endverbatim
00287 *>
00288 *> \param[out] RCONDV
00289 *> \verbatim
00290 *>          RCONDV is REAL array, dimension (N)
00291 *>          If SENSE = 'V' or 'B', the estimated reciprocal condition
00292 *>          numbers of the eigenvectors, stored in consecutive elements
00293 *>          of the array. For a complex eigenvector two consecutive
00294 *>          elements of RCONDV are set to the same value. If the
00295 *>          eigenvalues cannot be reordered to compute RCONDV(j),
00296 *>          RCONDV(j) is set to 0; this can only occur when the true
00297 *>          value would be very small anyway.
00298 *>          If SENSE = 'N' or 'E', RCONDV is not referenced.
00299 *> \endverbatim
00300 *>
00301 *> \param[out] WORK
00302 *> \verbatim
00303 *>          WORK is REAL array, dimension (MAX(1,LWORK))
00304 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00305 *> \endverbatim
00306 *>
00307 *> \param[in] LWORK
00308 *> \verbatim
00309 *>          LWORK is INTEGER
00310 *>          The dimension of the array WORK. LWORK >= max(1,2*N).
00311 *>          If BALANC = 'S' or 'B', or JOBVL = 'V', or JOBVR = 'V',
00312 *>          LWORK >= max(1,6*N).
00313 *>          If SENSE = 'E', LWORK >= max(1,10*N).
00314 *>          If SENSE = 'V' or 'B', LWORK >= 2*N*N+8*N+16.
00315 *>
00316 *>          If LWORK = -1, then a workspace query is assumed; the routine
00317 *>          only calculates the optimal size of the WORK array, returns
00318 *>          this value as the first entry of the WORK array, and no error
00319 *>          message related to LWORK is issued by XERBLA.
00320 *> \endverbatim
00321 *>
00322 *> \param[out] IWORK
00323 *> \verbatim
00324 *>          IWORK is INTEGER array, dimension (N+6)
00325 *>          If SENSE = 'E', IWORK is not referenced.
00326 *> \endverbatim
00327 *>
00328 *> \param[out] BWORK
00329 *> \verbatim
00330 *>          BWORK is LOGICAL array, dimension (N)
00331 *>          If SENSE = 'N', BWORK is not referenced.
00332 *> \endverbatim
00333 *>
00334 *> \param[out] INFO
00335 *> \verbatim
00336 *>          INFO is INTEGER
00337 *>          = 0:  successful exit
00338 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00339 *>          = 1,...,N:
00340 *>                The QZ iteration failed.  No eigenvectors have been
00341 *>                calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
00342 *>                should be correct for j=INFO+1,...,N.
00343 *>          > N:  =N+1: other than QZ iteration failed in SHGEQZ.
00344 *>                =N+2: error return from STGEVC.
00345 *> \endverbatim
00346 *
00347 *  Authors:
00348 *  ========
00349 *
00350 *> \author Univ. of Tennessee 
00351 *> \author Univ. of California Berkeley 
00352 *> \author Univ. of Colorado Denver 
00353 *> \author NAG Ltd. 
00354 *
00355 *> \date April 2012
00356 *
00357 *> \ingroup realGEeigen
00358 *
00359 *> \par Further Details:
00360 *  =====================
00361 *>
00362 *> \verbatim
00363 *>
00364 *>  Balancing a matrix pair (A,B) includes, first, permuting rows and
00365 *>  columns to isolate eigenvalues, second, applying diagonal similarity
00366 *>  transformation to the rows and columns to make the rows and columns
00367 *>  as close in norm as possible. The computed reciprocal condition
00368 *>  numbers correspond to the balanced matrix. Permuting rows and columns
00369 *>  will not change the condition numbers (in exact arithmetic) but
00370 *>  diagonal scaling will.  For further explanation of balancing, see
00371 *>  section 4.11.1.2 of LAPACK Users' Guide.
00372 *>
00373 *>  An approximate error bound on the chordal distance between the i-th
00374 *>  computed generalized eigenvalue w and the corresponding exact
00375 *>  eigenvalue lambda is
00376 *>
00377 *>       chord(w, lambda) <= EPS * norm(ABNRM, BBNRM) / RCONDE(I)
00378 *>
00379 *>  An approximate error bound for the angle between the i-th computed
00380 *>  eigenvector VL(i) or VR(i) is given by
00381 *>
00382 *>       EPS * norm(ABNRM, BBNRM) / DIF(i).
00383 *>
00384 *>  For further explanation of the reciprocal condition numbers RCONDE
00385 *>  and RCONDV, see section 4.11 of LAPACK User's Guide.
00386 *> \endverbatim
00387 *>
00388 *  =====================================================================
00389       SUBROUTINE SGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
00390      $                   ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
00391      $                   IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
00392      $                   RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
00393 *
00394 *  -- LAPACK driver routine (version 3.4.1) --
00395 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00396 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00397 *     April 2012
00398 *
00399 *     .. Scalar Arguments ..
00400       CHARACTER          BALANC, JOBVL, JOBVR, SENSE
00401       INTEGER            IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
00402       REAL               ABNRM, BBNRM
00403 *     ..
00404 *     .. Array Arguments ..
00405       LOGICAL            BWORK( * )
00406       INTEGER            IWORK( * )
00407       REAL               A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
00408      $                   B( LDB, * ), BETA( * ), LSCALE( * ),
00409      $                   RCONDE( * ), RCONDV( * ), RSCALE( * ),
00410      $                   VL( LDVL, * ), VR( LDVR, * ), WORK( * )
00411 *     ..
00412 *
00413 *  =====================================================================
00414 *
00415 *     .. Parameters ..
00416       REAL               ZERO, ONE
00417       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00418 *     ..
00419 *     .. Local Scalars ..
00420       LOGICAL            ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY, NOSCL,
00421      $                   PAIR, WANTSB, WANTSE, WANTSN, WANTSV
00422       CHARACTER          CHTEMP
00423       INTEGER            I, ICOLS, IERR, IJOBVL, IJOBVR, IN, IROWS,
00424      $                   ITAU, IWRK, IWRK1, J, JC, JR, M, MAXWRK,
00425      $                   MINWRK, MM
00426       REAL               ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
00427      $                   SMLNUM, TEMP
00428 *     ..
00429 *     .. Local Arrays ..
00430       LOGICAL            LDUMMA( 1 )
00431 *     ..
00432 *     .. External Subroutines ..
00433       EXTERNAL           SGEQRF, SGGBAK, SGGBAL, SGGHRD, SHGEQZ, SLABAD,
00434      $                   SLACPY, SLASCL, SLASET, SORGQR, SORMQR, STGEVC,
00435      $                   STGSNA, XERBLA
00436 *     ..
00437 *     .. External Functions ..
00438       LOGICAL            LSAME
00439       INTEGER            ILAENV
00440       REAL               SLAMCH, SLANGE
00441       EXTERNAL           LSAME, ILAENV, SLAMCH, SLANGE
00442 *     ..
00443 *     .. Intrinsic Functions ..
00444       INTRINSIC          ABS, MAX, SQRT
00445 *     ..
00446 *     .. Executable Statements ..
00447 *
00448 *     Decode the input arguments
00449 *
00450       IF( LSAME( JOBVL, 'N' ) ) THEN
00451          IJOBVL = 1
00452          ILVL = .FALSE.
00453       ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
00454          IJOBVL = 2
00455          ILVL = .TRUE.
00456       ELSE
00457          IJOBVL = -1
00458          ILVL = .FALSE.
00459       END IF
00460 *
00461       IF( LSAME( JOBVR, 'N' ) ) THEN
00462          IJOBVR = 1
00463          ILVR = .FALSE.
00464       ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
00465          IJOBVR = 2
00466          ILVR = .TRUE.
00467       ELSE
00468          IJOBVR = -1
00469          ILVR = .FALSE.
00470       END IF
00471       ILV = ILVL .OR. ILVR
00472 *
00473       NOSCL  = LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'P' )
00474       WANTSN = LSAME( SENSE, 'N' )
00475       WANTSE = LSAME( SENSE, 'E' )
00476       WANTSV = LSAME( SENSE, 'V' )
00477       WANTSB = LSAME( SENSE, 'B' )
00478 *
00479 *     Test the input arguments
00480 *
00481       INFO = 0
00482       LQUERY = ( LWORK.EQ.-1 )
00483       IF( .NOT.( NOSCL .OR. LSAME( BALANC, 'S' ) .OR.
00484      $    LSAME( BALANC, 'B' ) ) ) THEN
00485          INFO = -1
00486       ELSE IF( IJOBVL.LE.0 ) THEN
00487          INFO = -2
00488       ELSE IF( IJOBVR.LE.0 ) THEN
00489          INFO = -3
00490       ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSB .OR. WANTSV ) )
00491      $          THEN
00492          INFO = -4
00493       ELSE IF( N.LT.0 ) THEN
00494          INFO = -5
00495       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00496          INFO = -7
00497       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00498          INFO = -9
00499       ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
00500          INFO = -14
00501       ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
00502          INFO = -16
00503       END IF
00504 *
00505 *     Compute workspace
00506 *      (Note: Comments in the code beginning "Workspace:" describe the
00507 *       minimal amount of workspace needed at that point in the code,
00508 *       as well as the preferred amount for good performance.
00509 *       NB refers to the optimal block size for the immediately
00510 *       following subroutine, as returned by ILAENV. The workspace is
00511 *       computed assuming ILO = 1 and IHI = N, the worst case.)
00512 *
00513       IF( INFO.EQ.0 ) THEN
00514          IF( N.EQ.0 ) THEN
00515             MINWRK = 1
00516             MAXWRK = 1
00517          ELSE
00518             IF( NOSCL .AND. .NOT.ILV ) THEN
00519                MINWRK = 2*N
00520             ELSE
00521                MINWRK = 6*N
00522             END IF
00523             IF( WANTSE ) THEN
00524                MINWRK = 10*N
00525             ELSE IF( WANTSV .OR. WANTSB ) THEN
00526                MINWRK = 2*N*( N + 4 ) + 16
00527             END IF
00528             MAXWRK = MINWRK
00529             MAXWRK = MAX( MAXWRK,
00530      $                    N + N*ILAENV( 1, 'SGEQRF', ' ', N, 1, N, 0 ) )
00531             MAXWRK = MAX( MAXWRK,
00532      $                    N + N*ILAENV( 1, 'SORMQR', ' ', N, 1, N, 0 ) )
00533             IF( ILVL ) THEN
00534                MAXWRK = MAX( MAXWRK, N +
00535      $                       N*ILAENV( 1, 'SORGQR', ' ', N, 1, N, 0 ) )
00536             END IF
00537          END IF
00538          WORK( 1 ) = MAXWRK
00539 *
00540          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
00541             INFO = -26
00542          END IF
00543       END IF
00544 *
00545       IF( INFO.NE.0 ) THEN
00546          CALL XERBLA( 'SGGEVX', -INFO )
00547          RETURN
00548       ELSE IF( LQUERY ) THEN
00549          RETURN
00550       END IF
00551 *
00552 *     Quick return if possible
00553 *
00554       IF( N.EQ.0 )
00555      $   RETURN
00556 *
00557 *
00558 *     Get machine constants
00559 *
00560       EPS = SLAMCH( 'P' )
00561       SMLNUM = SLAMCH( 'S' )
00562       BIGNUM = ONE / SMLNUM
00563       CALL SLABAD( SMLNUM, BIGNUM )
00564       SMLNUM = SQRT( SMLNUM ) / EPS
00565       BIGNUM = ONE / SMLNUM
00566 *
00567 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00568 *
00569       ANRM = SLANGE( 'M', N, N, A, LDA, WORK )
00570       ILASCL = .FALSE.
00571       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00572          ANRMTO = SMLNUM
00573          ILASCL = .TRUE.
00574       ELSE IF( ANRM.GT.BIGNUM ) THEN
00575          ANRMTO = BIGNUM
00576          ILASCL = .TRUE.
00577       END IF
00578       IF( ILASCL )
00579      $   CALL SLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
00580 *
00581 *     Scale B if max element outside range [SMLNUM,BIGNUM]
00582 *
00583       BNRM = SLANGE( 'M', N, N, B, LDB, WORK )
00584       ILBSCL = .FALSE.
00585       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
00586          BNRMTO = SMLNUM
00587          ILBSCL = .TRUE.
00588       ELSE IF( BNRM.GT.BIGNUM ) THEN
00589          BNRMTO = BIGNUM
00590          ILBSCL = .TRUE.
00591       END IF
00592       IF( ILBSCL )
00593      $   CALL SLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
00594 *
00595 *     Permute and/or balance the matrix pair (A,B)
00596 *     (Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
00597 *
00598       CALL SGGBAL( BALANC, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE,
00599      $             WORK, IERR )
00600 *
00601 *     Compute ABNRM and BBNRM
00602 *
00603       ABNRM = SLANGE( '1', N, N, A, LDA, WORK( 1 ) )
00604       IF( ILASCL ) THEN
00605          WORK( 1 ) = ABNRM
00606          CALL SLASCL( 'G', 0, 0, ANRMTO, ANRM, 1, 1, WORK( 1 ), 1,
00607      $                IERR )
00608          ABNRM = WORK( 1 )
00609       END IF
00610 *
00611       BBNRM = SLANGE( '1', N, N, B, LDB, WORK( 1 ) )
00612       IF( ILBSCL ) THEN
00613          WORK( 1 ) = BBNRM
00614          CALL SLASCL( 'G', 0, 0, BNRMTO, BNRM, 1, 1, WORK( 1 ), 1,
00615      $                IERR )
00616          BBNRM = WORK( 1 )
00617       END IF
00618 *
00619 *     Reduce B to triangular form (QR decomposition of B)
00620 *     (Workspace: need N, prefer N*NB )
00621 *
00622       IROWS = IHI + 1 - ILO
00623       IF( ILV .OR. .NOT.WANTSN ) THEN
00624          ICOLS = N + 1 - ILO
00625       ELSE
00626          ICOLS = IROWS
00627       END IF
00628       ITAU = 1
00629       IWRK = ITAU + IROWS
00630       CALL SGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
00631      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
00632 *
00633 *     Apply the orthogonal transformation to A
00634 *     (Workspace: need N, prefer N*NB)
00635 *
00636       CALL SORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
00637      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
00638      $             LWORK+1-IWRK, IERR )
00639 *
00640 *     Initialize VL and/or VR
00641 *     (Workspace: need N, prefer N*NB)
00642 *
00643       IF( ILVL ) THEN
00644          CALL SLASET( 'Full', N, N, ZERO, ONE, VL, LDVL )
00645          IF( IROWS.GT.1 ) THEN
00646             CALL SLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
00647      $                   VL( ILO+1, ILO ), LDVL )
00648          END IF
00649          CALL SORGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
00650      $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
00651       END IF
00652 *
00653       IF( ILVR )
00654      $   CALL SLASET( 'Full', N, N, ZERO, ONE, VR, LDVR )
00655 *
00656 *     Reduce to generalized Hessenberg form
00657 *     (Workspace: none needed)
00658 *
00659       IF( ILV .OR. .NOT.WANTSN ) THEN
00660 *
00661 *        Eigenvectors requested -- work on whole matrix.
00662 *
00663          CALL SGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
00664      $                LDVL, VR, LDVR, IERR )
00665       ELSE
00666          CALL SGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
00667      $                B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
00668       END IF
00669 *
00670 *     Perform QZ algorithm (Compute eigenvalues, and optionally, the
00671 *     Schur forms and Schur vectors)
00672 *     (Workspace: need N)
00673 *
00674       IF( ILV .OR. .NOT.WANTSN ) THEN
00675          CHTEMP = 'S'
00676       ELSE
00677          CHTEMP = 'E'
00678       END IF
00679 *
00680       CALL SHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
00681      $             ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, WORK,
00682      $             LWORK, IERR )
00683       IF( IERR.NE.0 ) THEN
00684          IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
00685             INFO = IERR
00686          ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
00687             INFO = IERR - N
00688          ELSE
00689             INFO = N + 1
00690          END IF
00691          GO TO 130
00692       END IF
00693 *
00694 *     Compute Eigenvectors and estimate condition numbers if desired
00695 *     (Workspace: STGEVC: need 6*N
00696 *                 STGSNA: need 2*N*(N+2)+16 if SENSE = 'V' or 'B',
00697 *                         need N otherwise )
00698 *
00699       IF( ILV .OR. .NOT.WANTSN ) THEN
00700          IF( ILV ) THEN
00701             IF( ILVL ) THEN
00702                IF( ILVR ) THEN
00703                   CHTEMP = 'B'
00704                ELSE
00705                   CHTEMP = 'L'
00706                END IF
00707             ELSE
00708                CHTEMP = 'R'
00709             END IF
00710 *
00711             CALL STGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL,
00712      $                   LDVL, VR, LDVR, N, IN, WORK, IERR )
00713             IF( IERR.NE.0 ) THEN
00714                INFO = N + 2
00715                GO TO 130
00716             END IF
00717          END IF
00718 *
00719          IF( .NOT.WANTSN ) THEN
00720 *
00721 *           compute eigenvectors (STGEVC) and estimate condition
00722 *           numbers (STGSNA). Note that the definition of the condition
00723 *           number is not invariant under transformation (u,v) to
00724 *           (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
00725 *           Schur form (S,T), Q and Z are orthogonal matrices. In order
00726 *           to avoid using extra 2*N*N workspace, we have to recalculate
00727 *           eigenvectors and estimate one condition numbers at a time.
00728 *
00729             PAIR = .FALSE.
00730             DO 20 I = 1, N
00731 *
00732                IF( PAIR ) THEN
00733                   PAIR = .FALSE.
00734                   GO TO 20
00735                END IF
00736                MM = 1
00737                IF( I.LT.N ) THEN
00738                   IF( A( I+1, I ).NE.ZERO ) THEN
00739                      PAIR = .TRUE.
00740                      MM = 2
00741                   END IF
00742                END IF
00743 *
00744                DO 10 J = 1, N
00745                   BWORK( J ) = .FALSE.
00746    10          CONTINUE
00747                IF( MM.EQ.1 ) THEN
00748                   BWORK( I ) = .TRUE.
00749                ELSE IF( MM.EQ.2 ) THEN
00750                   BWORK( I ) = .TRUE.
00751                   BWORK( I+1 ) = .TRUE.
00752                END IF
00753 *
00754                IWRK = MM*N + 1
00755                IWRK1 = IWRK + MM*N
00756 *
00757 *              Compute a pair of left and right eigenvectors.
00758 *              (compute workspace: need up to 4*N + 6*N)
00759 *
00760                IF( WANTSE .OR. WANTSB ) THEN
00761                   CALL STGEVC( 'B', 'S', BWORK, N, A, LDA, B, LDB,
00762      $                         WORK( 1 ), N, WORK( IWRK ), N, MM, M,
00763      $                         WORK( IWRK1 ), IERR )
00764                   IF( IERR.NE.0 ) THEN
00765                      INFO = N + 2
00766                      GO TO 130
00767                   END IF
00768                END IF
00769 *
00770                CALL STGSNA( SENSE, 'S', BWORK, N, A, LDA, B, LDB,
00771      $                      WORK( 1 ), N, WORK( IWRK ), N, RCONDE( I ),
00772      $                      RCONDV( I ), MM, M, WORK( IWRK1 ),
00773      $                      LWORK-IWRK1+1, IWORK, IERR )
00774 *
00775    20       CONTINUE
00776          END IF
00777       END IF
00778 *
00779 *     Undo balancing on VL and VR and normalization
00780 *     (Workspace: none needed)
00781 *
00782       IF( ILVL ) THEN
00783          CALL SGGBAK( BALANC, 'L', N, ILO, IHI, LSCALE, RSCALE, N, VL,
00784      $                LDVL, IERR )
00785 *
00786          DO 70 JC = 1, N
00787             IF( ALPHAI( JC ).LT.ZERO )
00788      $         GO TO 70
00789             TEMP = ZERO
00790             IF( ALPHAI( JC ).EQ.ZERO ) THEN
00791                DO 30 JR = 1, N
00792                   TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
00793    30          CONTINUE
00794             ELSE
00795                DO 40 JR = 1, N
00796                   TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
00797      $                   ABS( VL( JR, JC+1 ) ) )
00798    40          CONTINUE
00799             END IF
00800             IF( TEMP.LT.SMLNUM )
00801      $         GO TO 70
00802             TEMP = ONE / TEMP
00803             IF( ALPHAI( JC ).EQ.ZERO ) THEN
00804                DO 50 JR = 1, N
00805                   VL( JR, JC ) = VL( JR, JC )*TEMP
00806    50          CONTINUE
00807             ELSE
00808                DO 60 JR = 1, N
00809                   VL( JR, JC ) = VL( JR, JC )*TEMP
00810                   VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
00811    60          CONTINUE
00812             END IF
00813    70    CONTINUE
00814       END IF
00815       IF( ILVR ) THEN
00816          CALL SGGBAK( BALANC, 'R', N, ILO, IHI, LSCALE, RSCALE, N, VR,
00817      $                LDVR, IERR )
00818          DO 120 JC = 1, N
00819             IF( ALPHAI( JC ).LT.ZERO )
00820      $         GO TO 120
00821             TEMP = ZERO
00822             IF( ALPHAI( JC ).EQ.ZERO ) THEN
00823                DO 80 JR = 1, N
00824                   TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
00825    80          CONTINUE
00826             ELSE
00827                DO 90 JR = 1, N
00828                   TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
00829      $                   ABS( VR( JR, JC+1 ) ) )
00830    90          CONTINUE
00831             END IF
00832             IF( TEMP.LT.SMLNUM )
00833      $         GO TO 120
00834             TEMP = ONE / TEMP
00835             IF( ALPHAI( JC ).EQ.ZERO ) THEN
00836                DO 100 JR = 1, N
00837                   VR( JR, JC ) = VR( JR, JC )*TEMP
00838   100          CONTINUE
00839             ELSE
00840                DO 110 JR = 1, N
00841                   VR( JR, JC ) = VR( JR, JC )*TEMP
00842                   VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
00843   110          CONTINUE
00844             END IF
00845   120    CONTINUE
00846       END IF
00847 *
00848 *     Undo scaling if necessary
00849 *
00850   130 CONTINUE
00851 *
00852       IF( ILASCL ) THEN
00853          CALL SLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
00854          CALL SLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
00855       END IF
00856 *
00857       IF( ILBSCL ) THEN
00858          CALL SLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
00859       END IF
00860 *
00861       WORK( 1 ) = MAXWRK
00862       RETURN
00863 *
00864 *     End of SGGEVX
00865 *
00866       END
 All Files Functions