LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dgegv.f
Go to the documentation of this file.
00001 *> \brief <b> DGEEVX 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 DGEGV + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgegv.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgegv.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgegv.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
00022 *                         BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          JOBVL, JOBVR
00026 *       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
00030 *      $                   B( LDB, * ), BETA( * ), VL( LDVL, * ),
00031 *      $                   VR( LDVR, * ), WORK( * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> This routine is deprecated and has been replaced by routine DGGEV.
00041 *>
00042 *> DGEGV computes the eigenvalues and, optionally, the left and/or right
00043 *> eigenvectors of a real matrix pair (A,B).
00044 *> Given two square matrices A and B,
00045 *> the generalized nonsymmetric eigenvalue problem (GNEP) is to find the
00046 *> eigenvalues lambda and corresponding (non-zero) eigenvectors x such
00047 *> that
00048 *>
00049 *>    A*x = lambda*B*x.
00050 *>
00051 *> An alternate form is to find the eigenvalues mu and corresponding
00052 *> eigenvectors y such that
00053 *>
00054 *>    mu*A*y = B*y.
00055 *>
00056 *> These two forms are equivalent with mu = 1/lambda and x = y if
00057 *> neither lambda nor mu is zero.  In order to deal with the case that
00058 *> lambda or mu is zero or small, two values alpha and beta are returned
00059 *> for each eigenvalue, such that lambda = alpha/beta and
00060 *> mu = beta/alpha.
00061 *>
00062 *> The vectors x and y in the above equations are right eigenvectors of
00063 *> the matrix pair (A,B).  Vectors u and v satisfying
00064 *>
00065 *>    u**H*A = lambda*u**H*B  or  mu*v**H*A = v**H*B
00066 *>
00067 *> are left eigenvectors of (A,B).
00068 *>
00069 *> Note: this routine performs "full balancing" on A and B
00070 *> \endverbatim
00071 *
00072 *  Arguments:
00073 *  ==========
00074 *
00075 *> \param[in] JOBVL
00076 *> \verbatim
00077 *>          JOBVL is CHARACTER*1
00078 *>          = 'N':  do not compute the left generalized eigenvectors;
00079 *>          = 'V':  compute the left generalized eigenvectors (returned
00080 *>                  in VL).
00081 *> \endverbatim
00082 *>
00083 *> \param[in] JOBVR
00084 *> \verbatim
00085 *>          JOBVR is CHARACTER*1
00086 *>          = 'N':  do not compute the right generalized eigenvectors;
00087 *>          = 'V':  compute the right generalized eigenvectors (returned
00088 *>                  in VR).
00089 *> \endverbatim
00090 *>
00091 *> \param[in] N
00092 *> \verbatim
00093 *>          N is INTEGER
00094 *>          The order of the matrices A, B, VL, and VR.  N >= 0.
00095 *> \endverbatim
00096 *>
00097 *> \param[in,out] A
00098 *> \verbatim
00099 *>          A is DOUBLE PRECISION array, dimension (LDA, N)
00100 *>          On entry, the matrix A.
00101 *>          If JOBVL = 'V' or JOBVR = 'V', then on exit A
00102 *>          contains the real Schur form of A from the generalized Schur
00103 *>          factorization of the pair (A,B) after balancing.
00104 *>          If no eigenvectors were computed, then only the diagonal
00105 *>          blocks from the Schur form will be correct.  See DGGHRD and
00106 *>          DHGEQZ for details.
00107 *> \endverbatim
00108 *>
00109 *> \param[in] LDA
00110 *> \verbatim
00111 *>          LDA is INTEGER
00112 *>          The leading dimension of A.  LDA >= max(1,N).
00113 *> \endverbatim
00114 *>
00115 *> \param[in,out] B
00116 *> \verbatim
00117 *>          B is DOUBLE PRECISION array, dimension (LDB, N)
00118 *>          On entry, the matrix B.
00119 *>          If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the
00120 *>          upper triangular matrix obtained from B in the generalized
00121 *>          Schur factorization of the pair (A,B) after balancing.
00122 *>          If no eigenvectors were computed, then only those elements of
00123 *>          B corresponding to the diagonal blocks from the Schur form of
00124 *>          A will be correct.  See DGGHRD and DHGEQZ for details.
00125 *> \endverbatim
00126 *>
00127 *> \param[in] LDB
00128 *> \verbatim
00129 *>          LDB is INTEGER
00130 *>          The leading dimension of B.  LDB >= max(1,N).
00131 *> \endverbatim
00132 *>
00133 *> \param[out] ALPHAR
00134 *> \verbatim
00135 *>          ALPHAR is DOUBLE PRECISION array, dimension (N)
00136 *>          The real parts of each scalar alpha defining an eigenvalue of
00137 *>          GNEP.
00138 *> \endverbatim
00139 *>
00140 *> \param[out] ALPHAI
00141 *> \verbatim
00142 *>          ALPHAI is DOUBLE PRECISION array, dimension (N)
00143 *>          The imaginary parts of each scalar alpha defining an
00144 *>          eigenvalue of GNEP.  If ALPHAI(j) is zero, then the j-th
00145 *>          eigenvalue is real; if positive, then the j-th and
00146 *>          (j+1)-st eigenvalues are a complex conjugate pair, with
00147 *>          ALPHAI(j+1) = -ALPHAI(j).
00148 *> \endverbatim
00149 *>
00150 *> \param[out] BETA
00151 *> \verbatim
00152 *>          BETA is DOUBLE PRECISION array, dimension (N)
00153 *>          The scalars beta that define the eigenvalues of GNEP.
00154 *>          
00155 *>          Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
00156 *>          beta = BETA(j) represent the j-th eigenvalue of the matrix
00157 *>          pair (A,B), in one of the forms lambda = alpha/beta or
00158 *>          mu = beta/alpha.  Since either lambda or mu may overflow,
00159 *>          they should not, in general, be computed.
00160 *> \endverbatim
00161 *>
00162 *> \param[out] VL
00163 *> \verbatim
00164 *>          VL is DOUBLE PRECISION array, dimension (LDVL,N)
00165 *>          If JOBVL = 'V', the left eigenvectors u(j) are stored
00166 *>          in the columns of VL, in the same order as their eigenvalues.
00167 *>          If the j-th eigenvalue is real, then u(j) = VL(:,j).
00168 *>          If the j-th and (j+1)-st eigenvalues form a complex conjugate
00169 *>          pair, then
00170 *>             u(j) = VL(:,j) + i*VL(:,j+1)
00171 *>          and
00172 *>            u(j+1) = VL(:,j) - i*VL(:,j+1).
00173 *>
00174 *>          Each eigenvector is scaled so that its largest component has
00175 *>          abs(real part) + abs(imag. part) = 1, except for eigenvectors
00176 *>          corresponding to an eigenvalue with alpha = beta = 0, which
00177 *>          are set to zero.
00178 *>          Not referenced if JOBVL = 'N'.
00179 *> \endverbatim
00180 *>
00181 *> \param[in] LDVL
00182 *> \verbatim
00183 *>          LDVL is INTEGER
00184 *>          The leading dimension of the matrix VL. LDVL >= 1, and
00185 *>          if JOBVL = 'V', LDVL >= N.
00186 *> \endverbatim
00187 *>
00188 *> \param[out] VR
00189 *> \verbatim
00190 *>          VR is DOUBLE PRECISION array, dimension (LDVR,N)
00191 *>          If JOBVR = 'V', the right eigenvectors x(j) are stored
00192 *>          in the columns of VR, in the same order as their eigenvalues.
00193 *>          If the j-th eigenvalue is real, then x(j) = VR(:,j).
00194 *>          If the j-th and (j+1)-st eigenvalues form a complex conjugate
00195 *>          pair, then
00196 *>            x(j) = VR(:,j) + i*VR(:,j+1)
00197 *>          and
00198 *>            x(j+1) = VR(:,j) - i*VR(:,j+1).
00199 *>
00200 *>          Each eigenvector is scaled so that its largest component has
00201 *>          abs(real part) + abs(imag. part) = 1, except for eigenvalues
00202 *>          corresponding to an eigenvalue with alpha = beta = 0, which
00203 *>          are set to zero.
00204 *>          Not referenced if JOBVR = 'N'.
00205 *> \endverbatim
00206 *>
00207 *> \param[in] LDVR
00208 *> \verbatim
00209 *>          LDVR is INTEGER
00210 *>          The leading dimension of the matrix VR. LDVR >= 1, and
00211 *>          if JOBVR = 'V', LDVR >= N.
00212 *> \endverbatim
00213 *>
00214 *> \param[out] WORK
00215 *> \verbatim
00216 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
00217 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00218 *> \endverbatim
00219 *>
00220 *> \param[in] LWORK
00221 *> \verbatim
00222 *>          LWORK is INTEGER
00223 *>          The dimension of the array WORK.  LWORK >= max(1,8*N).
00224 *>          For good performance, LWORK must generally be larger.
00225 *>          To compute the optimal value of LWORK, call ILAENV to get
00226 *>          blocksizes (for DGEQRF, DORMQR, and DORGQR.)  Then compute:
00227 *>          NB  -- MAX of the blocksizes for DGEQRF, DORMQR, and DORGQR;
00228 *>          The optimal LWORK is:
00229 *>              2*N + MAX( 6*N, N*(NB+1) ).
00230 *>
00231 *>          If LWORK = -1, then a workspace query is assumed; the routine
00232 *>          only calculates the optimal size of the WORK array, returns
00233 *>          this value as the first entry of the WORK array, and no error
00234 *>          message related to LWORK is issued by XERBLA.
00235 *> \endverbatim
00236 *>
00237 *> \param[out] INFO
00238 *> \verbatim
00239 *>          INFO is INTEGER
00240 *>          = 0:  successful exit
00241 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00242 *>          = 1,...,N:
00243 *>                The QZ iteration failed.  No eigenvectors have been
00244 *>                calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
00245 *>                should be correct for j=INFO+1,...,N.
00246 *>          > N:  errors that usually indicate LAPACK problems:
00247 *>                =N+1: error return from DGGBAL
00248 *>                =N+2: error return from DGEQRF
00249 *>                =N+3: error return from DORMQR
00250 *>                =N+4: error return from DORGQR
00251 *>                =N+5: error return from DGGHRD
00252 *>                =N+6: error return from DHGEQZ (other than failed
00253 *>                                                iteration)
00254 *>                =N+7: error return from DTGEVC
00255 *>                =N+8: error return from DGGBAK (computing VL)
00256 *>                =N+9: error return from DGGBAK (computing VR)
00257 *>                =N+10: error return from DLASCL (various calls)
00258 *> \endverbatim
00259 *
00260 *  Authors:
00261 *  ========
00262 *
00263 *> \author Univ. of Tennessee 
00264 *> \author Univ. of California Berkeley 
00265 *> \author Univ. of Colorado Denver 
00266 *> \author NAG Ltd. 
00267 *
00268 *> \date November 2011
00269 *
00270 *> \ingroup doubleGEeigen
00271 *
00272 *> \par Further Details:
00273 *  =====================
00274 *>
00275 *> \verbatim
00276 *>
00277 *>  Balancing
00278 *>  ---------
00279 *>
00280 *>  This driver calls DGGBAL to both permute and scale rows and columns
00281 *>  of A and B.  The permutations PL and PR are chosen so that PL*A*PR
00282 *>  and PL*B*R will be upper triangular except for the diagonal blocks
00283 *>  A(i:j,i:j) and B(i:j,i:j), with i and j as close together as
00284 *>  possible.  The diagonal scaling matrices DL and DR are chosen so
00285 *>  that the pair  DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to
00286 *>  one (except for the elements that start out zero.)
00287 *>
00288 *>  After the eigenvalues and eigenvectors of the balanced matrices
00289 *>  have been computed, DGGBAK transforms the eigenvectors back to what
00290 *>  they would have been (in perfect arithmetic) if they had not been
00291 *>  balanced.
00292 *>
00293 *>  Contents of A and B on Exit
00294 *>  -------- -- - --- - -- ----
00295 *>
00296 *>  If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or
00297 *>  both), then on exit the arrays A and B will contain the real Schur
00298 *>  form[*] of the "balanced" versions of A and B.  If no eigenvectors
00299 *>  are computed, then only the diagonal blocks will be correct.
00300 *>
00301 *>  [*] See DHGEQZ, DGEGS, or read the book "Matrix Computations",
00302 *>      by Golub & van Loan, pub. by Johns Hopkins U. Press.
00303 *> \endverbatim
00304 *>
00305 *  =====================================================================
00306       SUBROUTINE DGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
00307      $                  BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
00308 *
00309 *  -- LAPACK driver routine (version 3.4.0) --
00310 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00311 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00312 *     November 2011
00313 *
00314 *     .. Scalar Arguments ..
00315       CHARACTER          JOBVL, JOBVR
00316       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
00317 *     ..
00318 *     .. Array Arguments ..
00319       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
00320      $                   B( LDB, * ), BETA( * ), VL( LDVL, * ),
00321      $                   VR( LDVR, * ), WORK( * )
00322 *     ..
00323 *
00324 *  =====================================================================
00325 *
00326 *     .. Parameters ..
00327       DOUBLE PRECISION   ZERO, ONE
00328       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00329 *     ..
00330 *     .. Local Scalars ..
00331       LOGICAL            ILIMIT, ILV, ILVL, ILVR, LQUERY
00332       CHARACTER          CHTEMP
00333       INTEGER            ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
00334      $                   IN, IRIGHT, IROWS, ITAU, IWORK, JC, JR, LOPT,
00335      $                   LWKMIN, LWKOPT, NB, NB1, NB2, NB3
00336       DOUBLE PRECISION   ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM,
00337      $                   BNRM1, BNRM2, EPS, ONEPLS, SAFMAX, SAFMIN,
00338      $                   SALFAI, SALFAR, SBETA, SCALE, TEMP
00339 *     ..
00340 *     .. Local Arrays ..
00341       LOGICAL            LDUMMA( 1 )
00342 *     ..
00343 *     .. External Subroutines ..
00344       EXTERNAL           DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLACPY,
00345      $                   DLASCL, DLASET, DORGQR, DORMQR, DTGEVC, XERBLA
00346 *     ..
00347 *     .. External Functions ..
00348       LOGICAL            LSAME
00349       INTEGER            ILAENV
00350       DOUBLE PRECISION   DLAMCH, DLANGE
00351       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANGE
00352 *     ..
00353 *     .. Intrinsic Functions ..
00354       INTRINSIC          ABS, INT, MAX
00355 *     ..
00356 *     .. Executable Statements ..
00357 *
00358 *     Decode the input arguments
00359 *
00360       IF( LSAME( JOBVL, 'N' ) ) THEN
00361          IJOBVL = 1
00362          ILVL = .FALSE.
00363       ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
00364          IJOBVL = 2
00365          ILVL = .TRUE.
00366       ELSE
00367          IJOBVL = -1
00368          ILVL = .FALSE.
00369       END IF
00370 *
00371       IF( LSAME( JOBVR, 'N' ) ) THEN
00372          IJOBVR = 1
00373          ILVR = .FALSE.
00374       ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
00375          IJOBVR = 2
00376          ILVR = .TRUE.
00377       ELSE
00378          IJOBVR = -1
00379          ILVR = .FALSE.
00380       END IF
00381       ILV = ILVL .OR. ILVR
00382 *
00383 *     Test the input arguments
00384 *
00385       LWKMIN = MAX( 8*N, 1 )
00386       LWKOPT = LWKMIN
00387       WORK( 1 ) = LWKOPT
00388       LQUERY = ( LWORK.EQ.-1 )
00389       INFO = 0
00390       IF( IJOBVL.LE.0 ) THEN
00391          INFO = -1
00392       ELSE IF( IJOBVR.LE.0 ) THEN
00393          INFO = -2
00394       ELSE IF( N.LT.0 ) THEN
00395          INFO = -3
00396       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00397          INFO = -5
00398       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00399          INFO = -7
00400       ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
00401          INFO = -12
00402       ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
00403          INFO = -14
00404       ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
00405          INFO = -16
00406       END IF
00407 *
00408       IF( INFO.EQ.0 ) THEN
00409          NB1 = ILAENV( 1, 'DGEQRF', ' ', N, N, -1, -1 )
00410          NB2 = ILAENV( 1, 'DORMQR', ' ', N, N, N, -1 )
00411          NB3 = ILAENV( 1, 'DORGQR', ' ', N, N, N, -1 )
00412          NB = MAX( NB1, NB2, NB3 )
00413          LOPT = 2*N + MAX( 6*N, N*( NB+1 ) )
00414          WORK( 1 ) = LOPT
00415       END IF
00416 *
00417       IF( INFO.NE.0 ) THEN
00418          CALL XERBLA( 'DGEGV ', -INFO )
00419          RETURN
00420       ELSE IF( LQUERY ) THEN
00421          RETURN
00422       END IF
00423 *
00424 *     Quick return if possible
00425 *
00426       IF( N.EQ.0 )
00427      $   RETURN
00428 *
00429 *     Get machine constants
00430 *
00431       EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
00432       SAFMIN = DLAMCH( 'S' )
00433       SAFMIN = SAFMIN + SAFMIN
00434       SAFMAX = ONE / SAFMIN
00435       ONEPLS = ONE + ( 4*EPS )
00436 *
00437 *     Scale A
00438 *
00439       ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
00440       ANRM1 = ANRM
00441       ANRM2 = ONE
00442       IF( ANRM.LT.ONE ) THEN
00443          IF( SAFMAX*ANRM.LT.ONE ) THEN
00444             ANRM1 = SAFMIN
00445             ANRM2 = SAFMAX*ANRM
00446          END IF
00447       END IF
00448 *
00449       IF( ANRM.GT.ZERO ) THEN
00450          CALL DLASCL( 'G', -1, -1, ANRM, ONE, N, N, A, LDA, IINFO )
00451          IF( IINFO.NE.0 ) THEN
00452             INFO = N + 10
00453             RETURN
00454          END IF
00455       END IF
00456 *
00457 *     Scale B
00458 *
00459       BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
00460       BNRM1 = BNRM
00461       BNRM2 = ONE
00462       IF( BNRM.LT.ONE ) THEN
00463          IF( SAFMAX*BNRM.LT.ONE ) THEN
00464             BNRM1 = SAFMIN
00465             BNRM2 = SAFMAX*BNRM
00466          END IF
00467       END IF
00468 *
00469       IF( BNRM.GT.ZERO ) THEN
00470          CALL DLASCL( 'G', -1, -1, BNRM, ONE, N, N, B, LDB, IINFO )
00471          IF( IINFO.NE.0 ) THEN
00472             INFO = N + 10
00473             RETURN
00474          END IF
00475       END IF
00476 *
00477 *     Permute the matrix to make it more nearly triangular
00478 *     Workspace layout:  (8*N words -- "work" requires 6*N words)
00479 *        left_permutation, right_permutation, work...
00480 *
00481       ILEFT = 1
00482       IRIGHT = N + 1
00483       IWORK = IRIGHT + N
00484       CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
00485      $             WORK( IRIGHT ), WORK( IWORK ), IINFO )
00486       IF( IINFO.NE.0 ) THEN
00487          INFO = N + 1
00488          GO TO 120
00489       END IF
00490 *
00491 *     Reduce B to triangular form, and initialize VL and/or VR
00492 *     Workspace layout:  ("work..." must have at least N words)
00493 *        left_permutation, right_permutation, tau, work...
00494 *
00495       IROWS = IHI + 1 - ILO
00496       IF( ILV ) THEN
00497          ICOLS = N + 1 - ILO
00498       ELSE
00499          ICOLS = IROWS
00500       END IF
00501       ITAU = IWORK
00502       IWORK = ITAU + IROWS
00503       CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
00504      $             WORK( IWORK ), LWORK+1-IWORK, IINFO )
00505       IF( IINFO.GE.0 )
00506      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00507       IF( IINFO.NE.0 ) THEN
00508          INFO = N + 2
00509          GO TO 120
00510       END IF
00511 *
00512       CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
00513      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ),
00514      $             LWORK+1-IWORK, IINFO )
00515       IF( IINFO.GE.0 )
00516      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00517       IF( IINFO.NE.0 ) THEN
00518          INFO = N + 3
00519          GO TO 120
00520       END IF
00521 *
00522       IF( ILVL ) THEN
00523          CALL DLASET( 'Full', N, N, ZERO, ONE, VL, LDVL )
00524          CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
00525      $                VL( ILO+1, ILO ), LDVL )
00526          CALL DORGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
00527      $                WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK,
00528      $                IINFO )
00529          IF( IINFO.GE.0 )
00530      $      LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00531          IF( IINFO.NE.0 ) THEN
00532             INFO = N + 4
00533             GO TO 120
00534          END IF
00535       END IF
00536 *
00537       IF( ILVR )
00538      $   CALL DLASET( 'Full', N, N, ZERO, ONE, VR, LDVR )
00539 *
00540 *     Reduce to generalized Hessenberg form
00541 *
00542       IF( ILV ) THEN
00543 *
00544 *        Eigenvectors requested -- work on whole matrix.
00545 *
00546          CALL DGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
00547      $                LDVL, VR, LDVR, IINFO )
00548       ELSE
00549          CALL DGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
00550      $                B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IINFO )
00551       END IF
00552       IF( IINFO.NE.0 ) THEN
00553          INFO = N + 5
00554          GO TO 120
00555       END IF
00556 *
00557 *     Perform QZ algorithm
00558 *     Workspace layout:  ("work..." must have at least 1 word)
00559 *        left_permutation, right_permutation, work...
00560 *
00561       IWORK = ITAU
00562       IF( ILV ) THEN
00563          CHTEMP = 'S'
00564       ELSE
00565          CHTEMP = 'E'
00566       END IF
00567       CALL DHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
00568      $             ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
00569      $             WORK( IWORK ), LWORK+1-IWORK, IINFO )
00570       IF( IINFO.GE.0 )
00571      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00572       IF( IINFO.NE.0 ) THEN
00573          IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN
00574             INFO = IINFO
00575          ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN
00576             INFO = IINFO - N
00577          ELSE
00578             INFO = N + 6
00579          END IF
00580          GO TO 120
00581       END IF
00582 *
00583       IF( ILV ) THEN
00584 *
00585 *        Compute Eigenvectors  (DTGEVC requires 6*N words of workspace)
00586 *
00587          IF( ILVL ) THEN
00588             IF( ILVR ) THEN
00589                CHTEMP = 'B'
00590             ELSE
00591                CHTEMP = 'L'
00592             END IF
00593          ELSE
00594             CHTEMP = 'R'
00595          END IF
00596 *
00597          CALL DTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
00598      $                VR, LDVR, N, IN, WORK( IWORK ), IINFO )
00599          IF( IINFO.NE.0 ) THEN
00600             INFO = N + 7
00601             GO TO 120
00602          END IF
00603 *
00604 *        Undo balancing on VL and VR, rescale
00605 *
00606          IF( ILVL ) THEN
00607             CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
00608      $                   WORK( IRIGHT ), N, VL, LDVL, IINFO )
00609             IF( IINFO.NE.0 ) THEN
00610                INFO = N + 8
00611                GO TO 120
00612             END IF
00613             DO 50 JC = 1, N
00614                IF( ALPHAI( JC ).LT.ZERO )
00615      $            GO TO 50
00616                TEMP = ZERO
00617                IF( ALPHAI( JC ).EQ.ZERO ) THEN
00618                   DO 10 JR = 1, N
00619                      TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
00620    10             CONTINUE
00621                ELSE
00622                   DO 20 JR = 1, N
00623                      TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
00624      $                      ABS( VL( JR, JC+1 ) ) )
00625    20             CONTINUE
00626                END IF
00627                IF( TEMP.LT.SAFMIN )
00628      $            GO TO 50
00629                TEMP = ONE / TEMP
00630                IF( ALPHAI( JC ).EQ.ZERO ) THEN
00631                   DO 30 JR = 1, N
00632                      VL( JR, JC ) = VL( JR, JC )*TEMP
00633    30             CONTINUE
00634                ELSE
00635                   DO 40 JR = 1, N
00636                      VL( JR, JC ) = VL( JR, JC )*TEMP
00637                      VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
00638    40             CONTINUE
00639                END IF
00640    50       CONTINUE
00641          END IF
00642          IF( ILVR ) THEN
00643             CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
00644      $                   WORK( IRIGHT ), N, VR, LDVR, IINFO )
00645             IF( IINFO.NE.0 ) THEN
00646                INFO = N + 9
00647                GO TO 120
00648             END IF
00649             DO 100 JC = 1, N
00650                IF( ALPHAI( JC ).LT.ZERO )
00651      $            GO TO 100
00652                TEMP = ZERO
00653                IF( ALPHAI( JC ).EQ.ZERO ) THEN
00654                   DO 60 JR = 1, N
00655                      TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
00656    60             CONTINUE
00657                ELSE
00658                   DO 70 JR = 1, N
00659                      TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
00660      $                      ABS( VR( JR, JC+1 ) ) )
00661    70             CONTINUE
00662                END IF
00663                IF( TEMP.LT.SAFMIN )
00664      $            GO TO 100
00665                TEMP = ONE / TEMP
00666                IF( ALPHAI( JC ).EQ.ZERO ) THEN
00667                   DO 80 JR = 1, N
00668                      VR( JR, JC ) = VR( JR, JC )*TEMP
00669    80             CONTINUE
00670                ELSE
00671                   DO 90 JR = 1, N
00672                      VR( JR, JC ) = VR( JR, JC )*TEMP
00673                      VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
00674    90             CONTINUE
00675                END IF
00676   100       CONTINUE
00677          END IF
00678 *
00679 *        End of eigenvector calculation
00680 *
00681       END IF
00682 *
00683 *     Undo scaling in alpha, beta
00684 *
00685 *     Note: this does not give the alpha and beta for the unscaled
00686 *     problem.
00687 *
00688 *     Un-scaling is limited to avoid underflow in alpha and beta
00689 *     if they are significant.
00690 *
00691       DO 110 JC = 1, N
00692          ABSAR = ABS( ALPHAR( JC ) )
00693          ABSAI = ABS( ALPHAI( JC ) )
00694          ABSB = ABS( BETA( JC ) )
00695          SALFAR = ANRM*ALPHAR( JC )
00696          SALFAI = ANRM*ALPHAI( JC )
00697          SBETA = BNRM*BETA( JC )
00698          ILIMIT = .FALSE.
00699          SCALE = ONE
00700 *
00701 *        Check for significant underflow in ALPHAI
00702 *
00703          IF( ABS( SALFAI ).LT.SAFMIN .AND. ABSAI.GE.
00704      $       MAX( SAFMIN, EPS*ABSAR, EPS*ABSB ) ) THEN
00705             ILIMIT = .TRUE.
00706             SCALE = ( ONEPLS*SAFMIN / ANRM1 ) /
00707      $              MAX( ONEPLS*SAFMIN, ANRM2*ABSAI )
00708 *
00709          ELSE IF( SALFAI.EQ.ZERO ) THEN
00710 *
00711 *           If insignificant underflow in ALPHAI, then make the
00712 *           conjugate eigenvalue real.
00713 *
00714             IF( ALPHAI( JC ).LT.ZERO .AND. JC.GT.1 ) THEN
00715                ALPHAI( JC-1 ) = ZERO
00716             ELSE IF( ALPHAI( JC ).GT.ZERO .AND. JC.LT.N ) THEN
00717                ALPHAI( JC+1 ) = ZERO
00718             END IF
00719          END IF
00720 *
00721 *        Check for significant underflow in ALPHAR
00722 *
00723          IF( ABS( SALFAR ).LT.SAFMIN .AND. ABSAR.GE.
00724      $       MAX( SAFMIN, EPS*ABSAI, EPS*ABSB ) ) THEN
00725             ILIMIT = .TRUE.
00726             SCALE = MAX( SCALE, ( ONEPLS*SAFMIN / ANRM1 ) /
00727      $              MAX( ONEPLS*SAFMIN, ANRM2*ABSAR ) )
00728          END IF
00729 *
00730 *        Check for significant underflow in BETA
00731 *
00732          IF( ABS( SBETA ).LT.SAFMIN .AND. ABSB.GE.
00733      $       MAX( SAFMIN, EPS*ABSAR, EPS*ABSAI ) ) THEN
00734             ILIMIT = .TRUE.
00735             SCALE = MAX( SCALE, ( ONEPLS*SAFMIN / BNRM1 ) /
00736      $              MAX( ONEPLS*SAFMIN, BNRM2*ABSB ) )
00737          END IF
00738 *
00739 *        Check for possible overflow when limiting scaling
00740 *
00741          IF( ILIMIT ) THEN
00742             TEMP = ( SCALE*SAFMIN )*MAX( ABS( SALFAR ), ABS( SALFAI ),
00743      $             ABS( SBETA ) )
00744             IF( TEMP.GT.ONE )
00745      $         SCALE = SCALE / TEMP
00746             IF( SCALE.LT.ONE )
00747      $         ILIMIT = .FALSE.
00748          END IF
00749 *
00750 *        Recompute un-scaled ALPHAR, ALPHAI, BETA if necessary.
00751 *
00752          IF( ILIMIT ) THEN
00753             SALFAR = ( SCALE*ALPHAR( JC ) )*ANRM
00754             SALFAI = ( SCALE*ALPHAI( JC ) )*ANRM
00755             SBETA = ( SCALE*BETA( JC ) )*BNRM
00756          END IF
00757          ALPHAR( JC ) = SALFAR
00758          ALPHAI( JC ) = SALFAI
00759          BETA( JC ) = SBETA
00760   110 CONTINUE
00761 *
00762   120 CONTINUE
00763       WORK( 1 ) = LWKOPT
00764 *
00765       RETURN
00766 *
00767 *     End of DGEGV
00768 *
00769       END
 All Files Functions