Universidad Complutense de Madrid _ Ingeniería Matemática _ asignatura_ Calculo Científico _Programacion en Fortran y C_ _ LAPACK_ZIP

Esta es una vista previa del archivo. Inicie sesión para ver el archivo original

cgeev.for
 SUBROUTINE CGEEV( JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR,
 $ WORK, LWORK, RWORK, INFO )
*
* -- LAPACK driver routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* June 30, 1999
*
* .. Scalar Arguments ..
 CHARACTER JOBVL, JOBVR
 INTEGER INFO, LDA, LDVL, LDVR, LWORK, N
* ..
* .. Array Arguments ..
 REAL RWORK( * )
 COMPLEX A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
 $ W( * ), WORK( * )
* ..
*
* Purpose
* =======
*
* CGEEV computes for an N-by-N complex nonsymmetric matrix A, the
* eigenvalues and, optionally, the left and/or right eigenvectors.
*
* The right eigenvector v(j) of A satisfies
* A * v(j) = lambda(j) * v(j)
* where lambda(j) is its eigenvalue.
* The left eigenvector u(j) of A satisfies
* u(j)**H * A = lambda(j) * u(j)**H
* where u(j)**H denotes the conjugate transpose of u(j).
*
* The computed eigenvectors are normalized to have Euclidean norm
* equal to 1 and largest component real.
*
* Arguments
* =========
*
* JOBVL (input) CHARACTER*1
* = 'N': left eigenvectors of A are not computed;
* = 'V': left eigenvectors of are computed.
*
* JOBVR (input) CHARACTER*1
* = 'N': right eigenvectors of A are not computed;
* = 'V': right eigenvectors of A are computed.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* A (input/output) COMPLEX array, dimension (LDA,N)
* On entry, the N-by-N matrix A.
* On exit, A has been overwritten.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* W (output) COMPLEX array, dimension (N)
* W contains the computed eigenvalues.
*
* VL (output) COMPLEX array, dimension (LDVL,N)
* If JOBVL = 'V', the left eigenvectors u(j) are stored one
* after another in the columns of VL, in the same order
* as their eigenvalues.
* If JOBVL = 'N', VL is not referenced.
* u(j) = VL(:,j), the j-th column of VL.
*
* LDVL (input) INTEGER
* The leading dimension of the array VL. LDVL >= 1; if
* JOBVL = 'V', LDVL >= N.
*
* VR (output) COMPLEX array, dimension (LDVR,N)
* If JOBVR = 'V', the right eigenvectors v(j) are stored one
* after another in the columns of VR, in the same order
* as their eigenvalues.
* If JOBVR = 'N', VR is not referenced.
* v(j) = VR(:,j), the j-th column of VR.
*
* LDVR (input) INTEGER
* The leading dimension of the array VR. LDVR >= 1; if
* JOBVR = 'V', LDVR >= N.
*
* WORK (workspace/output) COMPLEX array, dimension (LWORK)
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
* LWORK (input) INTEGER
* The dimension of the array WORK. LWORK >= max(1,2*N).
* For good performance, LWORK must generally be larger.
*
* If LWORK = -1, then a workspace query is assumed; the routine
* only calculates the optimal size of the WORK array, returns
* this value as the first entry of the WORK array, and no error
* message related to LWORK is issued by XERBLA.
*
* RWORK (workspace) REAL array, dimension (2*N)
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value.
* > 0: if INFO = i, the QR algorithm failed to compute all the
* eigenvalues, and no eigenvectors have been computed;
* elements and i+1:N of W contain eigenvalues which have
* converged.
*
* =====================================================================
*
* .. Parameters ..
 REAL ZERO, ONE
 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
* ..
* .. Local Scalars ..
 LOGICAL LQUERY, SCALEA, WANTVL, WANTVR
 CHARACTER SIDE
 INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, IRWORK, ITAU,
 $ IWRK, K, MAXB, MAXWRK, MINWRK, NOUT
 REAL ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM
 COMPLEX TMP
* ..
* .. Local Arrays ..
 LOGICAL SELECT( 1 )
 REAL DUM( 1 )
* ..
* .. External Subroutines ..
 EXTERNAL CGEBAK, CGEBAL, CGEHRD, CHSEQR, CLACPY, CLASCL,
 $ CSCAL, CSSCAL, CTREVC, CUNGHR, SLABAD, XERBLA
* ..
* .. External Functions ..
 LOGICAL LSAME
 INTEGER ILAENV, ISAMAX
 REAL CLANGE, SCNRM2, SLAMCH
 EXTERNAL LSAME, ILAENV, ISAMAX, CLANGE, SCNRM2, SLAMCH
* ..
* .. Intrinsic Functions ..
 INTRINSIC AIMAG, CMPLX, CONJG, MAX, MIN, REAL, SQRT
* ..
* .. Executable Statements ..
*
* Test the input arguments
*
 INFO = 0
 LQUERY = ( LWORK.EQ.-1 )
 WANTVL = LSAME( JOBVL, 'V' )
 WANTVR = LSAME( JOBVR, 'V' )
 IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
 INFO = -1
 ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
 INFO = -2
 ELSE IF( N.LT.0 ) THEN
 INFO = -3
 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 INFO = -5
 ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
 INFO = -8
 ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
 INFO = -10
 END IF
*
* Compute workspace
* (Note: Comments in the code beginning "Workspace:" describe the
* minimal amount of workspace needed at that point in the code,
* as well as the preferred amount for good performance.
* CWorkspace refers to complex workspace, and RWorkspace to real
* workspace. NB refers to the optimal block size for the
* immediately following subroutine, as returned by ILAENV.
* HSWORK refers to the workspace preferred by CHSEQR, as
* calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
* the worst case.)
*
 MINWRK = 1
 IF( INFO.EQ.0 .AND. ( LWORK.GE.1 .OR. LQUERY ) ) THEN
 MAXWRK = N + N*ILAENV( 1, 'CGEHRD', ' ', N, 1, N, 0 )
 IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN
 MINWRK = MAX( 1, 2*N )
 MAXB = MAX( ILAENV( 8, 'CHSEQR', 'EN', N, 1, N, -1 ), 2 )
 K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'CHSEQR', 'EN', N, 1,
 $ N, -1 ) ) )
 HSWORK = MAX( K*( K+2 ), 2*N )
 MAXWRK = MAX( MAXWRK, HSWORK )
 ELSE
 MINWRK = MAX( 1, 2*N )
 MAXWRK = MAX( MAXWRK, N+( N-1 )*
 $ ILAENV( 1, 'CUNGHR', ' ', N, 1, N, -1 ) )
 MAXB = MAX( ILAENV( 8, 'CHSEQR', 'SV', N, 1, N, -1 ), 2 )
 K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'CHSEQR', 'SV', N, 1,
 $ N, -1 ) ) )
 HSWORK = MAX( K*( K+2 ), 2*N )
 MAXWRK = MAX( MAXWRK, HSWORK, 2*N )
 END IF
 WORK( 1 ) = MAXWRK
 END IF
 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
 INFO = -12
 END IF
 IF( INFO.NE.0 ) THEN
 CALL XERBLA( 'CGEEV ', -INFO )
 RETURN
 ELSE IF( LQUERY ) THEN
 RETURN
 END IF
*
* Quick return if possible
*
 IF( N.EQ.0 )
 $ RETURN
*
* Get machine constants
*
 EPS = SLAMCH( 'P' )
 SMLNUM = SLAMCH( 'S' )
 BIGNUM = ONE / SMLNUM
 CALL SLABAD( SMLNUM, BIGNUM )
 SMLNUM = SQRT( SMLNUM ) / EPS
 BIGNUM = ONE / SMLNUM
*
* Scale A if max element outside range [SMLNUM,BIGNUM]
*
 ANRM = CLANGE( 'M', N, N,
A, LDA, DUM )
 SCALEA = .FALSE.
 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 SCALEA = .TRUE.
 CSCALE = SMLNUM
 ELSE IF( ANRM.GT.BIGNUM ) THEN
 SCALEA = .TRUE.
 CSCALE = BIGNUM
 END IF
 IF( SCALEA )
 $ CALL CLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
*
* Balance the matrix
* (CWorkspace: none)
* (RWorkspace: need N)
*
 IBAL = 1
 CALL CGEBAL( 'B', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR )
*
* Reduce to upper Hessenberg form
* (CWorkspace: need 2*N, prefer N+N*NB)
* (RWorkspace: none)
*
 ITAU = 1
 IWRK = ITAU + N
 CALL CGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
 $ LWORK-IWRK+1, IERR )
*
 IF( WANTVL ) THEN
*
* Want left eigenvectors
* Copy Householder vectors to VL
*
 SIDE = 'L'
 CALL CLACPY( 'L', N, N, A, LDA, VL, LDVL )
*
* Generate unitary matrix in VL
* (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
* (RWorkspace: none)
*
 CALL CUNGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
 $ LWORK-IWRK+1, IERR )
*
* Perform QR iteration, accumulating Schur vectors in VL
* (CWorkspace: need 1, prefer HSWORK (see comments) )
* (RWorkspace: none)
*
 IWRK = ITAU
 CALL CHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VL, LDVL,
 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
*
 IF( WANTVR ) THEN
*
* Want left and right eigenvectors
* Copy Schur vectors to VR
*
 SIDE = 'B'
 CALL CLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
 END IF
*
 ELSE IF( WANTVR ) THEN
*
* Want right eigenvectors
* Copy Householder vectors to VR
*
 SIDE = 'R'
 CALL CLACPY( 'L', N, N, A, LDA, VR, LDVR )
