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