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