*
* Generate unitary matrix in VR
* (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
* (RWorkspace: none)
*
 CALL CUNGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
 $ LWORK-IWRK+1, IERR )
*
* Perform QR iteration, accumulating Schur vectors in VR
* (CWorkspace: need 1, prefer HSWORK (see comments) )
* (RWorkspace: none)
*
 IWRK = ITAU
 CALL CHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VR, LDVR,
 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
*
 ELSE
*
* Compute eigenvalues only
* (CWorkspace: need 1, prefer HSWORK (see comments) )
* (RWorkspace: none)
*
 IWRK = ITAU
 CALL CHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, W, VR, LDVR,
 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
 END IF
*
* If INFO > 0 from CHSEQR, then quit
*
 IF( INFO.GT.0 )
 $ GO TO 50
*
 IF( WANTVL .OR. WANTVR ) THEN
*
* Compute left and/or right eigenvectors
* (CWorkspace: need 2*N)
* (RWorkspace: need 2*N)
*
 IRWORK = IBAL + N
 CALL CTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
 $ N, NOUT, WORK( IWRK ), RWORK( IRWORK ), IERR )
 END IF
*
 IF( WANTVL ) THEN
*
* Undo balancing of left eigenvectors
* (CWorkspace: none)
* (RWorkspace: need N)
*
 CALL CGEBAK( 'B', 'L', N, ILO, IHI, RWORK( IBAL ), N, VL, LDVL,
 $ IERR )
*
* Normalize left eigenvectors and make largest component real
*
 DO 20 I = 1, N
 SCL = ONE / SCNRM2( N, VL( 1, I ), 1 )
 CALL CSSCAL( N, SCL, VL( 1, I ), 1 )
 DO 10 K = 1, N
 RWORK( IRWORK+K-1 ) = REAL( VL( K, I ) )**2 +
 $ AIMAG( VL( K, I ) )**2
 10 CONTINUE
 K = ISAMAX( N, RWORK( IRWORK ), 1 )
 TMP = CONJG( VL( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) )
 CALL CSCAL( N, TMP, VL( 1, I ), 1 )
 VL( K, I ) = CMPLX( REAL( VL( K, I ) ), ZERO )
 20 CONTINUE
 END IF
*
 IF( WANTVR ) THEN
*
* Undo balancing of right eigenvectors
* (CWorkspace: none)
* (RWorkspace: need N)
*
 CALL CGEBAK( 'B', 'R', N, ILO, IHI, RWORK( IBAL ), N, VR, LDVR,
 $ IERR )
*
* Normalize right eigenvectors and make largest component real
*
 DO 40 I = 1, N
 SCL = ONE / SCNRM2( N, VR( 1, I ), 1 )
 CALL CSSCAL( N, SCL, VR( 1, I ), 1 )
 DO 30 K = 1, N
 RWORK( IRWORK+K-1 ) = REAL( VR( K, I ) )**2 +
 $ AIMAG( VR( K, I ) )**2
 30 CONTINUE
 K = ISAMAX( N, RWORK( IRWORK ), 1 )
 TMP = CONJG( VR( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) )
 CALL CSCAL( N, TMP, VR( 1, I ), 1 )
 VR( K, I ) = CMPLX( REAL( VR( K, I ) ), ZERO )
 40 CONTINUE
 END IF
*
* Undo scaling if necessary
*
 50 CONTINUE
 IF( SCALEA ) THEN
 CALL CLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, W( INFO+1 ),
 $ MAX( N-INFO, 1 ), IERR )
 IF( INFO.GT.0 ) THEN
 CALL CLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, W, N, IERR )
 END IF
 END IF
*
 WORK( 1 ) = MAXWRK
 RETURN
*
* End of CGEEV
*
 END
cgeevx.for
 SUBROUTINE CGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, W, VL,
 $ LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE,
 $ RCONDV, WORK, LWORK, RWORK, INFO )
*
* -- LAPACK driver routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* June 30, 1999
*
* .. Scalar Arguments ..
 CHARACTER BALANC, JOBVL, JOBVR, SENSE
 INTEGER IHI, ILO, INFO, LDA, LDVL, LDVR, LWORK, N
 REAL ABNRM
* ..
* .. Array Arguments ..
 REAL RCONDE( * ), RCONDV( * ), RWORK( * ),
 $ SCALE( * )
 COMPLEX A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
 $ W( * ), WORK( * )
* ..
*
* Purpose
* =======
*
* CGEEVX computes for an N-by-N complex nonsymmetric matrix A, the
* eigenvalues and, optionally, the left and/or right eigenvectors.
*
* Optionally also, it computes a balancing transformation to improve
* the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
* SCALE, and ABNRM), reciprocal condition numbers for the eigenvalues
* (RCONDE), and reciprocal condition numbers for the right
* eigenvectors (RCONDV).
*
* The right eigenvector v(j) of A satisfies
* A * v(j) = lambda(j) * v(j)
* where lambda(j) is its eigenvalue.
* The left eigenvector u(j) of A satisfies
* u(j)**H * A = lambda(j) * u(j)**H
* where u(j)**H denotes the conjugate transpose of u(j).
*
* The computed eigenvectors are normalized to have Euclidean norm
* equal to 1 and largest component real.
*
* Balancing a matrix means permuting the rows and columns to make it
* more nearly upper triangular, and applying a diagonal similarity
* transformation D * A * D**(-1), where D is a diagonal matrix, to
* make its rows and columns closer in norm and the condition numbers
* of its eigenvalues and eigenvectors smaller. The computed
* reciprocal condition numbers correspond to the balanced matrix.
* Permuting rows and columns will not change the condition numbers
* (in exact arithmetic) but diagonal scaling will. For further
* explanation of balancing, see section 4.10.2 of the LAPACK
* Users' Guide.
*
* Arguments
* =========
*
* BALANC (input) CHARACTER*1
* Indicates how the input matrix should be diagonally scaled
* and/or permuted to improve the conditioning of its
* eigenvalues.
* = 'N': Do not diagonally scale
or permute;
* = 'P': Perform permutations to make the matrix more nearly
* upper triangular. Do not diagonally scale;
* = 'S': Diagonally scale the matrix, ie. replace A by
* D*A*D**(-1), where D is a diagonal matrix chosen
* to make the rows and columns of A more equal in
* norm. Do not permute;
* = 'B': Both diagonally scale and permute A.
*
* Computed reciprocal condition numbers will be for the matrix
* after balancing and/or permuting. Permuting does not change
* condition numbers (in exact arithmetic), but balancing does.
*
* JOBVL (input) CHARACTER*1
* = 'N': left eigenvectors of A are not computed;
* = 'V': left eigenvectors of A are computed.
* If SENSE = 'E' or 'B', JOBVL must = 'V'.
*
* JOBVR (input) CHARACTER*1
* = 'N': right eigenvectors of A are not computed;
* = 'V': right eigenvectors of A are computed.
* If SENSE = 'E' or 'B', JOBVR must = 'V'.
*
* SENSE (input) CHARACTER*1
* Determines which reciprocal condition numbers are computed.
* = 'N': None are computed;
* = 'E': Computed for eigenvalues only;
* = 'V': Computed for right eigenvectors only;
* = 'B': Computed for eigenvalues and right eigenvectors.
*
* If SENSE = 'E' or 'B', both left and right eigenvectors
* must also be computed (JOBVL = 'V' and JOBVR = 'V').
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* A (input/output) COMPLEX array, dimension (LDA,N)
* On entry, the N-by-N matrix A.
* On exit, A has been overwritten. If JOBVL = 'V' or
* JOBVR = 'V', A contains the Schur form of the balanced 
* version of the matrix A.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* W (output) COMPLEX array, dimension (N)
* W contains the computed eigenvalues.
*
* VL (output) COMPLEX array, dimension (LDVL,N)
* If JOBVL = 'V', the left eigenvectors u(j) are stored one
* after another in the columns of VL, in the same order
* as their eigenvalues.
* If JOBVL = 'N', VL is not referenced.
* u(j) = VL(:,j), the j-th column of VL.
*
* LDVL (input) INTEGER
* The leading dimension of the array VL. LDVL >= 1; if
* JOBVL = 'V', LDVL >= N.
*
* VR (output) COMPLEX array, dimension (LDVR,N)
* If JOBVR = 'V', the right eigenvectors v(j) are stored one
* after another in the columns of VR, in the same order
* as their eigenvalues.
* If JOBVR = 'N', VR is not referenced.
* v(j) = VR(:,j), the j-th column of VR.
*
* LDVR (input) INTEGER
* The leading dimension of the array VR. LDVR >= 1; if
* JOBVR = 'V', LDVR >= N.
*
* ILO,IHI (output) INTEGER
* ILO and IHI are integer values determined when A was
* balanced. The balanced A(i,j) = 0 if I > J and
* J = 1,...,ILO-1 or I = IHI+1,...,N.
*
* SCALE (output) REAL array, dimension (N)
* Details of the permutations and scaling factors applied
* when balancing A. If P(j) is the index of the row and column
* interchanged with row and column j, and D(j) is the scaling
* factor applied to row and column j, then
* SCALE(J) = P(J), for J = 1,...,ILO-1
* = D(J), for J = ILO,...,IHI
* = P(J) for J = IHI+1,...,N.
* The order in which the interchanges are made is N to IHI+1,
* then 1 to ILO-1.
*
* ABNRM (output) REAL
* The one-norm of the balanced matrix (the maximum
* of the sum of absolute values of elements of any column).
*
* RCONDE (output) REAL array, dimension (N)
* RCONDE(j) is the reciprocal condition number of the j-th
* eigenvalue.
*
* RCONDV (output) REAL array, dimension (N)
* RCONDV(j) is the reciprocal condition number of the j-th
* right eigenvector.
*
* WORK (workspace/output) COMPLEX array, dimension (LWORK)
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
* LWORK (input) INTEGER
* The dimension of the array WORK. If SENSE = 'N' or 'E',
* LWORK >= max(1,2*N), and if SENSE = 'V' or 'B',
* LWORK >= N*N+2*N.
* For good performance, LWORK must generally be larger.
*
* If LWORK = -1, then a workspace query is assumed; the routine
* only calculates the optimal size of the WORK array, returns
* this value as the first entry of the WORK array, and no error
* message related to LWORK is issued by XERBLA.
*
* RWORK (workspace) REAL array, dimension (2*N)
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value.
* > 0: if INFO = i, the QR algorithm failed to compute all the
* eigenvalues, and no eigenvectors or condition numbers
* have been computed; elements 1:ILO-1 and i+1:N of W
* contain eigenvalues which have converged.
*
* =====================================================================
*
* .. Parameters ..
 REAL ZERO, ONE
 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
