LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dggevx.f
Go to the documentation of this file.
00001 *> \brief <b> DGGEVX 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 DGGEVX + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggevx.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggevx.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggevx.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DGGEVX( 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 *       DOUBLE PRECISION   ABNRM, BBNRM
00030 *       ..
00031 *       .. Array Arguments ..
00032 *       LOGICAL            BWORK( * )
00033 *       INTEGER            IWORK( * )
00034 *       DOUBLE PRECISION   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 *> DGGEVX 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
00156 *> \endverbatim
00157 *>
00158 *> \param[out] ALPHAI
00159 *> \verbatim
00160 *>          ALPHAI is DOUBLE PRECISION array, dimension (N)
00161 *> \endverbatim
00162 *>
00163 *> \param[out] BETA
00164 *> \verbatim
00165 *>          BETA is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
00267 *>          The one-norm of the balanced matrix A.
00268 *> \endverbatim
00269 *>
00270 *> \param[out] BBNRM
00271 *> \verbatim
00272 *>          BBNRM is DOUBLE PRECISION
00273 *>          The one-norm of the balanced matrix B.
00274 *> \endverbatim
00275 *>
00276 *> \param[out] RCONDE
00277 *> \verbatim
00278 *>          RCONDE is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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' or 'B', 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 DHGEQZ.
00344 *>                =N+2: error return from DTGEVC.
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 doubleGEeigen
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 DGGEVX( 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       DOUBLE PRECISION   ABNRM, BBNRM
00403 *     ..
00404 *     .. Array Arguments ..
00405       LOGICAL            BWORK( * )
00406       INTEGER            IWORK( * )
00407       DOUBLE PRECISION   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       DOUBLE PRECISION   ZERO, ONE
00417       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+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       DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
00427      $                   SMLNUM, TEMP
00428 *     ..
00429 *     .. Local Arrays ..
00430       LOGICAL            LDUMMA( 1 )
00431 *     ..
00432 *     .. External Subroutines ..
00433       EXTERNAL           DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLABAD,
00434      $                   DLACPY, DLASCL, DLASET, DORGQR, DORMQR, DTGEVC,
00435      $                   DTGSNA, XERBLA 
00436 *     ..
00437 *     .. External Functions ..
00438       LOGICAL            LSAME
00439       INTEGER            ILAENV
00440       DOUBLE PRECISION   DLAMCH, DLANGE
00441       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANGE
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.( LSAME( BALANC, 'N' ) .OR. LSAME( BALANC,
00484      $    'S' ) .OR. LSAME( BALANC, 'P' ) .OR. LSAME( BALANC, 'B' ) ) )
00485      $     THEN
00486          INFO = -1
00487       ELSE IF( IJOBVL.LE.0 ) THEN
00488          INFO = -2
00489       ELSE IF( IJOBVR.LE.0 ) THEN
00490          INFO = -3
00491       ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSB .OR. WANTSV ) )
00492      $          THEN
00493          INFO = -4
00494       ELSE IF( N.LT.0 ) THEN
00495          INFO = -5
00496       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00497          INFO = -7
00498       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00499          INFO = -9
00500       ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
00501          INFO = -14
00502       ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
00503          INFO = -16
00504       END IF
00505 *
00506 *     Compute workspace
00507 *      (Note: Comments in the code beginning "Workspace:" describe the
00508 *       minimal amount of workspace needed at that point in the code,
00509 *       as well as the preferred amount for good performance.
00510 *       NB refers to the optimal block size for the immediately
00511 *       following subroutine, as returned by ILAENV. The workspace is
00512 *       computed assuming ILO = 1 and IHI = N, the worst case.)
00513 *
00514       IF( INFO.EQ.0 ) THEN
00515          IF( N.EQ.0 ) THEN
00516             MINWRK = 1
00517             MAXWRK = 1
00518          ELSE
00519             IF( NOSCL .AND. .NOT.ILV ) THEN
00520                MINWRK = 2*N
00521             ELSE
00522                MINWRK = 6*N
00523             END IF
00524             IF( WANTSE .OR. WANTSB ) THEN
00525                MINWRK = 10*N
00526             END IF
00527             IF( WANTSV .OR. WANTSB ) THEN
00528                MINWRK = MAX( MINWRK, 2*N*( N + 4 ) + 16 )
00529             END IF
00530             MAXWRK = MINWRK
00531             MAXWRK = MAX( MAXWRK,
00532      $                    N + N*ILAENV( 1, 'DGEQRF', ' ', N, 1, N, 0 ) )
00533             MAXWRK = MAX( MAXWRK,
00534      $                    N + N*ILAENV( 1, 'DORMQR', ' ', N, 1, N, 0 ) )
00535             IF( ILVL ) THEN
00536                MAXWRK = MAX( MAXWRK, N +
00537      $                       N*ILAENV( 1, 'DORGQR', ' ', N, 1, N, 0 ) )
00538             END IF
00539          END IF
00540          WORK( 1 ) = MAXWRK
00541 *
00542          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
00543             INFO = -26
00544          END IF
00545       END IF
00546 *
00547       IF( INFO.NE.0 ) THEN
00548          CALL XERBLA( 'DGGEVX', -INFO )
00549          RETURN
00550       ELSE IF( LQUERY ) THEN
00551          RETURN
00552       END IF
00553 *
00554 *     Quick return if possible
00555 *
00556       IF( N.EQ.0 )
00557      $   RETURN
00558 *
00559 *
00560 *     Get machine constants
00561 *
00562       EPS = DLAMCH( 'P' )
00563       SMLNUM = DLAMCH( 'S' )
00564       BIGNUM = ONE / SMLNUM
00565       CALL DLABAD( SMLNUM, BIGNUM )
00566       SMLNUM = SQRT( SMLNUM ) / EPS
00567       BIGNUM = ONE / SMLNUM
00568 *
00569 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00570 *
00571       ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
00572       ILASCL = .FALSE.
00573       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00574          ANRMTO = SMLNUM
00575          ILASCL = .TRUE.
00576       ELSE IF( ANRM.GT.BIGNUM ) THEN
00577          ANRMTO = BIGNUM
00578          ILASCL = .TRUE.
00579       END IF
00580       IF( ILASCL )
00581      $   CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
00582 *
00583 *     Scale B if max element outside range [SMLNUM,BIGNUM]
00584 *
00585       BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
00586       ILBSCL = .FALSE.
00587       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
00588          BNRMTO = SMLNUM
00589          ILBSCL = .TRUE.
00590       ELSE IF( BNRM.GT.BIGNUM ) THEN
00591          BNRMTO = BIGNUM
00592          ILBSCL = .TRUE.
00593       END IF
00594       IF( ILBSCL )
00595      $   CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
00596 *
00597 *     Permute and/or balance the matrix pair (A,B)
00598 *     (Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
00599 *
00600       CALL DGGBAL( BALANC, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE,
00601      $             WORK, IERR )
00602 *
00603 *     Compute ABNRM and BBNRM
00604 *
00605       ABNRM = DLANGE( '1', N, N, A, LDA, WORK( 1 ) )
00606       IF( ILASCL ) THEN
00607          WORK( 1 ) = ABNRM
00608          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, 1, 1, WORK( 1 ), 1,
00609      $                IERR )
00610          ABNRM = WORK( 1 )
00611       END IF
00612 *
00613       BBNRM = DLANGE( '1', N, N, B, LDB, WORK( 1 ) )
00614       IF( ILBSCL ) THEN
00615          WORK( 1 ) = BBNRM
00616          CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, 1, 1, WORK( 1 ), 1,
00617      $                IERR )
00618          BBNRM = WORK( 1 )
00619       END IF
00620 *
00621 *     Reduce B to triangular form (QR decomposition of B)
00622 *     (Workspace: need N, prefer N*NB )
00623 *
00624       IROWS = IHI + 1 - ILO
00625       IF( ILV .OR. .NOT.WANTSN ) THEN
00626          ICOLS = N + 1 - ILO
00627       ELSE
00628          ICOLS = IROWS
00629       END IF
00630       ITAU = 1
00631       IWRK = ITAU + IROWS
00632       CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
00633      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
00634 *
00635 *     Apply the orthogonal transformation to A
00636 *     (Workspace: need N, prefer N*NB)
00637 *
00638       CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
00639      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
00640      $             LWORK+1-IWRK, IERR )
00641 *
00642 *     Initialize VL and/or VR
00643 *     (Workspace: need N, prefer N*NB)
00644 *
00645       IF( ILVL ) THEN
00646          CALL DLASET( 'Full', N, N, ZERO, ONE, VL, LDVL )
00647          IF( IROWS.GT.1 ) THEN
00648             CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
00649      $                   VL( ILO+1, ILO ), LDVL )
00650          END IF
00651          CALL DORGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
00652      $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
00653       END IF
00654 *
00655       IF( ILVR )
00656      $   CALL DLASET( 'Full', N, N, ZERO, ONE, VR, LDVR )
00657 *
00658 *     Reduce to generalized Hessenberg form
00659 *     (Workspace: none needed)
00660 *
00661       IF( ILV .OR. .NOT.WANTSN ) THEN
00662 *
00663 *        Eigenvectors requested -- work on whole matrix.
00664 *
00665          CALL DGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
00666      $                LDVL, VR, LDVR, IERR )
00667       ELSE
00668          CALL DGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
00669      $                B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
00670       END IF
00671 *
00672 *     Perform QZ algorithm (Compute eigenvalues, and optionally, the
00673 *     Schur forms and Schur vectors)
00674 *     (Workspace: need N)
00675 *
00676       IF( ILV .OR. .NOT.WANTSN ) THEN
00677          CHTEMP = 'S'
00678       ELSE
00679          CHTEMP = 'E'
00680       END IF
00681 *
00682       CALL DHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
00683      $             ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, WORK,
00684      $             LWORK, IERR )
00685       IF( IERR.NE.0 ) THEN
00686          IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
00687             INFO = IERR
00688          ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
00689             INFO = IERR - N
00690          ELSE
00691             INFO = N + 1
00692          END IF
00693          GO TO 130
00694       END IF
00695 *
00696 *     Compute Eigenvectors and estimate condition numbers if desired
00697 *     (Workspace: DTGEVC: need 6*N
00698 *                 DTGSNA: need 2*N*(N+2)+16 if SENSE = 'V' or 'B',
00699 *                         need N otherwise )
00700 *
00701       IF( ILV .OR. .NOT.WANTSN ) THEN
00702          IF( ILV ) THEN
00703             IF( ILVL ) THEN
00704                IF( ILVR ) THEN
00705                   CHTEMP = 'B'
00706                ELSE
00707                   CHTEMP = 'L'
00708                END IF
00709             ELSE
00710                CHTEMP = 'R'
00711             END IF
00712 *
00713             CALL DTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL,
00714      $                   LDVL, VR, LDVR, N, IN, WORK, IERR )
00715             IF( IERR.NE.0 ) THEN
00716                INFO = N + 2
00717                GO TO 130
00718             END IF
00719          END IF
00720 *
00721          IF( .NOT.WANTSN ) THEN
00722 *
00723 *           compute eigenvectors (DTGEVC) and estimate condition
00724 *           numbers (DTGSNA). Note that the definition of the condition
00725 *           number is not invariant under transformation (u,v) to
00726 *           (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
00727 *           Schur form (S,T), Q and Z are orthogonal matrices. In order
00728 *           to avoid using extra 2*N*N workspace, we have to recalculate
00729 *           eigenvectors and estimate one condition numbers at a time.
00730 *
00731             PAIR = .FALSE.
00732             DO 20 I = 1, N
00733 *
00734                IF( PAIR ) THEN
00735                   PAIR = .FALSE.
00736                   GO TO 20
00737                END IF
00738                MM = 1
00739                IF( I.LT.N ) THEN
00740                   IF( A( I+1, I ).NE.ZERO ) THEN
00741                      PAIR = .TRUE.
00742                      MM = 2
00743                   END IF
00744                END IF
00745 *
00746                DO 10 J = 1, N
00747                   BWORK( J ) = .FALSE.
00748    10          CONTINUE
00749                IF( MM.EQ.1 ) THEN
00750                   BWORK( I ) = .TRUE.
00751                ELSE IF( MM.EQ.2 ) THEN
00752                   BWORK( I ) = .TRUE.
00753                   BWORK( I+1 ) = .TRUE.
00754                END IF
00755 *
00756                IWRK = MM*N + 1
00757                IWRK1 = IWRK + MM*N
00758 *
00759 *              Compute a pair of left and right eigenvectors.
00760 *              (compute workspace: need up to 4*N + 6*N)
00761 *
00762                IF( WANTSE .OR. WANTSB ) THEN
00763                   CALL DTGEVC( 'B', 'S', BWORK, N, A, LDA, B, LDB,
00764      $                         WORK( 1 ), N, WORK( IWRK ), N, MM, M,
00765      $                         WORK( IWRK1 ), IERR )
00766                   IF( IERR.NE.0 ) THEN
00767                      INFO = N + 2
00768                      GO TO 130
00769                   END IF
00770                END IF
00771 *
00772                CALL DTGSNA( SENSE, 'S', BWORK, N, A, LDA, B, LDB,
00773      $                      WORK( 1 ), N, WORK( IWRK ), N, RCONDE( I ),
00774      $                      RCONDV( I ), MM, M, WORK( IWRK1 ),
00775      $                      LWORK-IWRK1+1, IWORK, IERR )
00776 *
00777    20       CONTINUE
00778          END IF
00779       END IF
00780 *
00781 *     Undo balancing on VL and VR and normalization
00782 *     (Workspace: none needed)
00783 *
00784       IF( ILVL ) THEN
00785          CALL DGGBAK( BALANC, 'L', N, ILO, IHI, LSCALE, RSCALE, N, VL,
00786      $                LDVL, IERR )
00787 *
00788          DO 70 JC = 1, N
00789             IF( ALPHAI( JC ).LT.ZERO )
00790      $         GO TO 70
00791             TEMP = ZERO
00792             IF( ALPHAI( JC ).EQ.ZERO ) THEN
00793                DO 30 JR = 1, N
00794                   TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
00795    30          CONTINUE
00796             ELSE
00797                DO 40 JR = 1, N
00798                   TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
00799      $                   ABS( VL( JR, JC+1 ) ) )
00800    40          CONTINUE
00801             END IF
00802             IF( TEMP.LT.SMLNUM )
00803      $         GO TO 70
00804             TEMP = ONE / TEMP
00805             IF( ALPHAI( JC ).EQ.ZERO ) THEN
00806                DO 50 JR = 1, N
00807                   VL( JR, JC ) = VL( JR, JC )*TEMP
00808    50          CONTINUE
00809             ELSE
00810                DO 60 JR = 1, N
00811                   VL( JR, JC ) = VL( JR, JC )*TEMP
00812                   VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
00813    60          CONTINUE
00814             END IF
00815    70    CONTINUE
00816       END IF
00817       IF( ILVR ) THEN
00818          CALL DGGBAK( BALANC, 'R', N, ILO, IHI, LSCALE, RSCALE, N, VR,
00819      $                LDVR, IERR )
00820          DO 120 JC = 1, N
00821             IF( ALPHAI( JC ).LT.ZERO )
00822      $         GO TO 120
00823             TEMP = ZERO
00824             IF( ALPHAI( JC ).EQ.ZERO ) THEN
00825                DO 80 JR = 1, N
00826                   TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
00827    80          CONTINUE
00828             ELSE
00829                DO 90 JR = 1, N
00830                   TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
00831      $                   ABS( VR( JR, JC+1 ) ) )
00832    90          CONTINUE
00833             END IF
00834             IF( TEMP.LT.SMLNUM )
00835      $         GO TO 120
00836             TEMP = ONE / TEMP
00837             IF( ALPHAI( JC ).EQ.ZERO ) THEN
00838                DO 100 JR = 1, N
00839                   VR( JR, JC ) = VR( JR, JC )*TEMP
00840   100          CONTINUE
00841             ELSE
00842                DO 110 JR = 1, N
00843                   VR( JR, JC ) = VR( JR, JC )*TEMP
00844                   VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
00845   110          CONTINUE
00846             END IF
00847   120    CONTINUE
00848       END IF
00849 *
00850 *     Undo scaling if necessary
00851 *
00852   130 CONTINUE
00853 *
00854       IF( ILASCL ) THEN
00855          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
00856          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
00857       END IF
00858 *
00859       IF( ILBSCL ) THEN
00860          CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
00861       END IF
00862 *
00863       WORK( 1 ) = MAXWRK
00864       RETURN
00865 *
00866 *     End of DGGEVX
00867 *
00868       END
 All Files Functions