LAPACK 3.3.0

zggevx.f

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