* ..
* .. Local Scalars ..
 LOGICAL LQUERY, SCALEA, WANTVL, WANTVR, WNTSNB, WNTSNE,
 $ WNTSNN, WNTSNV
 CHARACTER JOB, SIDE
 INTEGER HSWORK, I, ICOND, IERR, ITAU, IWRK, K, MAXB,
 $ MAXWRK, MINWRK, NOUT
 REAL ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM
 COMPLEX TMP
* ..
* .. Local Arrays ..
 LOGICAL SELECT( 1 )
 REAL DUM( 1 )
* ..
* .. External Subroutines ..
 EXTERNAL CGEBAK, CGEBAL, CGEHRD, CHSEQR, CLACPY, CLASCL,
 $ CSCAL, CSSCAL, CTREVC, CTRSNA, CUNGHR, SLABAD,
 $ SLASCL, XERBLA
* ..
* .. External Functions ..
 LOGICAL LSAME
 INTEGER ILAENV, ISAMAX
 REAL CLANGE, SCNRM2, SLAMCH
 EXTERNAL LSAME, ILAENV, ISAMAX, CLANGE, SCNRM2, SLAMCH
* ..
* .. Intrinsic Functions ..
 INTRINSIC AIMAG, CMPLX, CONJG, MAX, MIN, REAL, SQRT
* ..
* .. Executable Statements ..
*
* Test the input arguments
*
 INFO = 0
 LQUERY = ( LWORK.EQ.-1 )
 WANTVL = LSAME( JOBVL, 'V' )
 WANTVR = LSAME( JOBVR, 'V' )
 WNTSNN = LSAME( SENSE, 'N' )
 WNTSNE = LSAME( SENSE, 'E' )
 WNTSNV = LSAME( SENSE, 'V' )
 WNTSNB = LSAME( SENSE, 'B' )
 IF( .NOT.( LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'S' ) .OR.
 $ LSAME( BALANC, 'P' ) .OR. LSAME( BALANC, 'B' ) ) ) THEN
 INFO = -1
 ELSE IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
 INFO = -2
 ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
 INFO = -3
 ELSE IF( .NOT.( WNTSNN .OR. WNTSNE .OR. WNTSNB .OR. WNTSNV ) .OR.
 $ ( ( WNTSNE .OR. WNTSNB ) .AND. .NOT.( WANTVL .AND.
 $ WANTVR ) ) ) THEN
 INFO = -4
 ELSE IF( N.LT.0 ) THEN
 INFO = -5
 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 INFO = -7
 ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
 INFO = -10
 ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
 INFO = -12
 END IF
*
* Compute workspace
* (Note: Comments in the code beginning "Workspace:" describe the
* minimal amount of workspace needed at that point in the code,
* as well as the
preferred amount for good performance.
* CWorkspace refers to complex workspace, and RWorkspace to real
* workspace. NB refers to the optimal block size for the
* immediately following subroutine, as returned by ILAENV.
* HSWORK refers to the workspace preferred by CHSEQR, as
* calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
* the worst case.)
*
 MINWRK = 1
 IF( INFO.EQ.0 .AND. ( LWORK.GE.1 .OR. LQUERY ) ) THEN
 MAXWRK = N + N*ILAENV( 1, 'CGEHRD', ' ', N, 1, N, 0 )
 IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN
 MINWRK = MAX( 1, 2*N )
 IF( .NOT.( WNTSNN .OR. WNTSNE ) )
 $ MINWRK = MAX( MINWRK, N*N+2*N )
 MAXB = MAX( ILAENV( 8, 'CHSEQR', 'SN', N, 1, N, -1 ), 2 )
 IF( WNTSNN ) THEN
 K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'CHSEQR', 'EN', N,
 $ 1, N, -1 ) ) )
 ELSE
 K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'CHSEQR', 'SN', N,
 $ 1, N, -1 ) ) )
 END IF
 HSWORK = MAX( K*( K+2 ), 2*N )
 MAXWRK = MAX( MAXWRK, 1, HSWORK )
 IF( .NOT.( WNTSNN .OR. WNTSNE ) )
 $ MAXWRK = MAX( MAXWRK, N*N+2*N )
 ELSE
 MINWRK = MAX( 1, 2*N )
 IF( .NOT.( WNTSNN .OR. WNTSNE ) )
 $ MINWRK = MAX( MINWRK, N*N+2*N )
 MAXB = MAX( ILAENV( 8, 'CHSEQR', 'SN', N, 1, N, -1 ), 2 )
 K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'CHSEQR', 'EN', N, 1,
 $ N, -1 ) ) )
 HSWORK = MAX( K*( K+2 ), 2*N )
 MAXWRK = MAX( MAXWRK, 1, HSWORK )
 MAXWRK = MAX( MAXWRK, N+( N-1 )*
 $ ILAENV( 1, 'CUNGHR', ' ', N, 1, N, -1 ) )
 IF( .NOT.( WNTSNN .OR. WNTSNE ) )
 $ MAXWRK = MAX( MAXWRK, N*N+2*N )
 MAXWRK = MAX( MAXWRK, 2*N, 1 )
 END IF
 WORK( 1 ) = MAXWRK
 END IF
 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
 INFO = -20
 END IF
 IF( INFO.NE.0 ) THEN
 CALL XERBLA( 'CGEEVX', -INFO )
 RETURN
 ELSE IF( LQUERY ) THEN
 RETURN
 END IF
*
* Quick return if possible
*
 IF( N.EQ.0 )
 $ RETURN
*
* Get machine constants
*
 EPS = SLAMCH( 'P' )
 SMLNUM = SLAMCH( 'S' )
 BIGNUM = ONE / SMLNUM
 CALL SLABAD( SMLNUM, BIGNUM )
 SMLNUM = SQRT( SMLNUM ) / EPS
 BIGNUM = ONE / SMLNUM
*
* Scale A if max element outside range [SMLNUM,BIGNUM]
*
 ICOND = 0
 ANRM = CLANGE( 'M', N, N, A, LDA, DUM )
 SCALEA = .FALSE.
 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 SCALEA = .TRUE.
 CSCALE = SMLNUM
 ELSE IF( ANRM.GT.BIGNUM ) THEN
 SCALEA = .TRUE.
 CSCALE = BIGNUM
 END IF
 IF( SCALEA )
 $ CALL CLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
*
* Balance the matrix and compute ABNRM
*
 CALL CGEBAL( BALANC, N, A, LDA, ILO, IHI, SCALE, IERR )
 ABNRM = CLANGE( '1', N, N, A, LDA, DUM )
 IF( SCALEA ) THEN
 DUM( 1 ) = ABNRM
 CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR )
 ABNRM = DUM( 1 )
 END IF
*
* Reduce to upper Hessenberg form
* (CWorkspace: need 2*N, prefer N+N*NB)
* (RWorkspace: none)
*
 ITAU = 1
 IWRK = ITAU + N
 CALL CGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
 $ LWORK-IWRK+1, IERR )
*
 IF( WANTVL ) THEN
*
* Want left eigenvectors
* Copy Householder vectors to VL
*
 SIDE = 'L'
 CALL CLACPY( 'L', N, N, A, LDA, VL, LDVL )
*
* Generate unitary matrix in VL
* (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
* (RWorkspace: none)
*
 CALL CUNGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
 $ LWORK-IWRK+1, IERR )
*
* Perform QR iteration, accumulating Schur vectors in VL
* (CWorkspace: need 1, prefer HSWORK (see comments) )
* (RWorkspace: none)
*
 IWRK = ITAU
 CALL CHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VL, LDVL,
 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
*
 IF( WANTVR ) THEN
*
* Want left and right eigenvectors
* Copy Schur vectors to VR
*
 SIDE = 'B'
 CALL CLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
 END IF
*
 ELSE IF( WANTVR ) THEN
*
* Want right eigenvectors
* Copy Householder vectors to VR
*
 SIDE = 'R'
 CALL CLACPY( 'L', N, N, A, LDA, VR, LDVR )
*
* Generate unitary matrix in VR
* (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
* (RWorkspace: none)
*
 CALL CUNGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
 $ LWORK-IWRK+1, IERR )
*
* Perform QR iteration, accumulating Schur vectors in VR
* (CWorkspace: need 1, prefer HSWORK (see comments) )
* (RWorkspace: none)
*
 IWRK = ITAU
 CALL CHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VR, LDVR,
 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
*
 ELSE
*
* Compute eigenvalues only
* If condition numbers desired, compute Schur form
*
 IF( WNTSNN ) THEN
 JOB = 'E'
 ELSE
 JOB = 'S'
 END IF
*
* (CWorkspace: need 1, prefer HSWORK (see comments) )
* (RWorkspace: none)
*
 IWRK = ITAU
 CALL CHSEQR( JOB, 'N', N, ILO, IHI, A, LDA, W, VR, LDVR,
 $ WORK( IWRK ), LWORK-IWRK+1, INFO )
 END IF
*
* If INFO > 0 from CHSEQR, then quit
*
 IF( INFO.GT.0 )
 $ GO TO 50
*
 IF( WANTVL .OR. WANTVR ) THEN
*
* Compute left and/or right eigenvectors
* (CWorkspace: need 2*N)
* (RWorkspace: need N)
*
 CALL CTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
 $ N, NOUT, WORK( IWRK ), RWORK, IERR )
 END IF
*
* Compute condition numbers if desired
* (CWorkspace: need N*N+2*N unless SENSE = 'E')
* (RWorkspace: need 2*N unless SENSE = 'E')
*
 IF( .NOT.WNTSNN ) THEN
 CALL CTRSNA( SENSE, 'A', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
 $ RCONDE, RCONDV, N, NOUT, WORK( IWRK ), N, RWORK,
 $ ICOND )
 END IF
*
 IF( WANTVL ) THEN
*
* Undo balancing of left eigenvectors
*
 CALL CGEBAK( BALANC, 'L', N, ILO, IHI, SCALE, N, VL, LDVL,
 $ IERR )
*
* Normalize left eigenvectors and make largest component real
*
 DO 20 I = 1, N
 SCL = ONE / SCNRM2( N, VL( 1, I ), 1 )
 CALL CSSCAL( N, SCL, VL( 1, I ), 1 )
 DO 10 K = 1, N
 RWORK( K ) = REAL( VL( K, I ) )**2 +
 $ AIMAG( VL( K, I ) )**2
 10 CONTINUE
 K = ISAMAX( N, RWORK, 1 )
 TMP = CONJG( VL( K, I ) ) / SQRT( RWORK( K ) )
 CALL CSCAL( N, TMP, VL( 1, I ), 1 )
 VL( K, I ) = CMPLX( REAL( VL( K, I ) ), ZERO )
 20 CONTINUE
 END IF
*
 IF( WANTVR ) THEN
*
* Undo balancing of right eigenvectors
*
 CALL CGEBAK( BALANC, 'R', N, ILO, IHI, SCALE, N, VR, LDVR,
 $ IERR )
*
* Normalize right eigenvectors and make largest component real
*
 DO 40 I = 1, N
 SCL = ONE / SCNRM2( N, VR( 1, I ), 1 )
 CALL CSSCAL( N, SCL, VR( 1, I ), 1 )
 DO 30 K = 1, N
 RWORK( K ) = REAL( VR( K, I ) )**2 +
 $ AIMAG( VR( K, I ) )**2
 30 CONTINUE
 K = ISAMAX( N, RWORK, 1 )
 TMP = CONJG( VR( K, I ) ) / SQRT( RWORK( K ) )
 CALL CSCAL( N, TMP, VR( 1, I ), 1 )
 VR( K, I ) = CMPLX( REAL( VR( K, I ) ), ZERO )
40 CONTINUE
 END IF
*
* Undo scaling if necessary
*
 50 CONTINUE
 IF( SCALEA ) THEN
 CALL CLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, W( INFO+1 ),
 $ MAX( N-INFO, 1 ), IERR )
 IF( INFO.EQ.0 ) THEN
 IF( ( WNTSNV .OR. WNTSNB ) .AND. ICOND.EQ.0 )
 $ CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, RCONDV, N,
 $ IERR )
 ELSE
 CALL CLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, W, N, IERR )
 END IF
 END IF
*
 WORK( 1 ) = MAXWRK
 RETURN
*
* End of CGEEVX
*
 END
cgegs.for
 SUBROUTINE CGEGS( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHA, BETA,
 $ VSL, LDVSL, VSR, LDVSR, WORK, LWORK, RWORK,
 $ INFO )
*
* -- LAPACK driver routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* June 30, 1999
*
* .. Scalar Arguments ..
 CHARACTER JOBVSL, JOBVSR
 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
* ..
* .. Array Arguments ..
 REAL RWORK( * )
 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ),
 $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
 $ WORK( * )
* ..
*
* Purpose
* =======
*
* This routine is deprecated and has been replaced by routine CGGES.
*
* CGEGS computes for a pair of N-by-N complex nonsymmetric matrices A,
* B: the generalized eigenvalues (alpha, beta), the complex Schur
* form (A, B), and optionally left and/or right Schur vectors
* (VSL and VSR).
*
* (If only the generalized eigenvalues are needed, use the driver CGEGV
* instead.)
*
* A generalized eigenvalue for a pair of matrices (A,B) is, roughly
* speaking, a scalar w or a ratio alpha/beta = w, such that A - w*B
* is singular. It is usually represented as the pair (alpha,beta),
* as there is a reasonable interpretation for beta=0, and even for
* both being zero. A good beginning reference is the book, "Matrix
* Computations", by G. Golub & C. van Loan (Johns Hopkins U. Press)
*
* The (generalized) Schur form of a pair of matrices is the result of
* multiplying both matrices on the left by one unitary matrix and
* both on the right by another unitary matrix, these two unitary
* matrices being chosen so as to bring the pair of matrices into
* upper triangular form with the diagonal elements of B being
* non-negative real numbers (this is also called complex Schur form.)
*
* The left and right Schur vectors are the columns of VSL and VSR,
* respectively, where VSL and VSR are the unitary matrices
* which reduce A and B to Schur form:
*
* Schur form of (A,B) = ( (VSL)**H A (VSR), (VSL)**H B (VSR) )
*
* Arguments
* =========
*
* JOBVSL (input) CHARACTER*1
* = 'N': do not compute the left Schur vectors;
* = 'V': compute the left Schur vectors.
*
* JOBVSR (input) CHARACTER*1
* = 'N': do not compute the right Schur vectors;
* = 'V': compute the right Schur vectors.
*
* N (input) INTEGER
* The order of the matrices A, B, VSL, and VSR. N >= 0.
*
* A (input/output) COMPLEX array, dimension (LDA, N)
* On entry, the first of the pair of matrices whose generalized
* eigenvalues and (optionally) Schur vectors are to be
* computed.
* On exit, the generalized Schur form of A.
*
* LDA (input) INTEGER
* The leading dimension of A. LDA >= max(1,N).
*
* B (input/output) COMPLEX array, dimension (LDB, N)
* On entry, the second of the pair of matrices whose
* generalized eigenvalues and (optionally) Schur vectors are
* to be computed.
* On exit, the generalized Schur form of B.
*
* LDB (input) INTEGER
* The leading dimension of B. LDB >= max(1,N).
*
* ALPHA (output) COMPLEX array, dimension (N)
* BETA (output) COMPLEX array, dimension (N)
* On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
* generalized eigenvalues. ALPHA(j), j=1,...,N and BETA(j),
* j=1,...,N are the diagonals of the complex Schur form (A,B)
* output by CGEGS. The BETA(j) will be non-negative real.
*
* Note: the quotients ALPHA(j)/BETA(j) may easily over- or
* underflow, and BETA(j) may even be zero. Thus, the user
* should avoid naively computing the ratio alpha/beta.
* However, ALPHA will be always less than and usually
* comparable with norm(A) in magnitude, and BETA always less
* than and usually comparable with norm(B).
*
* VSL (output) COMPLEX array, dimension (LDVSL,N)
* If JOBVSL = 'V', VSL will contain the left Schur vectors.
* (See "Purpose", above.)
* Not referenced if JOBVSL = 'N'.
*
* LDVSL (input) INTEGER
* The leading dimension of the matrix VSL. LDVSL >= 1, and
* if JOBVSL = 'V', LDVSL >= N.
*
* VSR (output) COMPLEX array, dimension (LDVSR,N)
* If JOBVSR = 'V', VSR will contain the right Schur vectors.
* (See "Purpose", above.)
* Not referenced if JOBVSR = 'N'.
*
* LDVSR (input) INTEGER
* The leading dimension of the matrix VSR. LDVSR >= 1, and
* if JOBVSR = 'V', LDVSR >= N.
*
* WORK (workspace/output) COMPLEX array, dimension (LWORK)
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
* LWORK (input) INTEGER
* The dimension of the array WORK. LWORK >= max(1,2*N).
* For good performance, LWORK must generally be larger.
* To compute the optimal value of LWORK, call ILAENV to get
* blocksizes (for CGEQRF, CUNMQR, and CUNGQR.) Then compute:
* NB -- MAX of the blocksizes for CGEQRF, CUNMQR, and CUNGQR;
* the optimal LWORK is N*(NB+1).
*
* If LWORK = -1, then a workspace query is assumed; the routine
* only calculates the optimal size of the WORK array, returns
* this value as the first entry of the WORK array, and no error
* message related to LWORK is issued by XERBLA.
*
* RWORK (workspace) REAL array, dimension (3*N)
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value.
* =1,...,N:
* The QZ iteration failed. (A,B) are not in Schur
* form, but ALPHA(j) and BETA(j) should be correct for
* j=INFO+1,...,N.
* > N: errors that usually indicate LAPACK problems:
* =N+1: error return from CGGBAL
* =N+2: error return from CGEQRF
* =N+3: error return from CUNMQR
* =N+4: error return from CUNGQR
* =N+5: error return from CGGHRD
* =N+6: error return from CHGEQZ (other than failed
* iteration)
* =N+7: error return from CGGBAK (computing VSL)
* =N+8: error return from CGGBAK (computing VSR)
* =N+9: error return from CLASCL (various places)
*
* =====================================================================
*
* .. Parameters ..
 REAL ZERO, ONE
 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
 COMPLEX CZERO, CONE
 PARAMETER ( CZERO = ( 0.0E0, 0.0E0 ),
 $ CONE = ( 1.0E0, 0.0E0 ) )
* ..
* .. Local Scalars ..
 LOGICAL ILASCL, ILBSCL, ILVSL, ILVSR, LQUERY
 INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT,
 $ ILO, IRIGHT, IROWS, IRWORK, ITAU, IWORK,
 $ LOPT, LWKMIN, LWKOPT, NB, NB1, NB2,
NB3
 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
 $ SAFMIN, SMLNUM
* ..
* .. External Subroutines ..
 EXTERNAL CGEQRF, CGGBAK, CGGBAL, CGGHRD, CHGEQZ, CLACPY,
 $ CLASCL, CLASET, CUNGQR, CUNMQR, XERBLA
* ..
* .. External Functions ..
 LOGICAL LSAME
 INTEGER ILAENV
 REAL CLANGE, SLAMCH
 EXTERNAL ILAENV, LSAME, CLANGE, SLAMCH
* ..
* .. Intrinsic Functions ..
 INTRINSIC INT, MAX
* ..
* .. Executable Statements ..
*
* Decode the input arguments
*
 IF( LSAME( JOBVSL, 'N' ) ) THEN
 IJOBVL = 1
 ILVSL = .FALSE.
 ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
 IJOBVL = 2
 ILVSL = .TRUE.
 ELSE
 IJOBVL = -1
 ILVSL = .FALSE.
 END IF
*
 IF( LSAME( JOBVSR, 'N' ) ) THEN
 IJOBVR = 1
 ILVSR = .FALSE.
 ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
 IJOBVR = 2
 ILVSR = .TRUE.
 ELSE
 IJOBVR = -1
 ILVSR = .FALSE.
 END IF
*
* Test the input arguments
*
 LWKMIN = MAX( 2*N, 1 )
 LWKOPT = LWKMIN
 WORK( 1 ) = LWKOPT
 LQUERY = ( LWORK.EQ.-1 )
 INFO = 0
 IF( IJOBVL.LE.0 ) THEN
 INFO = -1
 ELSE IF( IJOBVR.LE.0 ) THEN
 INFO = -2
 ELSE IF( N.LT.0 ) THEN
 INFO = -3
 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 INFO = -5
 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 INFO = -7
 ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
 INFO = -11
 ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
 INFO = -13
 ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
 INFO = -15
 END IF
*
 IF( INFO.EQ.0 ) THEN
 NB1 = ILAENV( 1, 'CGEQRF', ' ', N, N, -1, -1 )
 NB2 = ILAENV( 1, 'CUNMQR', ' ', N, N, N, -1 )
 NB3 = ILAENV( 1, 'CUNGQR', ' ', N, N, N, -1 )
 NB = MAX( NB1, NB2, NB3 )
 LOPT = N*(NB+1)
 WORK( 1 ) = LOPT
 END IF
*
 IF( INFO.NE.0 ) THEN
 CALL XERBLA( 'CGEGS ', -INFO )
 RETURN
 ELSE IF( LQUERY ) THEN
 RETURN
 END IF
*
* Quick return if possible
*
 IF( N.EQ.0 )
 $ RETURN
*
* Get machine constants
*
 EPS = SLAMCH( 'E' )*SLAMCH( 'B' )
 SAFMIN = SLAMCH( 'S' )
 SMLNUM = N*SAFMIN / EPS
 BIGNUM = ONE / SMLNUM
*
* Scale A if max element outside range [SMLNUM,BIGNUM]
*
 ANRM = CLANGE( 'M', N, N, A, LDA, RWORK )
 ILASCL = .FALSE.
 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 ANRMTO = SMLNUM
 ILASCL = .TRUE.
 ELSE IF( ANRM.GT.BIGNUM ) THEN
 ANRMTO = BIGNUM
 ILASCL = .TRUE.
 END IF
*
 IF( ILASCL ) THEN
 CALL CLASCL( 'G', -1, -1, ANRM, ANRMTO, N, N, A, LDA, IINFO )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 9
 RETURN
 END IF
 END IF
*
* Scale B if max element outside range [SMLNUM,BIGNUM]
*
 BNRM = CLANGE( 'M', N, N, B, LDB, RWORK )
 ILBSCL = .FALSE.
 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
 BNRMTO = SMLNUM
 ILBSCL = .TRUE.
 ELSE IF( BNRM.GT.BIGNUM ) THEN
 BNRMTO = BIGNUM
 ILBSCL = .TRUE.
 END IF
*
 IF( ILBSCL ) THEN
 CALL CLASCL( 'G', -1, -1, BNRM, BNRMTO, N, N, B, LDB, IINFO )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 9
 RETURN
 END IF
 END IF
*
* Permute the matrix to make it more nearly triangular
*
 ILEFT = 1
 IRIGHT = N + 1
 IRWORK = IRIGHT + N
 IWORK = 1
 CALL CGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
 $ RWORK( IRIGHT ), RWORK( IRWORK ), IINFO )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 1
 GO TO 10
 END IF
*
* Reduce B to triangular form, and initialize VSL and/or VSR
*
 IROWS = IHI + 1 - ILO
 ICOLS = N + 1 - ILO
 ITAU = IWORK
 IWORK = ITAU + IROWS
 CALL CGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
 $ WORK( IWORK ), LWORK+1-IWORK, IINFO )
 IF( IINFO.GE.0 )
 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 2
 GO TO 10
 END IF
*
 CALL CUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ),
 $ LWORK+1-IWORK, IINFO )
 IF( IINFO.GE.0 )
 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 3
 GO TO 10
 END IF
*
 IF( ILVSL ) THEN
 CALL CLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL )
 CALL CLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
 $ VSL( ILO+1, ILO ), LDVSL )
 CALL CUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
 $ WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK,
 $ IINFO )
 IF( IINFO.GE.0 )
 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 4
 GO TO 10
 END IF
 END IF
*
 IF( ILVSR )
 $ CALL CLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR )
*
* Reduce to generalized Hessenberg form
*
 CALL CGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
 $ LDVSL, VSR, LDVSR, IINFO )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 5
 GO TO 10
 END IF
*
* Perform QZ algorithm, computing Schur vectors if desired
*
 IWORK = ITAU
 CALL CHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
 $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWORK ),
 $ LWORK+1-IWORK, RWORK( IRWORK ), IINFO )
 IF( IINFO.GE.0 )
 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
 IF( IINFO.NE.0 ) THEN
 IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN
 INFO = IINFO
 ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN
 INFO = IINFO - N
 ELSE
 INFO = N + 6
 END IF
 GO TO 10
 END IF
*
* Apply permutation to VSL and VSR
*
 IF( ILVSL ) THEN
 CALL CGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
 $ RWORK( IRIGHT ), N, VSL, LDVSL, IINFO )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 7
 GO TO 10
 END IF
 END IF
 IF( ILVSR ) THEN
 CALL CGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
 $ RWORK( IRIGHT ), N, VSR, LDVSR, IINFO )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 8
 GO TO 10
 END IF
 END IF
*
* Undo scaling
*
 IF( ILASCL ) THEN
 CALL CLASCL( 'U', -1, -1, ANRMTO, ANRM, N, N, A, LDA, IINFO )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 9
 RETURN
 END IF
 CALL CLASCL( 'G', -1, -1, ANRMTO, ANRM, N, 1, ALPHA, N, IINFO )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 9
 RETURN
 END IF
 END IF
*
 IF( ILBSCL ) THEN
 CALL CLASCL( 'U', -1, -1, BNRMTO, BNRM, N, N, B, LDB, IINFO )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 9
 RETURN
 END IF
 CALL CLASCL( 'G', -1, -1, BNRMTO, BNRM, N, 1, BETA, N, IINFO )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 9
 RETURN
 END IF
 END IF
*
 10 CONTINUE
 WORK( 1 ) = LWKOPT
*
 RETURN
*
* End of CGEGS
*
 END
cgegv.for
 SUBROUTINE CGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
 $ VL, LDVL, VR, LDVR, WORK, LWORK, RWORK,
INFO )
*
* -- LAPACK driver routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* June 30, 1999
*
* .. Scalar Arguments ..
 CHARACTER JOBVL, JOBVR
 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
* ..
* .. Array Arguments ..
 REAL RWORK( * )
 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ),
 $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
 $ WORK( * )
* ..
*
* Purpose
* =======
*
* This routine is deprecated and has been replaced by routine CGGEV.
*
* CGEGV computes for a pair of N-by-N complex nonsymmetric matrices A
* and B, the generalized eigenvalues (alpha, beta), and optionally,
* the left and/or right generalized eigenvectors (VL and VR).
*
* A generalized eigenvalue for a pair of matrices (A,B) is, roughly
* speaking, a scalar w or a ratio alpha/beta = w, such that A - w*B
* is singular. It is usually represented as the pair (alpha,beta),
* as there is a reasonable interpretation for beta=0, and even for
* both being zero. A good beginning reference is the book, "Matrix
* Computations", by G. Golub & C. van Loan (Johns Hopkins U. Press)
*
* A right generalized eigenvector corresponding to a generalized
* eigenvalue w for a pair of matrices (A,B) is a vector r such
* that (A - w B) r = 0 . A left generalized eigenvector is a vector
* l such that l**H * (A - w B) = 0, where l**H is the
* conjugate-transpose of l.
*
* Note: this routine performs "full balancing" on A and B -- see
* "Further Details", below.
*
* Arguments
* =========
*
* JOBVL (input) CHARACTER*1
* = 'N': do not compute the left generalized eigenvectors;
* = 'V': compute the left generalized eigenvectors.
*
* JOBVR (input) CHARACTER*1
* = 'N': do not compute the right generalized eigenvectors;
* = 'V': compute the right generalized eigenvectors.
*
* N (input) INTEGER
* The order of the matrices A, B, VL, and VR. N >= 0.
*
* A (input/output) COMPLEX array, dimension (LDA, N)
* On entry, the first of the pair of matrices whose
* generalized eigenvalues and (optionally) generalized
* eigenvectors are to be computed.
* On exit, the contents will have been destroyed. (For a
* description of the contents of A on exit, see "Further
* Details", below.)
*
* LDA (input) INTEGER
* The leading dimension of A. LDA >= max(1,N).
*
* B (input/output) COMPLEX array, dimension (LDB, N)
* On entry, the second of the pair of matrices whose
* generalized eigenvalues and (optionally) generalized
* eigenvectors are to be computed.
* On exit, the contents will have been destroyed. (For a
* description of the contents of B on exit, see "Further
* Details", below.)
*
* LDB (input) INTEGER
* The leading dimension of B. LDB >= max(1,N).
*
* ALPHA (output) COMPLEX array, dimension (N)
* BETA (output) COMPLEX array, dimension (N)
* On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
* generalized eigenvalues.
*
* Note: the quotients ALPHA(j)/BETA(j) may easily over- or
* underflow, and BETA(j) may even be zero. Thus, the user
* should avoid naively computing the ratio alpha/beta.
* However, ALPHA will be always less than and usually
* comparable with norm(A) in magnitude, and BETA always less
* than and usually comparable with norm(B).
*
* VL (output) COMPLEX array, dimension (LDVL,N)
* If JOBVL = 'V', the left generalized eigenvectors. (See
* "Purpose", above.)
* Each eigenvector will be scaled so the largest component
* will have abs(real part) + abs(imag. part) = 1, *except*
* that for eigenvalues with alpha=beta=0, a zero vector will
* be returned as the corresponding eigenvector.
* Not referenced if JOBVL = 'N'.
*
* LDVL (input) INTEGER
* The leading dimension of the matrix VL. LDVL >= 1, and
* if JOBVL = 'V', LDVL >= N.
*
* VR (output) COMPLEX array, dimension (LDVR,N)
* If JOBVR = 'V', the right generalized eigenvectors. (See
* "Purpose", above.)
* Each eigenvector will be scaled so the largest component
* will have abs(real part) + abs(imag. part) = 1, *except*
* that for eigenvalues with alpha=beta=0, a zero vector will
* be returned as the corresponding eigenvector.
* Not referenced if JOBVR = 'N'.
*
* LDVR (input) INTEGER
* The leading dimension of the matrix VR. LDVR >= 1, and
* if JOBVR = 'V', LDVR >= N.
*
* WORK (workspace/output) COMPLEX array, dimension (LWORK)
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
* LWORK (input) INTEGER
* The dimension of the array WORK. LWORK >= max(1,2*N).
* For good performance, LWORK must generally be larger.
* To compute the optimal value of LWORK, call ILAENV to get
* blocksizes (for CGEQRF, CUNMQR, and CUNGQR.) Then compute:
* NB -- MAX of the blocksizes for CGEQRF, CUNMQR, and CUNGQR;
* The optimal LWORK is MAX( 2*N, N*(NB+1) ).
*
* If LWORK = -1, then a workspace query is assumed; the routine
* only calculates the optimal size of the WORK array, returns
* this value as the first entry of the WORK array, and no error
* message related to LWORK is issued by XERBLA.
*
* RWORK (workspace/output) REAL array, dimension (8*N)
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value.
* =1,...,N:
* The QZ iteration failed. No eigenvectors have been
* calculated, but ALPHA(j) and BETA(j) should be
* correct for j=INFO+1,...,N.
* > N: errors that usually indicate LAPACK problems:
* =N+1: error return from CGGBAL
* =N+2: error return from CGEQRF
* =N+3: error return from CUNMQR
* =N+4: error return from CUNGQR
* =N+5: error return from CGGHRD
* =N+6: error return from CHGEQZ (other than failed
* iteration)
* =N+7: error return from CTGEVC
* =N+8: error return from CGGBAK (computing VL)
* =N+9: error return from CGGBAK (computing VR)
* =N+10: error return from CLASCL (various calls)
*
* Further Details
* ===============
*
* Balancing
* ---------
*
* This driver calls CGGBAL to both permute and scale rows and columns
* of A and B. The permutations PL and PR are chosen so that PL*A*PR
* and PL*B*R will be upper triangular except for the diagonal blocks
* A(i:j,i:j) and B(i:j,i:j), with i and j as close together as
* possible. The diagonal scaling matrices DL and DR are chosen so
* that the pair DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to
* one (except for the elements that start out zero.)
*
* After the eigenvalues and eigenvectors of the balanced matrices
* have been computed, CGGBAK transforms the eigenvectors back to what
* they would have been (in perfect arithmetic) if they had not been
* balanced.
*
* Contents of A and B on Exit
* -------- -- - --- - -- ----
*
* If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or
* both), then on exit the arrays A and B will contain the complex Schur
* form[*] of the "balanced" versions of A and B. If no eigenvectors
* are computed, then only the diagonal blocks will be correct.
*
* [*] In other words, upper triangular form.
*
* =====================================================================
*
* .. Parameters ..
 REAL ZERO, ONE
 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
 COMPLEX CZERO, CONE
 PARAMETER ( CZERO = ( 0.0E0, 0.0E0 ),
 $ CONE = ( 1.0E0, 0.0E0 ) )
* ..
* .. Local Scalars ..
 LOGICAL ILIMIT, ILV, ILVL, ILVR, LQUERY
 CHARACTER CHTEMP
 INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
 $ IN, IRIGHT, IROWS, IRWORK, ITAU, IWORK, JC, JR,
 $ LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3
 REAL ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM,
 $ BNRM1, BNRM2, EPS, SAFMAX, SAFMIN, SALFAI,
 $ SALFAR, SBETA, SCALE, TEMP
 COMPLEX X
* ..
* .. Local Arrays ..
 LOGICAL LDUMMA( 1 )
* ..
* .. External Subroutines ..
 EXTERNAL CGEQRF, CGGBAK, CGGBAL, CGGHRD, CHGEQZ, CLACPY,
 $ CLASCL, CLASET, CTGEVC, CUNGQR, CUNMQR, XERBLA
* ..
* .. External Functions ..
 LOGICAL LSAME
 INTEGER ILAENV
 REAL CLANGE, SLAMCH
 EXTERNAL ILAENV, LSAME, CLANGE, SLAMCH
* ..
* .. Intrinsic Functions ..
 INTRINSIC ABS, AIMAG, CMPLX, INT, MAX, REAL
* ..
* .. Statement Functions ..
 REAL ABS1
* ..
* .. Statement Function definitions ..
 ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
* ..
* .. Executable Statements ..
*
* Decode the input arguments
*
 IF( LSAME( JOBVL, 'N' ) ) THEN
 IJOBVL = 1
 ILVL = .FALSE.
 ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
 IJOBVL = 2
 ILVL = .TRUE.
 ELSE
 IJOBVL = -1
 ILVL = .FALSE.
 END IF
*
 IF( LSAME( JOBVR, 'N' ) ) THEN
 IJOBVR = 1
 ILVR = .FALSE.
 ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
 IJOBVR = 2
 ILVR = .TRUE.
 ELSE
 IJOBVR = -1
 ILVR = .FALSE.
 END IF
 ILV = ILVL .OR. ILVR
*
* Test the input arguments
*
 LWKMIN = MAX( 2*N, 1 )
 LWKOPT = LWKMIN
 WORK( 1 ) = LWKOPT
 LQUERY = ( LWORK.EQ.-1 )
 INFO = 0
 IF( IJOBVL.LE.0 ) THEN
 INFO = -1
 ELSE IF( IJOBVR.LE.0 ) THEN
 INFO = -2
 ELSE IF( N.LT.0 ) THEN
 INFO = -3
 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 INFO = -5
 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 INFO = -7
 ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
 INFO = -11
 ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
 INFO = -13
 ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
 INFO = -15
 END IF
*
 IF( INFO.EQ.0 ) THEN
 NB1 = ILAENV( 1, 'CGEQRF', ' ', N, N, -1, -1 )
 NB2 = ILAENV( 1, 'CUNMQR', ' ', N, N, N, -1 )
 NB3 = ILAENV( 1, 'CUNGQR', ' ', N, N, N, -1 )
 NB = MAX( NB1, NB2, NB3 )
 LOPT = MAX( 2*N, N*(NB+1) )
 WORK( 1 ) = LOPT
 END IF
*
 IF( INFO.NE.0 ) THEN
 CALL XERBLA( 'CGEGV ', -INFO )
 RETURN
 ELSE IF( LQUERY ) THEN
 RETURN
 END IF
*
* Quick return if possible
*
 IF( N.EQ.0 )
 $ RETURN
*
* Get machine constants
*
 EPS = SLAMCH( 'E' )*SLAMCH( 'B' )
 SAFMIN = SLAMCH( 'S' )
 SAFMIN = SAFMIN + SAFMIN
 SAFMAX = ONE / SAFMIN
*
* Scale A
*
 ANRM = CLANGE( 'M', N, N, A, LDA, RWORK )
 ANRM1 = ANRM
 ANRM2 = ONE
 IF( ANRM.LT.ONE ) THEN
 IF( SAFMAX*ANRM.LT.ONE ) THEN
 ANRM1 = SAFMIN
 ANRM2 = SAFMAX*ANRM
 END IF
 END IF
*
 IF( ANRM.GT.ZERO ) THEN
 CALL CLASCL( 'G', -1, -1, ANRM, ONE, N, N, A, LDA, IINFO )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 10
 RETURN
 END IF
 END IF
*
* Scale B
*
 BNRM = CLANGE( 'M', N, N, B, LDB, RWORK )
 BNRM1 = BNRM
 BNRM2 = ONE
 IF( BNRM.LT.ONE ) THEN
 IF( SAFMAX*BNRM.LT.ONE ) THEN
 BNRM1 = SAFMIN
 BNRM2 = SAFMAX*BNRM
 END IF
 END IF
*
 IF( BNRM.GT.ZERO ) THEN
 CALL CLASCL( 'G', -1, -1, BNRM, ONE, N, N, B, LDB, IINFO )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 10
 RETURN
 END IF
 END IF
*
* Permute the matrix to make it more nearly triangular
* Also "balance" the matrix.
*
 ILEFT = 1
 IRIGHT = N + 1
 IRWORK = IRIGHT + N
 CALL CGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
 $ RWORK( IRIGHT ), RWORK( IRWORK ), IINFO )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 1
 GO TO 80
 END IF
*
* Reduce B to triangular form, and initialize VL and/or VR
*
 IROWS = IHI + 1 - ILO
 IF( ILV ) THEN
 ICOLS = N + 1 - ILO
 ELSE
 ICOLS = IROWS
 END IF
 ITAU = 1
 IWORK = ITAU + IROWS
 CALL CGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
 $ WORK( IWORK ), LWORK+1-IWORK, IINFO )
 IF( IINFO.GE.0 )
 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 2
 GO TO 80
 END IF
*
 CALL CUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ),
 $ LWORK+1-IWORK, IINFO )
 IF( IINFO.GE.0 )
 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 3
 GO TO 80
 END IF
*
 IF( ILVL ) THEN
 CALL CLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
 CALL CLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
 $ VL( ILO+1, ILO ), LDVL )
 CALL CUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
 $ WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK,
 $ IINFO )
 IF( IINFO.GE.0 )
 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 4
 GO TO 80
 END IF
 END IF
*
 IF( ILVR )
 $ CALL CLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
*
* Reduce to generalized Hessenberg form
*
 IF( ILV ) THEN
*
* Eigenvectors requested -- work on whole matrix.
*
 CALL CGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
 $ LDVL, VR, LDVR, IINFO )
 ELSE
 CALL CGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
 $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IINFO )
 END IF
 IF( IINFO.NE.0 ) THEN
 INFO = N + 5
 GO TO 80
 END IF
*
* Perform QZ algorithm
*
 IWORK = ITAU
 IF( ILV ) THEN
 CHTEMP = 'S'
 ELSE
 CHTEMP = 'E'
 END IF
 CALL CHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
 $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWORK ),
 $ LWORK+1-IWORK, RWORK( IRWORK ), IINFO )
 IF( IINFO.GE.0 )
 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
 IF( IINFO.NE.0 ) THEN
 IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN
 INFO = IINFO
 ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN
 INFO = IINFO - N
 ELSE
 INFO = N + 6
 END IF
 GO TO 80
 END IF
*
 IF( ILV ) THEN
*
* Compute Eigenvectors
*
 IF( ILVL ) THEN
 IF( ILVR ) THEN
 CHTEMP = 'B'
 ELSE
 CHTEMP = 'L'
END IF
 ELSE
 CHTEMP = 'R'
 END IF
*
 CALL CTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
 $ VR, LDVR, N, IN, WORK( IWORK ), RWORK( IRWORK ),
 $ IINFO )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 7
 GO TO 80
 END IF
*
* Undo balancing on VL and VR, rescale
*
 IF( ILVL ) THEN
 CALL CGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
 $ RWORK( IRIGHT ), N, VL, LDVL, IINFO )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 8
 GO TO 80
 END IF
 DO 30 JC = 1, N
 TEMP = ZERO
 DO 10 JR = 1, N
 TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
 10 CONTINUE
 IF( TEMP.LT.SAFMIN )
 $ GO TO 30
 TEMP = ONE / TEMP
 DO 20 JR = 1, N
 VL( JR, JC ) = VL( JR, JC )*TEMP
 20 CONTINUE
 30 CONTINUE
 END IF
 IF( ILVR ) THEN
 CALL CGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
 $ RWORK( IRIGHT ), N, VR, LDVR, IINFO )
 IF( IINFO.NE.0 ) THEN
 INFO = N + 9
 GO TO 80
 END IF
 DO 60 JC = 1, N
 TEMP = ZERO
 DO 40 JR = 1, N
 TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
 40 CONTINUE
 IF( TEMP.LT.SAFMIN )
 $ GO TO 60
 TEMP = ONE / TEMP
 DO 50 JR = 1, N
 VR( JR, JC ) = VR( JR, JC )*TEMP
 50 CONTINUE
 60 CONTINUE
 END IF
*
* End of eigenvector calculation
*
 END IF
*
* Undo scaling in alpha, beta
*
* Note: this does not give the alpha and beta for the unscaled
* problem.
*
* Un-scaling is limited to avoid underflow in alpha and beta
* if they are significant.
*
 DO 70 JC = 1, N
 ABSAR = ABS( REAL( ALPHA( JC ) ) )
 ABSAI = ABS( AIMAG( ALPHA( JC ) ) )
 ABSB = ABS( REAL( BETA( JC ) ) )
 SALFAR = ANRM*REAL( ALPHA( JC ) )
 SALFAI = ANRM*AIMAG( ALPHA( JC ) )
 SBETA = BNRM*REAL( BETA( JC ) )
 ILIMIT = .FALSE.
 SCALE = ONE
*
* Check for significant underflow in imaginary part of ALPHA
*
 IF( ABS( SALFAI ).LT.SAFMIN .AND. ABSAI.GE.
 $ MAX( SAFMIN, EPS*ABSAR, EPS*ABSB ) ) THEN
 ILIMIT = .TRUE.
 SCALE = ( SAFMIN / ANRM1 ) / MAX( SAFMIN, ANRM2*ABSAI )
 END IF
*
* Check for significant underflow in real part of ALPHA
*
 IF( ABS( SALFAR ).LT.SAFMIN .AND. ABSAR.GE.
 $ MAX( SAFMIN, EPS*ABSAI, EPS*ABSB ) ) THEN
 ILIMIT = .TRUE.
 SCALE = MAX( SCALE, ( SAFMIN / ANRM1 ) /
 $ MAX( SAFMIN, ANRM2*ABSAR ) )
 END IF
*
* Check for significant underflow in BETA
*
 IF( ABS( SBETA ).LT.SAFMIN .AND. ABSB.GE.
 $ MAX( SAFMIN, EPS*ABSAR, EPS*ABSAI ) ) THEN
 ILIMIT = .TRUE.
 SCALE = MAX( SCALE, ( SAFMIN / BNRM1 ) /
 $ MAX( SAFMIN, BNRM2*ABSB ) )
 END IF
*
* Check for possible overflow when limiting scaling
*
 IF( ILIMIT ) THEN
 TEMP = ( SCALE*SAFMIN )*MAX( ABS( SALFAR ), ABS( SALFAI ),
 $ ABS( SBETA ) )
 IF( TEMP.GT.ONE )
 $ SCALE = SCALE / TEMP
 IF( SCALE.LT.ONE )
 $ ILIMIT = .FALSE.
 END IF
*
* Recompute un-scaled ALPHA, BETA if necessary.
*
 IF( ILIMIT ) THEN
 SALFAR = ( SCALE*REAL( ALPHA( JC ) ) )*ANRM
 SALFAI = ( SCALE*AIMAG( ALPHA( JC ) ) )*ANRM
 SBETA = ( SCALE*BETA( JC ) )*BNRM
 END IF
 ALPHA( JC ) = CMPLX( SALFAR, SALFAI )
 BETA( JC ) = SBETA
 70 CONTINUE
*
 80 CONTINUE
 WORK( 1 ) = LWKOPT
*
 RETURN
*
* End of CGEGV
*
 END
cgehd2.for
 SUBROUTINE CGEHD2( N, ILO, IHI, A, LDA, TAU, WORK, INFO )
*
* -- LAPACK routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* September 30, 1994
*
* .. Scalar Arguments ..
 INTEGER IHI, ILO, INFO, LDA, N
* ..
* .. Array Arguments ..
 COMPLEX A( LDA, * ), TAU( * ), WORK( * )
* ..
*
* Purpose
* =======
*
* CGEHD2 reduces a complex general matrix A to upper Hessenberg form H
* by a unitary similarity transformation: Q' * A * Q = H .
*
* Arguments
* =========
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* ILO (input) INTEGER
* IHI (input) INTEGER
* It is assumed that A is already upper triangular in rows
* and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
* set by a previous call to CGEBAL; otherwise they should be
* set to 1 and N respectively. See Further Details.
* 1 <= ILO <= IHI <= max(1,N).
*
* A (input/output) COMPLEX array, dimension (LDA,N)
* On entry, the n by n general matrix to be reduced.
* On exit, the upper triangle and the first subdiagonal of A
* are overwritten with the upper Hessenberg matrix H, and the
* elements below the first subdiagonal, with the array TAU,
* represent the unitary matrix Q as a product of elementary
* reflectors. See Further Details.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* TAU (output) COMPLEX array, dimension (N-1)
* The scalar factors of the elementary reflectors (see Further
* Details).
*
* WORK (workspace) COMPLEX array, dimension (N)
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value.
*
* Further Details
* ===============
*
* The matrix Q is represented as a product of (ihi-ilo) elementary
* reflectors
*
* Q = H(ilo) H(ilo+1) . . . H(ihi-1).
*
* Each H(i) has the form
*
* H(i) = I - tau * v * v'
*
* where tau is a complex scalar, and v is a complex vector with
* v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on
* exit in A(i+2:ihi,i), and tau in TAU(i).
*
* The contents of A are illustrated by the following example, with
* n = 7, ilo = 2 and ihi = 6:
*
* on entry, on exit,
*
* ( a a a a a a a ) ( a a h h h h a )
* ( a a a a a a ) ( a h h h h a )
* ( a a a a a a ) ( h h h h h h )
* ( a a a a a a ) ( v2 h h h h h )
* ( a a a a a a ) ( v2 v3 h h h h )
* ( a a a a a a ) ( v2 v3 v4 h h h )
* ( a ) ( a )
*
* where a denotes an element of the original matrix A, h denotes a
* modified element of the upper Hessenberg matrix H, and vi denotes an
* element of the vector defining H(i).
*
* =====================================================================
*
* .. Parameters ..
 COMPLEX ONE
 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) )
* ..
* .. Local Scalars ..
 INTEGER I
 COMPLEX ALPHA
* ..
* .. External Subroutines ..
 EXTERNAL CLARF, CLARFG, XERBLA
* ..
* .. Intrinsic Functions ..
 INTRINSIC CONJG, MAX, MIN
* ..
* .. Executable Statements ..
*
* Test the input parameters
*
 INFO = 0
 IF( N.LT.0 ) THEN
 INFO = -1
 ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX(
1, N ) ) THEN
 INFO = -2
 ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN
 INFO = -3
 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 INFO = -5
 END IF
 IF( INFO.NE.0 ) THEN
 CALL XERBLA( 'CGEHD2', -INFO )
 RETURN
 END IF
*
 DO 10 I = ILO, IHI - 1
*
* Compute elementary reflector H(i) to annihilate A(i+2:ihi,i)
*
 ALPHA = A( I+1, I )
 CALL CLARFG( IHI-I, ALPHA, A( MIN( I+2, N ), I ), 1, TAU( I ) )
 A( I+1, I ) = ONE
*
* Apply H(i) to A(1:ihi,i+1:ihi) from the right
*
 CALL CLARF( 'Right', IHI, IHI-I, A( I+1, I ), 1, TAU( I ),
 $ A( 1, I+1 ), LDA, WORK )
*
* Apply H(i)' to A(i+1:ihi,i+1:n) from the left
*
 CALL CLARF( 'Left', IHI-I, N-I, A( I+1, I ), 1,
 $ CONJG( TAU( I ) ), A( I+1, I+1 ), LDA, WORK )
*
 A( I+1, I ) = ALPHA
 10 CONTINUE
*
 RETURN
*
* End of CGEHD2
*
 END
cgehrd.for
 SUBROUTINE CGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
*
* -- LAPACK routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* June 30, 1999
*
* .. Scalar Arguments ..
 INTEGER IHI, ILO, INFO, LDA, LWORK, N
* ..
* .. Array Arguments ..
 COMPLEX A( LDA, * ), TAU( * ), WORK( * )
* ..
*
* Purpose
* =======
*
* CGEHRD reduces a complex general matrix A to upper Hessenberg form H
* by a unitary similarity transformation: Q' * A * Q = H .
*
* Arguments
* =========
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* ILO (input) INTEGER
* IHI (input) INTEGER
* It is assumed that A is already upper triangular in rows
* and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
* set by a previous call to CGEBAL; otherwise they should be
* set to 1 and N respectively. See Further Details.
* 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
*
* A (input/output) COMPLEX array, dimension (LDA,N)
* On entry, the N-by-N general matrix to be reduced.
* On exit, the upper triangle and the first subdiagonal of A
* are overwritten with the upper Hessenberg matrix H, and the
* elements below the first subdiagonal, with the array TAU,
* represent the unitary matrix Q as a product of elementary
* reflectors. See Further Details.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* TAU (output) COMPLEX array, dimension (N-1)
* The scalar factors of the elementary reflectors (see Further
* Details). Elements 1:ILO-1 and IHI:N-1 of TAU are set to
* zero.
*
* WORK (workspace/output) COMPLEX array, dimension (LWORK)
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
* LWORK (input) INTEGER
* The length of the array WORK. LWORK >= max(1,N).
* For optimum performance LWORK >= N*NB, where NB is the
* optimal blocksize.
*
* If LWORK = -1, then a workspace query is assumed; the routine
* only calculates the optimal size of the WORK array, returns
* this value as the first entry of the WORK array, and no error
* message related to LWORK is issued by XERBLA.
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value.
*
* Further Details
* ===============
*
* The matrix Q is represented as a product of (ihi-ilo) elementary
* reflectors
*
* Q = H(ilo) H(ilo+1) . . . H(ihi-1).
*
* Each H(i) has the form
*
* H(i) = I - tau * v * v'
*
* where tau is a complex scalar, and v is a complex vector with
* v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on
* exit in A(i+2:ihi,i), and tau in TAU(i).
*
* The contents of A are illustrated by the following example, with
* n = 7, ilo = 2 and ihi = 6:
*
* on entry, on exit,
*
* ( a a a a a a a ) ( a a h h h h a )
* ( a a a a a a ) ( a h h h h a )
* ( a a a a a a ) ( h h h h h h )
* ( a a a a a a ) ( v2 h h h h h )
* ( a a a a a a ) ( v2 v3 h h h h )
* ( a a a a a a ) ( v2 v3 v4 h h h )
* ( a ) ( a )
*
* where a denotes an element of the original matrix A, h denotes a
* modified element of the upper Hessenberg matrix H, and vi denotes an
* element of the vector defining H(i).
*
* =====================================================================
*
* .. Parameters ..
 INTEGER NBMAX, LDT
 PARAMETER ( NBMAX = 64, LDT = NBMAX+1 )
 COMPLEX ZERO, ONE
 PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ),
 $ ONE = ( 1.0E+0, 0.0E+0 ) )
* ..
* .. Local Scalars ..
 LOGICAL LQUERY
 INTEGER I, IB, IINFO, IWS, LDWORK, LWKOPT, NB, NBMIN,
 $ NH, NX
 COMPLEX EI
* ..
* .. Local Arrays ..
 COMPLEX T( LDT, NBMAX )
* ..
* .. External Subroutines ..
 EXTERNAL CGEHD2, CGEMM, CLAHRD, CLARFB, XERBLA
* ..
* .. Intrinsic Functions ..
 INTRINSIC MAX, MIN
* ..
* .. External Functions ..
 INTEGER ILAENV
 EXTERNAL ILAENV
* ..
* .. Executable Statements ..
*
* Test the input parameters
*
 INFO = 0
 NB = MIN( NBMAX, ILAENV( 1, 'CGEHRD', ' ', N, ILO, IHI, -1 ) )
 LWKOPT = N*NB
 WORK( 1 ) = LWKOPT
 LQUERY = ( LWORK.EQ.-1 )
 IF( N.LT.0 ) THEN
 INFO = -1
 ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN
 INFO = -2
 ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN
 INFO = -3
 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 INFO = -5
 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
 INFO = -8
 END IF
 IF( INFO.NE.0 ) THEN
 CALL XERBLA( 'CGEHRD', -INFO )
 RETURN
 ELSE IF( LQUERY ) THEN
 RETURN
 END IF
*
* Set elements 1:ILO-1 and IHI:N-1 of TAU to zero
*
 DO 10 I = 1, ILO - 1
 TAU( I ) = ZERO
 10 CONTINUE
 DO 20 I = MAX( 1, IHI ), N - 1
 TAU( I ) = ZERO
 20 CONTINUE
*
* Quick return if possible
*
 NH = IHI - ILO + 1
 IF( NH.LE.1 ) THEN
 WORK( 1 ) = 1
 RETURN
 END IF
*
 NBMIN = 2
 IWS = 1
 IF( NB.GT.1 .AND. NB.LT.NH ) THEN
*
* Determine when to cross over from blocked to unblocked code
* (last block is always handled by unblocked code).
*
 NX = MAX( NB, ILAENV( 3, 'CGEHRD', ' ', N, ILO, IHI, -1 ) )
 IF( NX.LT.NH ) THEN
*
* Determine if workspace is large enough for blocked code.
*
 IWS = N*NB
 IF( LWORK.LT.IWS ) THEN
*
* Not enough workspace to use optimal NB: determine the
* minimum value of NB, and reduce NB or force use of
* unblocked code.
*
 NBMIN = MAX( 2, ILAENV( 2, 'CGEHRD', ' ', N, ILO, IHI,
 $ -1 ) )
 IF( LWORK.GE.N*NBMIN ) THEN
 NB = LWORK / N
 ELSE
 NB = 1
 END IF
 END IF
 END IF
 END IF
 LDWORK = N
*
 IF( NB.LT.NBMIN .OR. NB.GE.NH ) THEN
*
* Use unblocked code below
*
 I = ILO
*
 ELSE
*
* Use blocked code
*