Update brief description of function purpose (Reference-LAPACK 831)

lapack-netlib/SRC/cgegs.f

@@ -0,0 +1,528 @@
+*> \brief <b> CGEGS computes the eigenvalues, Schur form, and, optionally, the left and or/right Schur vectors of a complex matrix pair (A,B)</b>
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download CGEGS + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgegs.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgegs.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgegs.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE CGEGS( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHA, BETA,
+*                         VSL, LDVSL, VSR, LDVSR, WORK, LWORK, RWORK,
+*                         INFO )
+*
+*       .. 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( * )
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> This routine is deprecated and has been replaced by routine CGGES.
+*>
+*> CGEGS computes the eigenvalues, Schur form, and, optionally, the
+*> left and or/right Schur vectors of a complex matrix pair (A,B).
+*> Given two square matrices A and B, the generalized Schur
+*> factorization has the form
+*>
+*>    A = Q*S*Z**H,  B = Q*T*Z**H
+*>
+*> where Q and Z are unitary matrices and S and T are upper triangular.
+*> The columns of Q are the left Schur vectors
+*> and the columns of Z are the right Schur vectors.
+*>
+*> If only the eigenvalues of (A,B) are needed, the driver routine
+*> CGEGV should be used instead.  See CGEGV for a description of the
+*> eigenvalues of the generalized nonsymmetric eigenvalue problem
+*> (GNEP).
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] JOBVSL
+*> \verbatim
+*>          JOBVSL is CHARACTER*1
+*>          = 'N':  do not compute the left Schur vectors;
+*>          = 'V':  compute the left Schur vectors (returned in VSL).
+*> \endverbatim
+*>
+*> \param[in] JOBVSR
+*> \verbatim
+*>          JOBVSR is CHARACTER*1
+*>          = 'N':  do not compute the right Schur vectors;
+*>          = 'V':  compute the right Schur vectors (returned in VSR).
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the matrices A, B, VSL, and VSR.  N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*>          A is COMPLEX array, dimension (LDA, N)
+*>          On entry, the matrix A.
+*>          On exit, the upper triangular matrix S from the generalized
+*>          Schur factorization.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>          The leading dimension of A.  LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] B
+*> \verbatim
+*>          B is COMPLEX array, dimension (LDB, N)
+*>          On entry, the matrix B.
+*>          On exit, the upper triangular matrix T from the generalized
+*>          Schur factorization.
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*>          LDB is INTEGER
+*>          The leading dimension of B.  LDB >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] ALPHA
+*> \verbatim
+*>          ALPHA is COMPLEX array, dimension (N)
+*>          The complex scalars alpha that define the eigenvalues of
+*>          GNEP.  ALPHA(j) = S(j,j), the diagonal element of the Schur
+*>          form of A.
+*> \endverbatim
+*>
+*> \param[out] BETA
+*> \verbatim
+*>          BETA is COMPLEX array, dimension (N)
+*>          The non-negative real scalars beta that define the
+*>          eigenvalues of GNEP.  BETA(j) = T(j,j), the diagonal element
+*>          of the triangular factor T.
+*>
+*>          Together, the quantities alpha = ALPHA(j) and beta = BETA(j)
+*>          represent the j-th eigenvalue of the matrix pair (A,B), in
+*>          one of the forms lambda = alpha/beta or mu = beta/alpha.
+*>          Since either lambda or mu may overflow, they should not,
+*>          in general, be computed.
+*> \endverbatim
+*>
+*> \param[out] VSL
+*> \verbatim
+*>          VSL is COMPLEX array, dimension (LDVSL,N)
+*>          If JOBVSL = 'V', the matrix of left Schur vectors Q.
+*>          Not referenced if JOBVSL = 'N'.
+*> \endverbatim
+*>
+*> \param[in] LDVSL
+*> \verbatim
+*>          LDVSL is INTEGER
+*>          The leading dimension of the matrix VSL. LDVSL >= 1, and
+*>          if JOBVSL = 'V', LDVSL >= N.
+*> \endverbatim
+*>
+*> \param[out] VSR
+*> \verbatim
+*>          VSR is COMPLEX array, dimension (LDVSR,N)
+*>          If JOBVSR = 'V', the matrix of right Schur vectors Z.
+*>          Not referenced if JOBVSR = 'N'.
+*> \endverbatim
+*>
+*> \param[in] LDVSR
+*> \verbatim
+*>          LDVSR is INTEGER
+*>          The leading dimension of the matrix VSR. LDVSR >= 1, and
+*>          if JOBVSR = 'V', LDVSR >= N.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*>          LWORK is 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.
+*> \endverbatim
+*>
+*> \param[out] RWORK
+*> \verbatim
+*>          RWORK is REAL array, dimension (3*N)
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*>          INFO is 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)
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \ingroup complexGEeigen
+*
+*  =====================================================================
+      SUBROUTINE CGEGS( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHA, BETA,
+     $                  VSL, LDVSL, VSR, LDVSR, WORK, LWORK, RWORK,
+     $                  INFO )
+*
+*  -- LAPACK driver routine --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. 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( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. 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

lapack-netlib/SRC/cgegv.f

@@ -0,0 +1,703 @@
+*> \brief <b> CGEGV computes the eigenvalues and, optionally, the left and/or right eigenvectors of a complex matrix pair (A,B).</b>
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download CGEGV + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgegv.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgegv.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgegv.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE CGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
+*                         VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
+*
+*       .. 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( * )
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> This routine is deprecated and has been replaced by routine CGGEV.
+*>
+*> CGEGV computes the eigenvalues and, optionally, the left and/or right
+*> eigenvectors of a complex matrix pair (A,B).
+*> Given two square matrices A and B,
+*> the generalized nonsymmetric eigenvalue problem (GNEP) is to find the
+*> eigenvalues lambda and corresponding (non-zero) eigenvectors x such
+*> that
+*>    A*x = lambda*B*x.
+*>
+*> An alternate form is to find the eigenvalues mu and corresponding
+*> eigenvectors y such that
+*>    mu*A*y = B*y.
+*>
+*> These two forms are equivalent with mu = 1/lambda and x = y if
+*> neither lambda nor mu is zero.  In order to deal with the case that
+*> lambda or mu is zero or small, two values alpha and beta are returned
+*> for each eigenvalue, such that lambda = alpha/beta and
+*> mu = beta/alpha.
+*>
+*> The vectors x and y in the above equations are right eigenvectors of
+*> the matrix pair (A,B).  Vectors u and v satisfying
+*>    u**H*A = lambda*u**H*B  or  mu*v**H*A = v**H*B
+*> are left eigenvectors of (A,B).
+*>
+*> Note: this routine performs "full balancing" on A and B
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] JOBVL
+*> \verbatim
+*>          JOBVL is CHARACTER*1
+*>          = 'N':  do not compute the left generalized eigenvectors;
+*>          = 'V':  compute the left generalized eigenvectors (returned
+*>                  in VL).
+*> \endverbatim
+*>
+*> \param[in] JOBVR
+*> \verbatim
+*>          JOBVR is CHARACTER*1
+*>          = 'N':  do not compute the right generalized eigenvectors;
+*>          = 'V':  compute the right generalized eigenvectors (returned
+*>                  in VR).
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the matrices A, B, VL, and VR.  N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*>          A is COMPLEX array, dimension (LDA, N)
+*>          On entry, the matrix A.
+*>          If JOBVL = 'V' or JOBVR = 'V', then on exit A
+*>          contains the Schur form of A from the generalized Schur
+*>          factorization of the pair (A,B) after balancing.  If no
+*>          eigenvectors were computed, then only the diagonal elements
+*>          of the Schur form will be correct.  See CGGHRD and CHGEQZ
+*>          for details.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>          The leading dimension of A.  LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] B
+*> \verbatim
+*>          B is COMPLEX array, dimension (LDB, N)
+*>          On entry, the matrix B.
+*>          If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the
+*>          upper triangular matrix obtained from B in the generalized
+*>          Schur factorization of the pair (A,B) after balancing.
+*>          If no eigenvectors were computed, then only the diagonal
+*>          elements of B will be correct.  See CGGHRD and CHGEQZ for
+*>          details.
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*>          LDB is INTEGER
+*>          The leading dimension of B.  LDB >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] ALPHA
+*> \verbatim
+*>          ALPHA is COMPLEX array, dimension (N)
+*>          The complex scalars alpha that define the eigenvalues of
+*>          GNEP.
+*> \endverbatim
+*>
+*> \param[out] BETA
+*> \verbatim
+*>          BETA is COMPLEX array, dimension (N)
+*>          The complex scalars beta that define the eigenvalues of GNEP.
+*>
+*>          Together, the quantities alpha = ALPHA(j) and beta = BETA(j)
+*>          represent the j-th eigenvalue of the matrix pair (A,B), in
+*>          one of the forms lambda = alpha/beta or mu = beta/alpha.
+*>          Since either lambda or mu may overflow, they should not,
+*>          in general, be computed.
+*> \endverbatim
+*>
+*> \param[out] VL
+*> \verbatim
+*>          VL is COMPLEX array, dimension (LDVL,N)
+*>          If JOBVL = 'V', the left eigenvectors u(j) are stored
+*>          in the columns of VL, in the same order as their eigenvalues.
+*>          Each eigenvector is scaled so that its largest component has
+*>          abs(real part) + abs(imag. part) = 1, except for eigenvectors
+*>          corresponding to an eigenvalue with alpha = beta = 0, which
+*>          are set to zero.
+*>          Not referenced if JOBVL = 'N'.
+*> \endverbatim
+*>
+*> \param[in] LDVL
+*> \verbatim
+*>          LDVL is INTEGER
+*>          The leading dimension of the matrix VL. LDVL >= 1, and
+*>          if JOBVL = 'V', LDVL >= N.
+*> \endverbatim
+*>
+*> \param[out] VR
+*> \verbatim
+*>          VR is COMPLEX array, dimension (LDVR,N)
+*>          If JOBVR = 'V', the right eigenvectors x(j) are stored
+*>          in the columns of VR, in the same order as their eigenvalues.
+*>          Each eigenvector is scaled so that its largest component has
+*>          abs(real part) + abs(imag. part) = 1, except for eigenvectors
+*>          corresponding to an eigenvalue with alpha = beta = 0, which
+*>          are set to zero.
+*>          Not referenced if JOBVR = 'N'.
+*> \endverbatim
+*>
+*> \param[in] LDVR
+*> \verbatim
+*>          LDVR is INTEGER
+*>          The leading dimension of the matrix VR. LDVR >= 1, and
+*>          if JOBVR = 'V', LDVR >= N.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*>          LWORK is 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.
+*> \endverbatim
+*>
+*> \param[out] RWORK
+*> \verbatim
+*>          RWORK is REAL array, dimension (8*N)
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*>          INFO is 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)
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \ingroup complexGEeigen
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>
+*>  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.
+*> \endverbatim
+*>
+*  =====================================================================
+      SUBROUTINE CGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
+     $                  VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
+*
+*  -- LAPACK driver routine --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. 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( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. 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

lapack-netlib/SRC/dgegs.f

@@ -0,0 +1,538 @@
+*> \brief <b> DGEGS computes the eigenvalues, real Schur form, and, optionally, the left and/or right Schur vectors of a real matrix pair (A,B)</b>
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download DGEGS + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgegs.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgegs.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgegs.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE DGEGS( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHAR,
+*                         ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK,
+*                         LWORK, INFO )
+*
+*       .. Scalar Arguments ..
+*       CHARACTER          JOBVSL, JOBVSR
+*       INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
+*       ..
+*       .. Array Arguments ..
+*       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
+*      $                   B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
+*      $                   VSR( LDVSR, * ), WORK( * )
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> This routine is deprecated and has been replaced by routine DGGES.
+*>
+*> DGEGS computes the eigenvalues, real Schur form, and, optionally,
+*> left and or/right Schur vectors of a real matrix pair (A,B).
+*> Given two square matrices A and B, the generalized real Schur
+*> factorization has the form
+*>
+*>   A = Q*S*Z**T,  B = Q*T*Z**T
+*>
+*> where Q and Z are orthogonal matrices, T is upper triangular, and S
+*> is an upper quasi-triangular matrix with 1-by-1 and 2-by-2 diagonal
+*> blocks, the 2-by-2 blocks corresponding to complex conjugate pairs
+*> of eigenvalues of (A,B).  The columns of Q are the left Schur vectors
+*> and the columns of Z are the right Schur vectors.
+*>
+*> If only the eigenvalues of (A,B) are needed, the driver routine
+*> DGEGV should be used instead.  See DGEGV for a description of the
+*> eigenvalues of the generalized nonsymmetric eigenvalue problem
+*> (GNEP).
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] JOBVSL
+*> \verbatim
+*>          JOBVSL is CHARACTER*1
+*>          = 'N':  do not compute the left Schur vectors;
+*>          = 'V':  compute the left Schur vectors (returned in VSL).
+*> \endverbatim
+*>
+*> \param[in] JOBVSR
+*> \verbatim
+*>          JOBVSR is CHARACTER*1
+*>          = 'N':  do not compute the right Schur vectors;
+*>          = 'V':  compute the right Schur vectors (returned in VSR).
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the matrices A, B, VSL, and VSR.  N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*>          A is DOUBLE PRECISION array, dimension (LDA, N)
+*>          On entry, the matrix A.
+*>          On exit, the upper quasi-triangular matrix S from the
+*>          generalized real Schur factorization.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>          The leading dimension of A.  LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] B
+*> \verbatim
+*>          B is DOUBLE PRECISION array, dimension (LDB, N)
+*>          On entry, the matrix B.
+*>          On exit, the upper triangular matrix T from the generalized
+*>          real Schur factorization.
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*>          LDB is INTEGER
+*>          The leading dimension of B.  LDB >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] ALPHAR
+*> \verbatim
+*>          ALPHAR is DOUBLE PRECISION array, dimension (N)
+*>          The real parts of each scalar alpha defining an eigenvalue
+*>          of GNEP.
+*> \endverbatim
+*>
+*> \param[out] ALPHAI
+*> \verbatim
+*>          ALPHAI is DOUBLE PRECISION array, dimension (N)
+*>          The imaginary parts of each scalar alpha defining an
+*>          eigenvalue of GNEP.  If ALPHAI(j) is zero, then the j-th
+*>          eigenvalue is real; if positive, then the j-th and (j+1)-st
+*>          eigenvalues are a complex conjugate pair, with
+*>          ALPHAI(j+1) = -ALPHAI(j).
+*> \endverbatim
+*>
+*> \param[out] BETA
+*> \verbatim
+*>          BETA is DOUBLE PRECISION array, dimension (N)
+*>          The scalars beta that define the eigenvalues of GNEP.
+*>          Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
+*>          beta = BETA(j) represent the j-th eigenvalue of the matrix
+*>          pair (A,B), in one of the forms lambda = alpha/beta or
+*>          mu = beta/alpha.  Since either lambda or mu may overflow,
+*>          they should not, in general, be computed.
+*> \endverbatim
+*>
+*> \param[out] VSL
+*> \verbatim
+*>          VSL is DOUBLE PRECISION array, dimension (LDVSL,N)
+*>          If JOBVSL = 'V', the matrix of left Schur vectors Q.
+*>          Not referenced if JOBVSL = 'N'.
+*> \endverbatim
+*>
+*> \param[in] LDVSL
+*> \verbatim
+*>          LDVSL is INTEGER
+*>          The leading dimension of the matrix VSL. LDVSL >=1, and
+*>          if JOBVSL = 'V', LDVSL >= N.
+*> \endverbatim
+*>
+*> \param[out] VSR
+*> \verbatim
+*>          VSR is DOUBLE PRECISION array, dimension (LDVSR,N)
+*>          If JOBVSR = 'V', the matrix of right Schur vectors Z.
+*>          Not referenced if JOBVSR = 'N'.
+*> \endverbatim
+*>
+*> \param[in] LDVSR
+*> \verbatim
+*>          LDVSR is INTEGER
+*>          The leading dimension of the matrix VSR. LDVSR >= 1, and
+*>          if JOBVSR = 'V', LDVSR >= N.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*>          LWORK is INTEGER
+*>          The dimension of the array WORK.  LWORK >= max(1,4*N).
+*>          For good performance, LWORK must generally be larger.
+*>          To compute the optimal value of LWORK, call ILAENV to get
+*>          blocksizes (for DGEQRF, DORMQR, and DORGQR.)  Then compute:
+*>          NB  -- MAX of the blocksizes for DGEQRF, DORMQR, and DORGQR
+*>          The optimal LWORK is  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.
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*>          INFO is 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 ALPHAR(j), ALPHAI(j), and BETA(j) should
+*>                be correct for j=INFO+1,...,N.
+*>          > N:  errors that usually indicate LAPACK problems:
+*>                =N+1: error return from DGGBAL
+*>                =N+2: error return from DGEQRF
+*>                =N+3: error return from DORMQR
+*>                =N+4: error return from DORGQR
+*>                =N+5: error return from DGGHRD
+*>                =N+6: error return from DHGEQZ (other than failed
+*>                                                iteration)
+*>                =N+7: error return from DGGBAK (computing VSL)
+*>                =N+8: error return from DGGBAK (computing VSR)
+*>                =N+9: error return from DLASCL (various places)
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \ingroup doubleGEeigen
+*
+*  =====================================================================
+      SUBROUTINE DGEGS( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHAR,
+     $                  ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK,
+     $                  LWORK, INFO )
+*
+*  -- LAPACK driver routine --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      CHARACTER          JOBVSL, JOBVSR
+      INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
+     $                   B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
+     $                   VSR( LDVSR, * ), WORK( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            ILASCL, ILBSCL, ILVSL, ILVSR, LQUERY
+      INTEGER            ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
+     $                   IRIGHT, IROWS, ITAU, IWORK, LOPT, LWKMIN,
+     $                   LWKOPT, NB, NB1, NB2, NB3
+      DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
+     $                   SAFMIN, SMLNUM
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLACPY,
+     $                   DLASCL, DLASET, DORGQR, DORMQR, XERBLA
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      INTEGER            ILAENV
+      DOUBLE PRECISION   DLAMCH, DLANGE
+      EXTERNAL           LSAME, ILAENV, DLAMCH, DLANGE
+*     ..
+*     .. 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( 4*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 = -12
+      ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
+         INFO = -14
+      ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
+         INFO = -16
+      END IF
+*
+      IF( INFO.EQ.0 ) THEN
+         NB1 = ILAENV( 1, 'DGEQRF', ' ', N, N, -1, -1 )
+         NB2 = ILAENV( 1, 'DORMQR', ' ', N, N, N, -1 )
+         NB3 = ILAENV( 1, 'DORGQR', ' ', N, N, N, -1 )
+         NB = MAX( NB1, NB2, NB3 )
+         LOPT = 2*N + N*( NB+1 )
+         WORK( 1 ) = LOPT
+      END IF
+*
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'DGEGS ', -INFO )
+         RETURN
+      ELSE IF( LQUERY ) THEN
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
+*     Get machine constants
+*
+      EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
+      SAFMIN = DLAMCH( 'S' )
+      SMLNUM = N*SAFMIN / EPS
+      BIGNUM = ONE / SMLNUM
+*
+*     Scale A if max element outside range [SMLNUM,BIGNUM]
+*
+      ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
+      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 DLASCL( '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 = DLANGE( 'M', N, N, B, LDB, WORK )
+      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 DLASCL( '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
+*     Workspace layout:  (2*N words -- "work..." not actually used)
+*        left_permutation, right_permutation, work...
+*
+      ILEFT = 1
+      IRIGHT = N + 1
+      IWORK = IRIGHT + N
+      CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
+     $             WORK( IRIGHT ), WORK( IWORK ), 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
+*     Workspace layout:  ("work..." must have at least N words)
+*        left_permutation, right_permutation, tau, work...
+*
+      IROWS = IHI + 1 - ILO
+      ICOLS = N + 1 - ILO
+      ITAU = IWORK
+      IWORK = ITAU + IROWS
+      CALL DGEQRF( 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 DORMQR( 'L', 'T', 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 DLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL )
+         CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
+     $                VSL( ILO+1, ILO ), LDVSL )
+         CALL DORGQR( 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 DLASET( 'Full', N, N, ZERO, ONE, VSR, LDVSR )
+*
+*     Reduce to generalized Hessenberg form
+*
+      CALL DGGHRD( 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
+*     Workspace layout:  ("work..." must have at least 1 word)
+*        left_permutation, right_permutation, work...
+*
+      IWORK = ITAU
+      CALL DHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
+     $             ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
+     $             WORK( IWORK ), LWORK+1-IWORK, 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 DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
+     $                WORK( IRIGHT ), N, VSL, LDVSL, IINFO )
+         IF( IINFO.NE.0 ) THEN
+            INFO = N + 7
+            GO TO 10
+         END IF
+      END IF
+      IF( ILVSR ) THEN
+         CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
+     $                WORK( 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 DLASCL( 'H', -1, -1, ANRMTO, ANRM, N, N, A, LDA, IINFO )
+         IF( IINFO.NE.0 ) THEN
+            INFO = N + 9
+            RETURN
+         END IF
+         CALL DLASCL( 'G', -1, -1, ANRMTO, ANRM, N, 1, ALPHAR, N,
+     $                IINFO )
+         IF( IINFO.NE.0 ) THEN
+            INFO = N + 9
+            RETURN
+         END IF
+         CALL DLASCL( 'G', -1, -1, ANRMTO, ANRM, N, 1, ALPHAI, N,
+     $                IINFO )
+         IF( IINFO.NE.0 ) THEN
+            INFO = N + 9
+            RETURN
+         END IF
+      END IF
+*
+      IF( ILBSCL ) THEN
+         CALL DLASCL( 'U', -1, -1, BNRMTO, BNRM, N, N, B, LDB, IINFO )
+         IF( IINFO.NE.0 ) THEN
+            INFO = N + 9
+            RETURN
+         END IF
+         CALL DLASCL( '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 DGEGS
+*
+      END

lapack-netlib/SRC/dgegv.f

@@ -0,0 +1,766 @@
+*> \brief <b> DGEGV computes the eigenvalues and, optionally, the left and/or right eigenvectors of a real matrix pair (A,B).</b>
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download DGEGV + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgegv.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgegv.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgegv.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE DGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
+*                         BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
+*
+*       .. Scalar Arguments ..
+*       CHARACTER          JOBVL, JOBVR
+*       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
+*       ..
+*       .. Array Arguments ..
+*       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
+*      $                   B( LDB, * ), BETA( * ), VL( LDVL, * ),
+*      $                   VR( LDVR, * ), WORK( * )
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> This routine is deprecated and has been replaced by routine DGGEV.
+*>
+*> DGEGV computes the eigenvalues and, optionally, the left and/or right
+*> eigenvectors of a real matrix pair (A,B).
+*> Given two square matrices A and B,
+*> the generalized nonsymmetric eigenvalue problem (GNEP) is to find the
+*> eigenvalues lambda and corresponding (non-zero) eigenvectors x such
+*> that
+*>
+*>    A*x = lambda*B*x.
+*>
+*> An alternate form is to find the eigenvalues mu and corresponding
+*> eigenvectors y such that
+*>
+*>    mu*A*y = B*y.
+*>
+*> These two forms are equivalent with mu = 1/lambda and x = y if
+*> neither lambda nor mu is zero.  In order to deal with the case that
+*> lambda or mu is zero or small, two values alpha and beta are returned
+*> for each eigenvalue, such that lambda = alpha/beta and
+*> mu = beta/alpha.
+*>
+*> The vectors x and y in the above equations are right eigenvectors of
+*> the matrix pair (A,B).  Vectors u and v satisfying
+*>
+*>    u**H*A = lambda*u**H*B  or  mu*v**H*A = v**H*B
+*>
+*> are left eigenvectors of (A,B).
+*>
+*> Note: this routine performs "full balancing" on A and B
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] JOBVL
+*> \verbatim
+*>          JOBVL is CHARACTER*1
+*>          = 'N':  do not compute the left generalized eigenvectors;
+*>          = 'V':  compute the left generalized eigenvectors (returned
+*>                  in VL).
+*> \endverbatim
+*>
+*> \param[in] JOBVR
+*> \verbatim
+*>          JOBVR is CHARACTER*1
+*>          = 'N':  do not compute the right generalized eigenvectors;
+*>          = 'V':  compute the right generalized eigenvectors (returned
+*>                  in VR).
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the matrices A, B, VL, and VR.  N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*>          A is DOUBLE PRECISION array, dimension (LDA, N)
+*>          On entry, the matrix A.
+*>          If JOBVL = 'V' or JOBVR = 'V', then on exit A
+*>          contains the real Schur form of A from the generalized Schur
+*>          factorization of the pair (A,B) after balancing.
+*>          If no eigenvectors were computed, then only the diagonal
+*>          blocks from the Schur form will be correct.  See DGGHRD and
+*>          DHGEQZ for details.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>          The leading dimension of A.  LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] B
+*> \verbatim
+*>          B is DOUBLE PRECISION array, dimension (LDB, N)
+*>          On entry, the matrix B.
+*>          If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the
+*>          upper triangular matrix obtained from B in the generalized
+*>          Schur factorization of the pair (A,B) after balancing.
+*>          If no eigenvectors were computed, then only those elements of
+*>          B corresponding to the diagonal blocks from the Schur form of
+*>          A will be correct.  See DGGHRD and DHGEQZ for details.
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*>          LDB is INTEGER
+*>          The leading dimension of B.  LDB >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] ALPHAR
+*> \verbatim
+*>          ALPHAR is DOUBLE PRECISION array, dimension (N)
+*>          The real parts of each scalar alpha defining an eigenvalue of
+*>          GNEP.
+*> \endverbatim
+*>
+*> \param[out] ALPHAI
+*> \verbatim
+*>          ALPHAI is DOUBLE PRECISION array, dimension (N)
+*>          The imaginary parts of each scalar alpha defining an
+*>          eigenvalue of GNEP.  If ALPHAI(j) is zero, then the j-th
+*>          eigenvalue is real; if positive, then the j-th and
+*>          (j+1)-st eigenvalues are a complex conjugate pair, with
+*>          ALPHAI(j+1) = -ALPHAI(j).
+*> \endverbatim
+*>
+*> \param[out] BETA
+*> \verbatim
+*>          BETA is DOUBLE PRECISION array, dimension (N)
+*>          The scalars beta that define the eigenvalues of GNEP.
+*>
+*>          Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
+*>          beta = BETA(j) represent the j-th eigenvalue of the matrix
+*>          pair (A,B), in one of the forms lambda = alpha/beta or
+*>          mu = beta/alpha.  Since either lambda or mu may overflow,
+*>          they should not, in general, be computed.
+*> \endverbatim
+*>
+*> \param[out] VL
+*> \verbatim
+*>          VL is DOUBLE PRECISION array, dimension (LDVL,N)
+*>          If JOBVL = 'V', the left eigenvectors u(j) are stored
+*>          in the columns of VL, in the same order as their eigenvalues.
+*>          If the j-th eigenvalue is real, then u(j) = VL(:,j).
+*>          If the j-th and (j+1)-st eigenvalues form a complex conjugate
+*>          pair, then
+*>             u(j) = VL(:,j) + i*VL(:,j+1)
+*>          and
+*>            u(j+1) = VL(:,j) - i*VL(:,j+1).
+*>
+*>          Each eigenvector is scaled so that its largest component has
+*>          abs(real part) + abs(imag. part) = 1, except for eigenvectors
+*>          corresponding to an eigenvalue with alpha = beta = 0, which
+*>          are set to zero.
+*>          Not referenced if JOBVL = 'N'.
+*> \endverbatim
+*>
+*> \param[in] LDVL
+*> \verbatim
+*>          LDVL is INTEGER
+*>          The leading dimension of the matrix VL. LDVL >= 1, and
+*>          if JOBVL = 'V', LDVL >= N.
+*> \endverbatim
+*>
+*> \param[out] VR
+*> \verbatim
+*>          VR is DOUBLE PRECISION array, dimension (LDVR,N)
+*>          If JOBVR = 'V', the right eigenvectors x(j) are stored
+*>          in the columns of VR, in the same order as their eigenvalues.
+*>          If the j-th eigenvalue is real, then x(j) = VR(:,j).
+*>          If the j-th and (j+1)-st eigenvalues form a complex conjugate
+*>          pair, then
+*>            x(j) = VR(:,j) + i*VR(:,j+1)
+*>          and
+*>            x(j+1) = VR(:,j) - i*VR(:,j+1).
+*>
+*>          Each eigenvector is scaled so that its largest component has
+*>          abs(real part) + abs(imag. part) = 1, except for eigenvalues
+*>          corresponding to an eigenvalue with alpha = beta = 0, which
+*>          are set to zero.
+*>          Not referenced if JOBVR = 'N'.
+*> \endverbatim
+*>
+*> \param[in] LDVR
+*> \verbatim
+*>          LDVR is INTEGER
+*>          The leading dimension of the matrix VR. LDVR >= 1, and
+*>          if JOBVR = 'V', LDVR >= N.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*>          LWORK is INTEGER
+*>          The dimension of the array WORK.  LWORK >= max(1,8*N).
+*>          For good performance, LWORK must generally be larger.
+*>          To compute the optimal value of LWORK, call ILAENV to get
+*>          blocksizes (for DGEQRF, DORMQR, and DORGQR.)  Then compute:
+*>          NB  -- MAX of the blocksizes for DGEQRF, DORMQR, and DORGQR;
+*>          The optimal LWORK is:
+*>              2*N + MAX( 6*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.
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*>          INFO is 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 ALPHAR(j), ALPHAI(j), and BETA(j)
+*>                should be correct for j=INFO+1,...,N.
+*>          > N:  errors that usually indicate LAPACK problems:
+*>                =N+1: error return from DGGBAL
+*>                =N+2: error return from DGEQRF
+*>                =N+3: error return from DORMQR
+*>                =N+4: error return from DORGQR
+*>                =N+5: error return from DGGHRD
+*>                =N+6: error return from DHGEQZ (other than failed
+*>                                                iteration)
+*>                =N+7: error return from DTGEVC
+*>                =N+8: error return from DGGBAK (computing VL)
+*>                =N+9: error return from DGGBAK (computing VR)
+*>                =N+10: error return from DLASCL (various calls)
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \ingroup doubleGEeigen
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>
+*>  Balancing
+*>  ---------
+*>
+*>  This driver calls DGGBAL 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, DGGBAK 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 real Schur
+*>  form[*] of the "balanced" versions of A and B.  If no eigenvectors
+*>  are computed, then only the diagonal blocks will be correct.
+*>
+*>  [*] See DHGEQZ, DGEGS, or read the book "Matrix Computations",
+*>      by Golub & van Loan, pub. by Johns Hopkins U. Press.
+*> \endverbatim
+*>
+*  =====================================================================
+      SUBROUTINE DGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
+     $                  BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
+*
+*  -- LAPACK driver routine --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      CHARACTER          JOBVL, JOBVR
+      INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
+     $                   B( LDB, * ), BETA( * ), VL( LDVL, * ),
+     $                   VR( LDVR, * ), WORK( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            ILIMIT, ILV, ILVL, ILVR, LQUERY
+      CHARACTER          CHTEMP
+      INTEGER            ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
+     $                   IN, IRIGHT, IROWS, ITAU, IWORK, JC, JR, LOPT,
+     $                   LWKMIN, LWKOPT, NB, NB1, NB2, NB3
+      DOUBLE PRECISION   ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM,
+     $                   BNRM1, BNRM2, EPS, ONEPLS, SAFMAX, SAFMIN,
+     $                   SALFAI, SALFAR, SBETA, SCALE, TEMP
+*     ..
+*     .. Local Arrays ..
+      LOGICAL            LDUMMA( 1 )
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLACPY,
+     $                   DLASCL, DLASET, DORGQR, DORMQR, DTGEVC, XERBLA
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      INTEGER            ILAENV
+      DOUBLE PRECISION   DLAMCH, DLANGE
+      EXTERNAL           LSAME, ILAENV, DLAMCH, DLANGE
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, INT, MAX
+*     ..
+*     .. 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( 8*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 = -12
+      ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
+         INFO = -14
+      ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
+         INFO = -16
+      END IF
+*
+      IF( INFO.EQ.0 ) THEN
+         NB1 = ILAENV( 1, 'DGEQRF', ' ', N, N, -1, -1 )
+         NB2 = ILAENV( 1, 'DORMQR', ' ', N, N, N, -1 )
+         NB3 = ILAENV( 1, 'DORGQR', ' ', N, N, N, -1 )
+         NB = MAX( NB1, NB2, NB3 )
+         LOPT = 2*N + MAX( 6*N, N*( NB+1 ) )
+         WORK( 1 ) = LOPT
+      END IF
+*
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'DGEGV ', -INFO )
+         RETURN
+      ELSE IF( LQUERY ) THEN
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
+*     Get machine constants
+*
+      EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
+      SAFMIN = DLAMCH( 'S' )
+      SAFMIN = SAFMIN + SAFMIN
+      SAFMAX = ONE / SAFMIN
+      ONEPLS = ONE + ( 4*EPS )
+*
+*     Scale A
+*
+      ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
+      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 DLASCL( '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 = DLANGE( 'M', N, N, B, LDB, WORK )
+      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 DLASCL( '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
+*     Workspace layout:  (8*N words -- "work" requires 6*N words)
+*        left_permutation, right_permutation, work...
+*
+      ILEFT = 1
+      IRIGHT = N + 1
+      IWORK = IRIGHT + N
+      CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
+     $             WORK( IRIGHT ), WORK( IWORK ), IINFO )
+      IF( IINFO.NE.0 ) THEN
+         INFO = N + 1
+         GO TO 120
+      END IF
+*
+*     Reduce B to triangular form, and initialize VL and/or VR
+*     Workspace layout:  ("work..." must have at least N words)
+*        left_permutation, right_permutation, tau, work...
+*
+      IROWS = IHI + 1 - ILO
+      IF( ILV ) THEN
+         ICOLS = N + 1 - ILO
+      ELSE
+         ICOLS = IROWS
+      END IF
+      ITAU = IWORK
+      IWORK = ITAU + IROWS
+      CALL DGEQRF( 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 120
+      END IF
+*
+      CALL DORMQR( 'L', 'T', 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 120
+      END IF
+*
+      IF( ILVL ) THEN
+         CALL DLASET( 'Full', N, N, ZERO, ONE, VL, LDVL )
+         CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
+     $                VL( ILO+1, ILO ), LDVL )
+         CALL DORGQR( 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 120
+         END IF
+      END IF
+*
+      IF( ILVR )
+     $   CALL DLASET( 'Full', N, N, ZERO, ONE, VR, LDVR )
+*
+*     Reduce to generalized Hessenberg form
+*
+      IF( ILV ) THEN
+*
+*        Eigenvectors requested -- work on whole matrix.
+*
+         CALL DGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
+     $                LDVL, VR, LDVR, IINFO )
+      ELSE
+         CALL DGGHRD( '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 120
+      END IF
+*
+*     Perform QZ algorithm
+*     Workspace layout:  ("work..." must have at least 1 word)
+*        left_permutation, right_permutation, work...
+*
+      IWORK = ITAU
+      IF( ILV ) THEN
+         CHTEMP = 'S'
+      ELSE
+         CHTEMP = 'E'
+      END IF
+      CALL DHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
+     $             ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
+     $             WORK( IWORK ), LWORK+1-IWORK, 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 120
+      END IF
+*
+      IF( ILV ) THEN
+*
+*        Compute Eigenvectors  (DTGEVC requires 6*N words of workspace)
+*
+         IF( ILVL ) THEN
+            IF( ILVR ) THEN
+               CHTEMP = 'B'
+            ELSE
+               CHTEMP = 'L'
+            END IF
+         ELSE
+            CHTEMP = 'R'
+         END IF
+*
+         CALL DTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
+     $                VR, LDVR, N, IN, WORK( IWORK ), IINFO )
+         IF( IINFO.NE.0 ) THEN
+            INFO = N + 7
+            GO TO 120
+         END IF
+*
+*        Undo balancing on VL and VR, rescale
+*
+         IF( ILVL ) THEN
+            CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
+     $                   WORK( IRIGHT ), N, VL, LDVL, IINFO )
+            IF( IINFO.NE.0 ) THEN
+               INFO = N + 8
+               GO TO 120
+            END IF
+            DO 50 JC = 1, N
+               IF( ALPHAI( JC ).LT.ZERO )
+     $            GO TO 50
+               TEMP = ZERO
+               IF( ALPHAI( JC ).EQ.ZERO ) THEN
+                  DO 10 JR = 1, N
+                     TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
+   10             CONTINUE
+               ELSE
+                  DO 20 JR = 1, N
+                     TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
+     $                      ABS( VL( JR, JC+1 ) ) )
+   20             CONTINUE
+               END IF
+               IF( TEMP.LT.SAFMIN )
+     $            GO TO 50
+               TEMP = ONE / TEMP
+               IF( ALPHAI( JC ).EQ.ZERO ) THEN
+                  DO 30 JR = 1, N
+                     VL( JR, JC ) = VL( JR, JC )*TEMP
+   30             CONTINUE
+               ELSE
+                  DO 40 JR = 1, N
+                     VL( JR, JC ) = VL( JR, JC )*TEMP
+                     VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
+   40             CONTINUE
+               END IF
+   50       CONTINUE
+         END IF
+         IF( ILVR ) THEN
+            CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
+     $                   WORK( IRIGHT ), N, VR, LDVR, IINFO )
+            IF( IINFO.NE.0 ) THEN
+               INFO = N + 9
+               GO TO 120
+            END IF
+            DO 100 JC = 1, N
+               IF( ALPHAI( JC ).LT.ZERO )
+     $            GO TO 100
+               TEMP = ZERO
+               IF( ALPHAI( JC ).EQ.ZERO ) THEN
+                  DO 60 JR = 1, N
+                     TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
+   60             CONTINUE
+               ELSE
+                  DO 70 JR = 1, N
+                     TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
+     $                      ABS( VR( JR, JC+1 ) ) )
+   70             CONTINUE
+               END IF
+               IF( TEMP.LT.SAFMIN )
+     $            GO TO 100
+               TEMP = ONE / TEMP
+               IF( ALPHAI( JC ).EQ.ZERO ) THEN
+                  DO 80 JR = 1, N
+                     VR( JR, JC ) = VR( JR, JC )*TEMP
+   80             CONTINUE
+               ELSE
+                  DO 90 JR = 1, N
+                     VR( JR, JC ) = VR( JR, JC )*TEMP
+                     VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
+   90             CONTINUE
+               END IF
+  100       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 110 JC = 1, N
+         ABSAR = ABS( ALPHAR( JC ) )
+         ABSAI = ABS( ALPHAI( JC ) )
+         ABSB = ABS( BETA( JC ) )
+         SALFAR = ANRM*ALPHAR( JC )
+         SALFAI = ANRM*ALPHAI( JC )
+         SBETA = BNRM*BETA( JC )
+         ILIMIT = .FALSE.
+         SCALE = ONE
+*
+*        Check for significant underflow in ALPHAI
+*
+         IF( ABS( SALFAI ).LT.SAFMIN .AND. ABSAI.GE.
+     $       MAX( SAFMIN, EPS*ABSAR, EPS*ABSB ) ) THEN
+            ILIMIT = .TRUE.
+            SCALE = ( ONEPLS*SAFMIN / ANRM1 ) /
+     $              MAX( ONEPLS*SAFMIN, ANRM2*ABSAI )
+*
+         ELSE IF( SALFAI.EQ.ZERO ) THEN
+*
+*           If insignificant underflow in ALPHAI, then make the
+*           conjugate eigenvalue real.
+*
+            IF( ALPHAI( JC ).LT.ZERO .AND. JC.GT.1 ) THEN
+               ALPHAI( JC-1 ) = ZERO
+            ELSE IF( ALPHAI( JC ).GT.ZERO .AND. JC.LT.N ) THEN
+               ALPHAI( JC+1 ) = ZERO
+            END IF
+         END IF
+*
+*        Check for significant underflow in ALPHAR
+*
+         IF( ABS( SALFAR ).LT.SAFMIN .AND. ABSAR.GE.
+     $       MAX( SAFMIN, EPS*ABSAI, EPS*ABSB ) ) THEN
+            ILIMIT = .TRUE.
+            SCALE = MAX( SCALE, ( ONEPLS*SAFMIN / ANRM1 ) /
+     $              MAX( ONEPLS*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, ( ONEPLS*SAFMIN / BNRM1 ) /
+     $              MAX( ONEPLS*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 ALPHAR, ALPHAI, BETA if necessary.
+*
+         IF( ILIMIT ) THEN
+            SALFAR = ( SCALE*ALPHAR( JC ) )*ANRM
+            SALFAI = ( SCALE*ALPHAI( JC ) )*ANRM
+            SBETA = ( SCALE*BETA( JC ) )*BNRM
+         END IF
+         ALPHAR( JC ) = SALFAR
+         ALPHAI( JC ) = SALFAI
+         BETA( JC ) = SBETA
+  110 CONTINUE
+*
+  120 CONTINUE
+      WORK( 1 ) = LWKOPT
+*
+      RETURN
+*
+*     End of DGEGV
+*
+      END

lapack-netlib/SRC/sgegs.f

@@ -0,0 +1,538 @@
+*> \brief <b> SGEGS computes the eigenvalues, real Schur form, and, optionally, the left and/or right Schur vectors of a real matrix pair (A,B)</b>
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download SGEGS + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgegs.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgegs.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgegs.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE SGEGS( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHAR,
+*                         ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK,
+*                         LWORK, INFO )
+*
+*       .. Scalar Arguments ..
+*       CHARACTER          JOBVSL, JOBVSR
+*       INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
+*       ..
+*       .. Array Arguments ..
+*       REAL               A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
+*      $                   B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
+*      $                   VSR( LDVSR, * ), WORK( * )
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> This routine is deprecated and has been replaced by routine SGGES.
+*>
+*> SGEGS computes the eigenvalues, real Schur form, and, optionally,
+*> left and or/right Schur vectors of a real matrix pair (A,B).
+*> Given two square matrices A and B, the generalized real Schur
+*> factorization has the form
+*>
+*>   A = Q*S*Z**T,  B = Q*T*Z**T
+*>
+*> where Q and Z are orthogonal matrices, T is upper triangular, and S
+*> is an upper quasi-triangular matrix with 1-by-1 and 2-by-2 diagonal
+*> blocks, the 2-by-2 blocks corresponding to complex conjugate pairs
+*> of eigenvalues of (A,B).  The columns of Q are the left Schur vectors
+*> and the columns of Z are the right Schur vectors.
+*>
+*> If only the eigenvalues of (A,B) are needed, the driver routine
+*> SGEGV should be used instead.  See SGEGV for a description of the
+*> eigenvalues of the generalized nonsymmetric eigenvalue problem
+*> (GNEP).
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] JOBVSL
+*> \verbatim
+*>          JOBVSL is CHARACTER*1
+*>          = 'N':  do not compute the left Schur vectors;
+*>          = 'V':  compute the left Schur vectors (returned in VSL).
+*> \endverbatim
+*>
+*> \param[in] JOBVSR
+*> \verbatim
+*>          JOBVSR is CHARACTER*1
+*>          = 'N':  do not compute the right Schur vectors;
+*>          = 'V':  compute the right Schur vectors (returned in VSR).
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the matrices A, B, VSL, and VSR.  N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*>          A is REAL array, dimension (LDA, N)
+*>          On entry, the matrix A.
+*>          On exit, the upper quasi-triangular matrix S from the
+*>          generalized real Schur factorization.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>          The leading dimension of A.  LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] B
+*> \verbatim
+*>          B is REAL array, dimension (LDB, N)
+*>          On entry, the matrix B.
+*>          On exit, the upper triangular matrix T from the generalized
+*>          real Schur factorization.
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*>          LDB is INTEGER
+*>          The leading dimension of B.  LDB >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] ALPHAR
+*> \verbatim
+*>          ALPHAR is REAL array, dimension (N)
+*>          The real parts of each scalar alpha defining an eigenvalue
+*>          of GNEP.
+*> \endverbatim
+*>
+*> \param[out] ALPHAI
+*> \verbatim
+*>          ALPHAI is REAL array, dimension (N)
+*>          The imaginary parts of each scalar alpha defining an
+*>          eigenvalue of GNEP.  If ALPHAI(j) is zero, then the j-th
+*>          eigenvalue is real; if positive, then the j-th and (j+1)-st
+*>          eigenvalues are a complex conjugate pair, with
+*>          ALPHAI(j+1) = -ALPHAI(j).
+*> \endverbatim
+*>
+*> \param[out] BETA
+*> \verbatim
+*>          BETA is REAL array, dimension (N)
+*>          The scalars beta that define the eigenvalues of GNEP.
+*>          Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
+*>          beta = BETA(j) represent the j-th eigenvalue of the matrix
+*>          pair (A,B), in one of the forms lambda = alpha/beta or
+*>          mu = beta/alpha.  Since either lambda or mu may overflow,
+*>          they should not, in general, be computed.
+*> \endverbatim
+*>
+*> \param[out] VSL
+*> \verbatim
+*>          VSL is REAL array, dimension (LDVSL,N)
+*>          If JOBVSL = 'V', the matrix of left Schur vectors Q.
+*>          Not referenced if JOBVSL = 'N'.
+*> \endverbatim
+*>
+*> \param[in] LDVSL
+*> \verbatim
+*>          LDVSL is INTEGER
+*>          The leading dimension of the matrix VSL. LDVSL >=1, and
+*>          if JOBVSL = 'V', LDVSL >= N.
+*> \endverbatim
+*>
+*> \param[out] VSR
+*> \verbatim
+*>          VSR is REAL array, dimension (LDVSR,N)
+*>          If JOBVSR = 'V', the matrix of right Schur vectors Z.
+*>          Not referenced if JOBVSR = 'N'.
+*> \endverbatim
+*>
+*> \param[in] LDVSR
+*> \verbatim
+*>          LDVSR is INTEGER
+*>          The leading dimension of the matrix VSR. LDVSR >= 1, and
+*>          if JOBVSR = 'V', LDVSR >= N.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is REAL array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*>          LWORK is INTEGER
+*>          The dimension of the array WORK.  LWORK >= max(1,4*N).
+*>          For good performance, LWORK must generally be larger.
+*>          To compute the optimal value of LWORK, call ILAENV to get
+*>          blocksizes (for SGEQRF, SORMQR, and SORGQR.)  Then compute:
+*>          NB  -- MAX of the blocksizes for SGEQRF, SORMQR, and SORGQR
+*>          The optimal LWORK is  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.
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*>          INFO is 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 ALPHAR(j), ALPHAI(j), and BETA(j) should
+*>                be correct for j=INFO+1,...,N.
+*>          > N:  errors that usually indicate LAPACK problems:
+*>                =N+1: error return from SGGBAL
+*>                =N+2: error return from SGEQRF
+*>                =N+3: error return from SORMQR
+*>                =N+4: error return from SORGQR
+*>                =N+5: error return from SGGHRD
+*>                =N+6: error return from SHGEQZ (other than failed
+*>                                                iteration)
+*>                =N+7: error return from SGGBAK (computing VSL)
+*>                =N+8: error return from SGGBAK (computing VSR)
+*>                =N+9: error return from SLASCL (various places)
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \ingroup realGEeigen
+*
+*  =====================================================================
+      SUBROUTINE SGEGS( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHAR,
+     $                  ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK,
+     $                  LWORK, INFO )
+*
+*  -- LAPACK driver routine --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      CHARACTER          JOBVSL, JOBVSR
+      INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
+*     ..
+*     .. Array Arguments ..
+      REAL               A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
+     $                   B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
+     $                   VSR( LDVSR, * ), WORK( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ZERO, ONE
+      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            ILASCL, ILBSCL, ILVSL, ILVSR, LQUERY
+      INTEGER            ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT,
+     $                   ILO, IRIGHT, IROWS, ITAU, IWORK, LOPT, LWKMIN,
+     $                   LWKOPT, NB, NB1, NB2, NB3
+      REAL               ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
+     $                   SAFMIN, SMLNUM
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           SGEQRF, SGGBAK, SGGBAL, SGGHRD, SHGEQZ, SLACPY,
+     $                   SLASCL, SLASET, SORGQR, SORMQR, XERBLA
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      INTEGER            ILAENV
+      REAL               SLAMCH, SLANGE
+      EXTERNAL           ILAENV, LSAME, SLAMCH, SLANGE
+*     ..
+*     .. 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( 4*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 = -12
+      ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
+         INFO = -14
+      ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
+         INFO = -16
+      END IF
+*
+      IF( INFO.EQ.0 ) THEN
+         NB1 = ILAENV( 1, 'SGEQRF', ' ', N, N, -1, -1 )
+         NB2 = ILAENV( 1, 'SORMQR', ' ', N, N, N, -1 )
+         NB3 = ILAENV( 1, 'SORGQR', ' ', N, N, N, -1 )
+         NB = MAX( NB1, NB2, NB3 )
+         LOPT = 2*N+N*(NB+1)
+         WORK( 1 ) = LOPT
+      END IF
+*
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'SGEGS ', -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 = SLANGE( 'M', N, N, A, LDA, WORK )
+      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 SLASCL( '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 = SLANGE( 'M', N, N, B, LDB, WORK )
+      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 SLASCL( '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
+*     Workspace layout:  (2*N words -- "work..." not actually used)
+*        left_permutation, right_permutation, work...
+*
+      ILEFT = 1
+      IRIGHT = N + 1
+      IWORK = IRIGHT + N
+      CALL SGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
+     $             WORK( IRIGHT ), WORK( IWORK ), 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
+*     Workspace layout:  ("work..." must have at least N words)
+*        left_permutation, right_permutation, tau, work...
+*
+      IROWS = IHI + 1 - ILO
+      ICOLS = N + 1 - ILO
+      ITAU = IWORK
+      IWORK = ITAU + IROWS
+      CALL SGEQRF( 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 SORMQR( 'L', 'T', 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 SLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL )
+         CALL SLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
+     $                VSL( ILO+1, ILO ), LDVSL )
+         CALL SORGQR( 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 SLASET( 'Full', N, N, ZERO, ONE, VSR, LDVSR )
+*
+*     Reduce to generalized Hessenberg form
+*
+      CALL SGGHRD( 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
+*     Workspace layout:  ("work..." must have at least 1 word)
+*        left_permutation, right_permutation, work...
+*
+      IWORK = ITAU
+      CALL SHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
+     $             ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
+     $             WORK( IWORK ), LWORK+1-IWORK, 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 SGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
+     $                WORK( IRIGHT ), N, VSL, LDVSL, IINFO )
+         IF( IINFO.NE.0 ) THEN
+            INFO = N + 7
+            GO TO 10
+         END IF
+      END IF
+      IF( ILVSR ) THEN
+         CALL SGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
+     $                WORK( 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 SLASCL( 'H', -1, -1, ANRMTO, ANRM, N, N, A, LDA, IINFO )
+         IF( IINFO.NE.0 ) THEN
+            INFO = N + 9
+            RETURN
+         END IF
+         CALL SLASCL( 'G', -1, -1, ANRMTO, ANRM, N, 1, ALPHAR, N,
+     $                IINFO )
+         IF( IINFO.NE.0 ) THEN
+            INFO = N + 9
+            RETURN
+         END IF
+         CALL SLASCL( 'G', -1, -1, ANRMTO, ANRM, N, 1, ALPHAI, N,
+     $                IINFO )
+         IF( IINFO.NE.0 ) THEN
+            INFO = N + 9
+            RETURN
+         END IF
+      END IF
+*
+      IF( ILBSCL ) THEN
+         CALL SLASCL( 'U', -1, -1, BNRMTO, BNRM, N, N, B, LDB, IINFO )
+         IF( IINFO.NE.0 ) THEN
+            INFO = N + 9
+            RETURN
+         END IF
+         CALL SLASCL( '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 SGEGS
+*
+      END

lapack-netlib/SRC/sgegv.f

@@ -0,0 +1,766 @@
+*> \brief <b> SGEGV computes the eigenvalues and, optionally, the left and/or right eigenvectors of a real matrix pair (A,B).</b>
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download SGEGV + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgegv.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgegv.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgegv.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE SGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
+*                         BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
+*
+*       .. Scalar Arguments ..
+*       CHARACTER          JOBVL, JOBVR
+*       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
+*       ..
+*       .. Array Arguments ..
+*       REAL               A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
+*      $                   B( LDB, * ), BETA( * ), VL( LDVL, * ),
+*      $                   VR( LDVR, * ), WORK( * )
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> This routine is deprecated and has been replaced by routine SGGEV.
+*>
+*> SGEGV computes the eigenvalues and, optionally, the left and/or right
+*> eigenvectors of a real matrix pair (A,B).
+*> Given two square matrices A and B,
+*> the generalized nonsymmetric eigenvalue problem (GNEP) is to find the
+*> eigenvalues lambda and corresponding (non-zero) eigenvectors x such
+*> that
+*>
+*>    A*x = lambda*B*x.
+*>
+*> An alternate form is to find the eigenvalues mu and corresponding
+*> eigenvectors y such that
+*>
+*>    mu*A*y = B*y.
+*>
+*> These two forms are equivalent with mu = 1/lambda and x = y if
+*> neither lambda nor mu is zero.  In order to deal with the case that
+*> lambda or mu is zero or small, two values alpha and beta are returned
+*> for each eigenvalue, such that lambda = alpha/beta and
+*> mu = beta/alpha.
+*>
+*> The vectors x and y in the above equations are right eigenvectors of
+*> the matrix pair (A,B).  Vectors u and v satisfying
+*>
+*>    u**H*A = lambda*u**H*B  or  mu*v**H*A = v**H*B
+*>
+*> are left eigenvectors of (A,B).
+*>
+*> Note: this routine performs "full balancing" on A and B
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] JOBVL
+*> \verbatim
+*>          JOBVL is CHARACTER*1
+*>          = 'N':  do not compute the left generalized eigenvectors;
+*>          = 'V':  compute the left generalized eigenvectors (returned
+*>                  in VL).
+*> \endverbatim
+*>
+*> \param[in] JOBVR
+*> \verbatim
+*>          JOBVR is CHARACTER*1
+*>          = 'N':  do not compute the right generalized eigenvectors;
+*>          = 'V':  compute the right generalized eigenvectors (returned
+*>                  in VR).
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the matrices A, B, VL, and VR.  N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*>          A is REAL array, dimension (LDA, N)
+*>          On entry, the matrix A.
+*>          If JOBVL = 'V' or JOBVR = 'V', then on exit A
+*>          contains the real Schur form of A from the generalized Schur
+*>          factorization of the pair (A,B) after balancing.
+*>          If no eigenvectors were computed, then only the diagonal
+*>          blocks from the Schur form will be correct.  See SGGHRD and
+*>          SHGEQZ for details.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>          The leading dimension of A.  LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] B
+*> \verbatim
+*>          B is REAL array, dimension (LDB, N)
+*>          On entry, the matrix B.
+*>          If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the
+*>          upper triangular matrix obtained from B in the generalized
+*>          Schur factorization of the pair (A,B) after balancing.
+*>          If no eigenvectors were computed, then only those elements of
+*>          B corresponding to the diagonal blocks from the Schur form of
+*>          A will be correct.  See SGGHRD and SHGEQZ for details.
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*>          LDB is INTEGER
+*>          The leading dimension of B.  LDB >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] ALPHAR
+*> \verbatim
+*>          ALPHAR is REAL array, dimension (N)
+*>          The real parts of each scalar alpha defining an eigenvalue of
+*>          GNEP.
+*> \endverbatim
+*>
+*> \param[out] ALPHAI
+*> \verbatim
+*>          ALPHAI is REAL array, dimension (N)
+*>          The imaginary parts of each scalar alpha defining an
+*>          eigenvalue of GNEP.  If ALPHAI(j) is zero, then the j-th
+*>          eigenvalue is real; if positive, then the j-th and
+*>          (j+1)-st eigenvalues are a complex conjugate pair, with
+*>          ALPHAI(j+1) = -ALPHAI(j).
+*> \endverbatim
+*>
+*> \param[out] BETA
+*> \verbatim
+*>          BETA is REAL array, dimension (N)
+*>          The scalars beta that define the eigenvalues of GNEP.
+*>
+*>          Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
+*>          beta = BETA(j) represent the j-th eigenvalue of the matrix
+*>          pair (A,B), in one of the forms lambda = alpha/beta or
+*>          mu = beta/alpha.  Since either lambda or mu may overflow,
+*>          they should not, in general, be computed.
+*> \endverbatim
+*>
+*> \param[out] VL
+*> \verbatim
+*>          VL is REAL array, dimension (LDVL,N)
+*>          If JOBVL = 'V', the left eigenvectors u(j) are stored
+*>          in the columns of VL, in the same order as their eigenvalues.
+*>          If the j-th eigenvalue is real, then u(j) = VL(:,j).
+*>          If the j-th and (j+1)-st eigenvalues form a complex conjugate
+*>          pair, then
+*>             u(j) = VL(:,j) + i*VL(:,j+1)
+*>          and
+*>            u(j+1) = VL(:,j) - i*VL(:,j+1).
+*>
+*>          Each eigenvector is scaled so that its largest component has
+*>          abs(real part) + abs(imag. part) = 1, except for eigenvectors
+*>          corresponding to an eigenvalue with alpha = beta = 0, which
+*>          are set to zero.
+*>          Not referenced if JOBVL = 'N'.
+*> \endverbatim
+*>
+*> \param[in] LDVL
+*> \verbatim
+*>          LDVL is INTEGER
+*>          The leading dimension of the matrix VL. LDVL >= 1, and
+*>          if JOBVL = 'V', LDVL >= N.
+*> \endverbatim
+*>
+*> \param[out] VR
+*> \verbatim
+*>          VR is REAL array, dimension (LDVR,N)
+*>          If JOBVR = 'V', the right eigenvectors x(j) are stored
+*>          in the columns of VR, in the same order as their eigenvalues.
+*>          If the j-th eigenvalue is real, then x(j) = VR(:,j).
+*>          If the j-th and (j+1)-st eigenvalues form a complex conjugate
+*>          pair, then
+*>            x(j) = VR(:,j) + i*VR(:,j+1)
+*>          and
+*>            x(j+1) = VR(:,j) - i*VR(:,j+1).
+*>
+*>          Each eigenvector is scaled so that its largest component has
+*>          abs(real part) + abs(imag. part) = 1, except for eigenvalues
+*>          corresponding to an eigenvalue with alpha = beta = 0, which
+*>          are set to zero.
+*>          Not referenced if JOBVR = 'N'.
+*> \endverbatim
+*>
+*> \param[in] LDVR
+*> \verbatim
+*>          LDVR is INTEGER
+*>          The leading dimension of the matrix VR. LDVR >= 1, and
+*>          if JOBVR = 'V', LDVR >= N.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is REAL array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*>          LWORK is INTEGER
+*>          The dimension of the array WORK.  LWORK >= max(1,8*N).
+*>          For good performance, LWORK must generally be larger.
+*>          To compute the optimal value of LWORK, call ILAENV to get
+*>          blocksizes (for SGEQRF, SORMQR, and SORGQR.)  Then compute:
+*>          NB  -- MAX of the blocksizes for SGEQRF, SORMQR, and SORGQR;
+*>          The optimal LWORK is:
+*>              2*N + MAX( 6*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.
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*>          INFO is 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 ALPHAR(j), ALPHAI(j), and BETA(j)
+*>                should be correct for j=INFO+1,...,N.
+*>          > N:  errors that usually indicate LAPACK problems:
+*>                =N+1: error return from SGGBAL
+*>                =N+2: error return from SGEQRF
+*>                =N+3: error return from SORMQR
+*>                =N+4: error return from SORGQR
+*>                =N+5: error return from SGGHRD
+*>                =N+6: error return from SHGEQZ (other than failed
+*>                                                iteration)
+*>                =N+7: error return from STGEVC
+*>                =N+8: error return from SGGBAK (computing VL)
+*>                =N+9: error return from SGGBAK (computing VR)
+*>                =N+10: error return from SLASCL (various calls)
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \ingroup realGEeigen
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>
+*>  Balancing
+*>  ---------
+*>
+*>  This driver calls SGGBAL 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, SGGBAK 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 real Schur
+*>  form[*] of the "balanced" versions of A and B.  If no eigenvectors
+*>  are computed, then only the diagonal blocks will be correct.
+*>
+*>  [*] See SHGEQZ, SGEGS, or read the book "Matrix Computations",
+*>      by Golub & van Loan, pub. by Johns Hopkins U. Press.
+*> \endverbatim
+*>
+*  =====================================================================
+      SUBROUTINE SGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
+     $                  BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
+*
+*  -- LAPACK driver routine --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      CHARACTER          JOBVL, JOBVR
+      INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
+*     ..
+*     .. Array Arguments ..
+      REAL               A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
+     $                   B( LDB, * ), BETA( * ), VL( LDVL, * ),
+     $                   VR( LDVR, * ), WORK( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ZERO, ONE
+      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            ILIMIT, ILV, ILVL, ILVR, LQUERY
+      CHARACTER          CHTEMP
+      INTEGER            ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
+     $                   IN, IRIGHT, IROWS, ITAU, IWORK, JC, JR, LOPT,
+     $                   LWKMIN, LWKOPT, NB, NB1, NB2, NB3
+      REAL               ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM,
+     $                   BNRM1, BNRM2, EPS, ONEPLS, SAFMAX, SAFMIN,
+     $                   SALFAI, SALFAR, SBETA, SCALE, TEMP
+*     ..
+*     .. Local Arrays ..
+      LOGICAL            LDUMMA( 1 )
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           SGEQRF, SGGBAK, SGGBAL, SGGHRD, SHGEQZ, SLACPY,
+     $                   SLASCL, SLASET, SORGQR, SORMQR, STGEVC, XERBLA
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      INTEGER            ILAENV
+      REAL               SLAMCH, SLANGE
+      EXTERNAL           ILAENV, LSAME, SLAMCH, SLANGE
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, INT, MAX
+*     ..
+*     .. 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( 8*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 = -12
+      ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
+         INFO = -14
+      ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
+         INFO = -16
+      END IF
+*
+      IF( INFO.EQ.0 ) THEN
+         NB1 = ILAENV( 1, 'SGEQRF', ' ', N, N, -1, -1 )
+         NB2 = ILAENV( 1, 'SORMQR', ' ', N, N, N, -1 )
+         NB3 = ILAENV( 1, 'SORGQR', ' ', N, N, N, -1 )
+         NB = MAX( NB1, NB2, NB3 )
+         LOPT = 2*N + MAX( 6*N, N*(NB+1) )
+         WORK( 1 ) = LOPT
+      END IF
+*
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'SGEGV ', -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
+      ONEPLS = ONE + ( 4*EPS )
+*
+*     Scale A
+*
+      ANRM = SLANGE( 'M', N, N, A, LDA, WORK )
+      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 SLASCL( '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 = SLANGE( 'M', N, N, B, LDB, WORK )
+      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 SLASCL( '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
+*     Workspace layout:  (8*N words -- "work" requires 6*N words)
+*        left_permutation, right_permutation, work...
+*
+      ILEFT = 1
+      IRIGHT = N + 1
+      IWORK = IRIGHT + N
+      CALL SGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
+     $             WORK( IRIGHT ), WORK( IWORK ), IINFO )
+      IF( IINFO.NE.0 ) THEN
+         INFO = N + 1
+         GO TO 120
+      END IF
+*
+*     Reduce B to triangular form, and initialize VL and/or VR
+*     Workspace layout:  ("work..." must have at least N words)
+*        left_permutation, right_permutation, tau, work...
+*
+      IROWS = IHI + 1 - ILO
+      IF( ILV ) THEN
+         ICOLS = N + 1 - ILO
+      ELSE
+         ICOLS = IROWS
+      END IF
+      ITAU = IWORK
+      IWORK = ITAU + IROWS
+      CALL SGEQRF( 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 120
+      END IF
+*
+      CALL SORMQR( 'L', 'T', 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 120
+      END IF
+*
+      IF( ILVL ) THEN
+         CALL SLASET( 'Full', N, N, ZERO, ONE, VL, LDVL )
+         CALL SLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
+     $                VL( ILO+1, ILO ), LDVL )
+         CALL SORGQR( 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 120
+         END IF
+      END IF
+*
+      IF( ILVR )
+     $   CALL SLASET( 'Full', N, N, ZERO, ONE, VR, LDVR )
+*
+*     Reduce to generalized Hessenberg form
+*
+      IF( ILV ) THEN
+*
+*        Eigenvectors requested -- work on whole matrix.
+*
+         CALL SGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
+     $                LDVL, VR, LDVR, IINFO )
+      ELSE
+         CALL SGGHRD( '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 120
+      END IF
+*
+*     Perform QZ algorithm
+*     Workspace layout:  ("work..." must have at least 1 word)
+*        left_permutation, right_permutation, work...
+*
+      IWORK = ITAU
+      IF( ILV ) THEN
+         CHTEMP = 'S'
+      ELSE
+         CHTEMP = 'E'
+      END IF
+      CALL SHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
+     $             ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
+     $             WORK( IWORK ), LWORK+1-IWORK, 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 120
+      END IF
+*
+      IF( ILV ) THEN
+*
+*        Compute Eigenvectors  (STGEVC requires 6*N words of workspace)
+*
+         IF( ILVL ) THEN
+            IF( ILVR ) THEN
+               CHTEMP = 'B'
+            ELSE
+               CHTEMP = 'L'
+            END IF
+         ELSE
+            CHTEMP = 'R'
+         END IF
+*
+         CALL STGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
+     $                VR, LDVR, N, IN, WORK( IWORK ), IINFO )
+         IF( IINFO.NE.0 ) THEN
+            INFO = N + 7
+            GO TO 120
+         END IF
+*
+*        Undo balancing on VL and VR, rescale
+*
+         IF( ILVL ) THEN
+            CALL SGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
+     $                   WORK( IRIGHT ), N, VL, LDVL, IINFO )
+            IF( IINFO.NE.0 ) THEN
+               INFO = N + 8
+               GO TO 120
+            END IF
+            DO 50 JC = 1, N
+               IF( ALPHAI( JC ).LT.ZERO )
+     $            GO TO 50
+               TEMP = ZERO
+               IF( ALPHAI( JC ).EQ.ZERO ) THEN
+                  DO 10 JR = 1, N
+                     TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
+   10             CONTINUE
+               ELSE
+                  DO 20 JR = 1, N
+                     TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
+     $                      ABS( VL( JR, JC+1 ) ) )
+   20             CONTINUE
+               END IF
+               IF( TEMP.LT.SAFMIN )
+     $            GO TO 50
+               TEMP = ONE / TEMP
+               IF( ALPHAI( JC ).EQ.ZERO ) THEN
+                  DO 30 JR = 1, N
+                     VL( JR, JC ) = VL( JR, JC )*TEMP
+   30             CONTINUE
+               ELSE
+                  DO 40 JR = 1, N
+                     VL( JR, JC ) = VL( JR, JC )*TEMP
+                     VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
+   40             CONTINUE
+               END IF
+   50       CONTINUE
+         END IF
+         IF( ILVR ) THEN
+            CALL SGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
+     $                   WORK( IRIGHT ), N, VR, LDVR, IINFO )
+            IF( IINFO.NE.0 ) THEN
+               INFO = N + 9
+               GO TO 120
+            END IF
+            DO 100 JC = 1, N
+               IF( ALPHAI( JC ).LT.ZERO )
+     $            GO TO 100
+               TEMP = ZERO
+               IF( ALPHAI( JC ).EQ.ZERO ) THEN
+                  DO 60 JR = 1, N
+                     TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
+   60             CONTINUE
+               ELSE
+                  DO 70 JR = 1, N
+                     TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
+     $                      ABS( VR( JR, JC+1 ) ) )
+   70             CONTINUE
+               END IF
+               IF( TEMP.LT.SAFMIN )
+     $            GO TO 100
+               TEMP = ONE / TEMP
+               IF( ALPHAI( JC ).EQ.ZERO ) THEN
+                  DO 80 JR = 1, N
+                     VR( JR, JC ) = VR( JR, JC )*TEMP
+   80             CONTINUE
+               ELSE
+                  DO 90 JR = 1, N
+                     VR( JR, JC ) = VR( JR, JC )*TEMP
+                     VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
+   90             CONTINUE
+               END IF
+  100       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 110 JC = 1, N
+         ABSAR = ABS( ALPHAR( JC ) )
+         ABSAI = ABS( ALPHAI( JC ) )
+         ABSB = ABS( BETA( JC ) )
+         SALFAR = ANRM*ALPHAR( JC )
+         SALFAI = ANRM*ALPHAI( JC )
+         SBETA = BNRM*BETA( JC )
+         ILIMIT = .FALSE.
+         SCALE = ONE
+*
+*        Check for significant underflow in ALPHAI
+*
+         IF( ABS( SALFAI ).LT.SAFMIN .AND. ABSAI.GE.
+     $       MAX( SAFMIN, EPS*ABSAR, EPS*ABSB ) ) THEN
+            ILIMIT = .TRUE.
+            SCALE = ( ONEPLS*SAFMIN / ANRM1 ) /
+     $              MAX( ONEPLS*SAFMIN, ANRM2*ABSAI )
+*
+         ELSE IF( SALFAI.EQ.ZERO ) THEN
+*
+*           If insignificant underflow in ALPHAI, then make the
+*           conjugate eigenvalue real.
+*
+            IF( ALPHAI( JC ).LT.ZERO .AND. JC.GT.1 ) THEN
+               ALPHAI( JC-1 ) = ZERO
+            ELSE IF( ALPHAI( JC ).GT.ZERO .AND. JC.LT.N ) THEN
+               ALPHAI( JC+1 ) = ZERO
+            END IF
+         END IF
+*
+*        Check for significant underflow in ALPHAR
+*
+         IF( ABS( SALFAR ).LT.SAFMIN .AND. ABSAR.GE.
+     $       MAX( SAFMIN, EPS*ABSAI, EPS*ABSB ) ) THEN
+            ILIMIT = .TRUE.
+            SCALE = MAX( SCALE, ( ONEPLS*SAFMIN / ANRM1 ) /
+     $              MAX( ONEPLS*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, ( ONEPLS*SAFMIN / BNRM1 ) /
+     $              MAX( ONEPLS*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 ALPHAR, ALPHAI, BETA if necessary.
+*
+         IF( ILIMIT ) THEN
+            SALFAR = ( SCALE*ALPHAR( JC ) )*ANRM
+            SALFAI = ( SCALE*ALPHAI( JC ) )*ANRM
+            SBETA = ( SCALE*BETA( JC ) )*BNRM
+         END IF
+         ALPHAR( JC ) = SALFAR
+         ALPHAI( JC ) = SALFAI
+         BETA( JC ) = SBETA
+  110 CONTINUE
+*
+  120 CONTINUE
+      WORK( 1 ) = LWKOPT
+*
+      RETURN
+*
+*     End of SGEGV
+*
+      END

lapack-netlib/SRC/zgegs.f

@@ -0,0 +1,528 @@
+*> \brief <b> ZGEGS computes the eigenvalues, Schur form, and, optionally, the left and or/right Schur vectors of a complex matrix pair (A,B)</b>
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download ZGEGS + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgegs.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgegs.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgegs.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE ZGEGS( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHA, BETA,
+*                         VSL, LDVSL, VSR, LDVSR, WORK, LWORK, RWORK,
+*                         INFO )
+*
+*       .. Scalar Arguments ..
+*       CHARACTER          JOBVSL, JOBVSR
+*       INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
+*       ..
+*       .. Array Arguments ..
+*       DOUBLE PRECISION   RWORK( * )
+*       COMPLEX*16         A( LDA, * ), ALPHA( * ), B( LDB, * ),
+*      $                   BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
+*      $                   WORK( * )
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> This routine is deprecated and has been replaced by routine ZGGES.
+*>
+*> ZGEGS computes the eigenvalues, Schur form, and, optionally, the
+*> left and or/right Schur vectors of a complex matrix pair (A,B).
+*> Given two square matrices A and B, the generalized Schur
+*> factorization has the form
+*>
+*>    A = Q*S*Z**H,  B = Q*T*Z**H
+*>
+*> where Q and Z are unitary matrices and S and T are upper triangular.
+*> The columns of Q are the left Schur vectors
+*> and the columns of Z are the right Schur vectors.
+*>
+*> If only the eigenvalues of (A,B) are needed, the driver routine
+*> ZGEGV should be used instead.  See ZGEGV for a description of the
+*> eigenvalues of the generalized nonsymmetric eigenvalue problem
+*> (GNEP).
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] JOBVSL
+*> \verbatim
+*>          JOBVSL is CHARACTER*1
+*>          = 'N':  do not compute the left Schur vectors;
+*>          = 'V':  compute the left Schur vectors (returned in VSL).
+*> \endverbatim
+*>
+*> \param[in] JOBVSR
+*> \verbatim
+*>          JOBVSR is CHARACTER*1
+*>          = 'N':  do not compute the right Schur vectors;
+*>          = 'V':  compute the right Schur vectors (returned in VSR).
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the matrices A, B, VSL, and VSR.  N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*>          A is COMPLEX*16 array, dimension (LDA, N)
+*>          On entry, the matrix A.
+*>          On exit, the upper triangular matrix S from the generalized
+*>          Schur factorization.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>          The leading dimension of A.  LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] B
+*> \verbatim
+*>          B is COMPLEX*16 array, dimension (LDB, N)
+*>          On entry, the matrix B.
+*>          On exit, the upper triangular matrix T from the generalized
+*>          Schur factorization.
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*>          LDB is INTEGER
+*>          The leading dimension of B.  LDB >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] ALPHA
+*> \verbatim
+*>          ALPHA is COMPLEX*16 array, dimension (N)
+*>          The complex scalars alpha that define the eigenvalues of
+*>          GNEP.  ALPHA(j) = S(j,j), the diagonal element of the Schur
+*>          form of A.
+*> \endverbatim
+*>
+*> \param[out] BETA
+*> \verbatim
+*>          BETA is COMPLEX*16 array, dimension (N)
+*>          The non-negative real scalars beta that define the
+*>          eigenvalues of GNEP.  BETA(j) = T(j,j), the diagonal element
+*>          of the triangular factor T.
+*>
+*>          Together, the quantities alpha = ALPHA(j) and beta = BETA(j)
+*>          represent the j-th eigenvalue of the matrix pair (A,B), in
+*>          one of the forms lambda = alpha/beta or mu = beta/alpha.
+*>          Since either lambda or mu may overflow, they should not,
+*>          in general, be computed.
+*> \endverbatim
+*>
+*> \param[out] VSL
+*> \verbatim
+*>          VSL is COMPLEX*16 array, dimension (LDVSL,N)
+*>          If JOBVSL = 'V', the matrix of left Schur vectors Q.
+*>          Not referenced if JOBVSL = 'N'.
+*> \endverbatim
+*>
+*> \param[in] LDVSL
+*> \verbatim
+*>          LDVSL is INTEGER
+*>          The leading dimension of the matrix VSL. LDVSL >= 1, and
+*>          if JOBVSL = 'V', LDVSL >= N.
+*> \endverbatim
+*>
+*> \param[out] VSR
+*> \verbatim
+*>          VSR is COMPLEX*16 array, dimension (LDVSR,N)
+*>          If JOBVSR = 'V', the matrix of right Schur vectors Z.
+*>          Not referenced if JOBVSR = 'N'.
+*> \endverbatim
+*>
+*> \param[in] LDVSR
+*> \verbatim
+*>          LDVSR is INTEGER
+*>          The leading dimension of the matrix VSR. LDVSR >= 1, and
+*>          if JOBVSR = 'V', LDVSR >= N.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*>          LWORK is 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 ZGEQRF, ZUNMQR, and CUNGQR.)  Then compute:
+*>          NB  -- MAX of the blocksizes for ZGEQRF, ZUNMQR, 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.
+*> \endverbatim
+*>
+*> \param[out] RWORK
+*> \verbatim
+*>          RWORK is DOUBLE PRECISION array, dimension (3*N)
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*>          INFO is 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 ZGGBAL
+*>                =N+2: error return from ZGEQRF
+*>                =N+3: error return from ZUNMQR
+*>                =N+4: error return from ZUNGQR
+*>                =N+5: error return from ZGGHRD
+*>                =N+6: error return from ZHGEQZ (other than failed
+*>                                               iteration)
+*>                =N+7: error return from ZGGBAK (computing VSL)
+*>                =N+8: error return from ZGGBAK (computing VSR)
+*>                =N+9: error return from ZLASCL (various places)
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \ingroup complex16GEeigen
+*
+*  =====================================================================
+      SUBROUTINE ZGEGS( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHA, BETA,
+     $                  VSL, LDVSL, VSR, LDVSR, WORK, LWORK, RWORK,
+     $                  INFO )
+*
+*  -- LAPACK driver routine --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      CHARACTER          JOBVSL, JOBVSR
+      INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   RWORK( * )
+      COMPLEX*16         A( LDA, * ), ALPHA( * ), B( LDB, * ),
+     $                   BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
+     $                   WORK( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
+      COMPLEX*16         CZERO, CONE
+      PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
+     $                   CONE = ( 1.0D0, 0.0D0 ) )
+*     ..
+*     .. 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
+      DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
+     $                   SAFMIN, SMLNUM
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD, ZHGEQZ,
+     $                   ZLACPY, ZLASCL, ZLASET, ZUNGQR, ZUNMQR
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      INTEGER            ILAENV
+      DOUBLE PRECISION   DLAMCH, ZLANGE
+      EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANGE
+*     ..
+*     .. 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, 'ZGEQRF', ' ', N, N, -1, -1 )
+         NB2 = ILAENV( 1, 'ZUNMQR', ' ', N, N, N, -1 )
+         NB3 = ILAENV( 1, 'ZUNGQR', ' ', 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( 'ZGEGS ', -INFO )
+         RETURN
+      ELSE IF( LQUERY ) THEN
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
+*     Get machine constants
+*
+      EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
+      SAFMIN = DLAMCH( 'S' )
+      SMLNUM = N*SAFMIN / EPS
+      BIGNUM = ONE / SMLNUM
+*
+*     Scale A if max element outside range [SMLNUM,BIGNUM]
+*
+      ANRM = ZLANGE( '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 ZLASCL( '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 = ZLANGE( '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 ZLASCL( '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 ZGGBAL( '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 ZGEQRF( 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 ZUNMQR( '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 ZLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL )
+         CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
+     $                VSL( ILO+1, ILO ), LDVSL )
+         CALL ZUNGQR( 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 ZLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR )
+*
+*     Reduce to generalized Hessenberg form
+*
+      CALL ZGGHRD( 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 ZHGEQZ( '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 ZGGBAK( '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 ZGGBAK( '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 ZLASCL( 'U', -1, -1, ANRMTO, ANRM, N, N, A, LDA, IINFO )
+         IF( IINFO.NE.0 ) THEN
+            INFO = N + 9
+            RETURN
+         END IF
+         CALL ZLASCL( '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 ZLASCL( 'U', -1, -1, BNRMTO, BNRM, N, N, B, LDB, IINFO )
+         IF( IINFO.NE.0 ) THEN
+            INFO = N + 9
+            RETURN
+         END IF
+         CALL ZLASCL( '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 ZGEGS
+*
+      END

lapack-netlib/SRC/zgegv.f

@@ -0,0 +1,703 @@
+*> \brief <b> ZGEGV computes the eigenvalues and, optionally, the left and/or right eigenvectors of a complex matrix pair (A,B).</b>
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download ZGEGV + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgegv.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgegv.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgegv.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE ZGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
+*                         VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
+*
+*       .. Scalar Arguments ..
+*       CHARACTER          JOBVL, JOBVR
+*       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
+*       ..
+*       .. Array Arguments ..
+*       DOUBLE PRECISION   RWORK( * )
+*       COMPLEX*16         A( LDA, * ), ALPHA( * ), B( LDB, * ),
+*      $                   BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
+*      $                   WORK( * )
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> This routine is deprecated and has been replaced by routine ZGGEV.
+*>
+*> ZGEGV computes the eigenvalues and, optionally, the left and/or right
+*> eigenvectors of a complex matrix pair (A,B).
+*> Given two square matrices A and B,
+*> the generalized nonsymmetric eigenvalue problem (GNEP) is to find the
+*> eigenvalues lambda and corresponding (non-zero) eigenvectors x such
+*> that
+*>    A*x = lambda*B*x.
+*>
+*> An alternate form is to find the eigenvalues mu and corresponding
+*> eigenvectors y such that
+*>    mu*A*y = B*y.
+*>
+*> These two forms are equivalent with mu = 1/lambda and x = y if
+*> neither lambda nor mu is zero.  In order to deal with the case that
+*> lambda or mu is zero or small, two values alpha and beta are returned
+*> for each eigenvalue, such that lambda = alpha/beta and
+*> mu = beta/alpha.
+*>
+*> The vectors x and y in the above equations are right eigenvectors of
+*> the matrix pair (A,B).  Vectors u and v satisfying
+*>    u**H*A = lambda*u**H*B  or  mu*v**H*A = v**H*B
+*> are left eigenvectors of (A,B).
+*>
+*> Note: this routine performs "full balancing" on A and B
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] JOBVL
+*> \verbatim
+*>          JOBVL is CHARACTER*1
+*>          = 'N':  do not compute the left generalized eigenvectors;
+*>          = 'V':  compute the left generalized eigenvectors (returned
+*>                  in VL).
+*> \endverbatim
+*>
+*> \param[in] JOBVR
+*> \verbatim
+*>          JOBVR is CHARACTER*1
+*>          = 'N':  do not compute the right generalized eigenvectors;
+*>          = 'V':  compute the right generalized eigenvectors (returned
+*>                  in VR).
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the matrices A, B, VL, and VR.  N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*>          A is COMPLEX*16 array, dimension (LDA, N)
+*>          On entry, the matrix A.
+*>          If JOBVL = 'V' or JOBVR = 'V', then on exit A
+*>          contains the Schur form of A from the generalized Schur
+*>          factorization of the pair (A,B) after balancing.  If no
+*>          eigenvectors were computed, then only the diagonal elements
+*>          of the Schur form will be correct.  See ZGGHRD and ZHGEQZ
+*>          for details.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>          The leading dimension of A.  LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] B
+*> \verbatim
+*>          B is COMPLEX*16 array, dimension (LDB, N)
+*>          On entry, the matrix B.
+*>          If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the
+*>          upper triangular matrix obtained from B in the generalized
+*>          Schur factorization of the pair (A,B) after balancing.
+*>          If no eigenvectors were computed, then only the diagonal
+*>          elements of B will be correct.  See ZGGHRD and ZHGEQZ for
+*>          details.
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*>          LDB is INTEGER
+*>          The leading dimension of B.  LDB >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] ALPHA
+*> \verbatim
+*>          ALPHA is COMPLEX*16 array, dimension (N)
+*>          The complex scalars alpha that define the eigenvalues of
+*>          GNEP.
+*> \endverbatim
+*>
+*> \param[out] BETA
+*> \verbatim
+*>          BETA is COMPLEX*16 array, dimension (N)
+*>          The complex scalars beta that define the eigenvalues of GNEP.
+*>
+*>          Together, the quantities alpha = ALPHA(j) and beta = BETA(j)
+*>          represent the j-th eigenvalue of the matrix pair (A,B), in
+*>          one of the forms lambda = alpha/beta or mu = beta/alpha.
+*>          Since either lambda or mu may overflow, they should not,
+*>          in general, be computed.
+*> \endverbatim
+*>
+*> \param[out] VL
+*> \verbatim
+*>          VL is COMPLEX*16 array, dimension (LDVL,N)
+*>          If JOBVL = 'V', the left eigenvectors u(j) are stored
+*>          in the columns of VL, in the same order as their eigenvalues.
+*>          Each eigenvector is scaled so that its largest component has
+*>          abs(real part) + abs(imag. part) = 1, except for eigenvectors
+*>          corresponding to an eigenvalue with alpha = beta = 0, which
+*>          are set to zero.
+*>          Not referenced if JOBVL = 'N'.
+*> \endverbatim
+*>
+*> \param[in] LDVL
+*> \verbatim
+*>          LDVL is INTEGER
+*>          The leading dimension of the matrix VL. LDVL >= 1, and
+*>          if JOBVL = 'V', LDVL >= N.
+*> \endverbatim
+*>
+*> \param[out] VR
+*> \verbatim
+*>          VR is COMPLEX*16 array, dimension (LDVR,N)
+*>          If JOBVR = 'V', the right eigenvectors x(j) are stored
+*>          in the columns of VR, in the same order as their eigenvalues.
+*>          Each eigenvector is scaled so that its largest component has
+*>          abs(real part) + abs(imag. part) = 1, except for eigenvectors
+*>          corresponding to an eigenvalue with alpha = beta = 0, which
+*>          are set to zero.
+*>          Not referenced if JOBVR = 'N'.
+*> \endverbatim
+*>
+*> \param[in] LDVR
+*> \verbatim
+*>          LDVR is INTEGER
+*>          The leading dimension of the matrix VR. LDVR >= 1, and
+*>          if JOBVR = 'V', LDVR >= N.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*>          LWORK is 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 ZGEQRF, ZUNMQR, and ZUNGQR.)  Then compute:
+*>          NB  -- MAX of the blocksizes for ZGEQRF, ZUNMQR, and ZUNGQR;
+*>          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.
+*> \endverbatim
+*>
+*> \param[out] RWORK
+*> \verbatim
+*>          RWORK is DOUBLE PRECISION array, dimension (8*N)
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*>          INFO is 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 ZGGBAL
+*>                =N+2: error return from ZGEQRF
+*>                =N+3: error return from ZUNMQR
+*>                =N+4: error return from ZUNGQR
+*>                =N+5: error return from ZGGHRD
+*>                =N+6: error return from ZHGEQZ (other than failed
+*>                                               iteration)
+*>                =N+7: error return from ZTGEVC
+*>                =N+8: error return from ZGGBAK (computing VL)
+*>                =N+9: error return from ZGGBAK (computing VR)
+*>                =N+10: error return from ZLASCL (various calls)
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \ingroup complex16GEeigen
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>
+*>  Balancing
+*>  ---------
+*>
+*>  This driver calls ZGGBAL 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, ZGGBAK 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.
+*> \endverbatim
+*>
+*  =====================================================================
+      SUBROUTINE ZGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
+     $                  VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
+*
+*  -- LAPACK driver routine --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      CHARACTER          JOBVL, JOBVR
+      INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   RWORK( * )
+      COMPLEX*16         A( LDA, * ), ALPHA( * ), B( LDB, * ),
+     $                   BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
+     $                   WORK( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
+      COMPLEX*16         CZERO, CONE
+      PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
+     $                   CONE = ( 1.0D0, 0.0D0 ) )
+*     ..
+*     .. 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
+      DOUBLE PRECISION   ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM,
+     $                   BNRM1, BNRM2, EPS, SAFMAX, SAFMIN, SALFAI,
+     $                   SALFAR, SBETA, SCALE, TEMP
+      COMPLEX*16         X
+*     ..
+*     .. Local Arrays ..
+      LOGICAL            LDUMMA( 1 )
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD, ZHGEQZ,
+     $                   ZLACPY, ZLASCL, ZLASET, ZTGEVC, ZUNGQR, ZUNMQR
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      INTEGER            ILAENV
+      DOUBLE PRECISION   DLAMCH, ZLANGE
+      EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANGE
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, DBLE, DCMPLX, DIMAG, INT, MAX
+*     ..
+*     .. Statement Functions ..
+      DOUBLE PRECISION   ABS1
+*     ..
+*     .. Statement Function definitions ..
+      ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( 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, 'ZGEQRF', ' ', N, N, -1, -1 )
+         NB2 = ILAENV( 1, 'ZUNMQR', ' ', N, N, N, -1 )
+         NB3 = ILAENV( 1, 'ZUNGQR', ' ', 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( 'ZGEGV ', -INFO )
+         RETURN
+      ELSE IF( LQUERY ) THEN
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
+*     Get machine constants
+*
+      EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
+      SAFMIN = DLAMCH( 'S' )
+      SAFMIN = SAFMIN + SAFMIN
+      SAFMAX = ONE / SAFMIN
+*
+*     Scale A
+*
+      ANRM = ZLANGE( '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 ZLASCL( '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 = ZLANGE( '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 ZLASCL( '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 ZGGBAL( '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 ZGEQRF( 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 ZUNMQR( '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 ZLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
+         CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
+     $                VL( ILO+1, ILO ), LDVL )
+         CALL ZUNGQR( 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 ZLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
+*
+*     Reduce to generalized Hessenberg form
+*
+      IF( ILV ) THEN
+*
+*        Eigenvectors requested -- work on whole matrix.
+*
+         CALL ZGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
+     $                LDVL, VR, LDVR, IINFO )
+      ELSE
+         CALL ZGGHRD( '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 ZHGEQZ( 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 ZTGEVC( 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 ZGGBAK( '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 ZGGBAK( '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( DBLE( ALPHA( JC ) ) )
+         ABSAI = ABS( DIMAG( ALPHA( JC ) ) )
+         ABSB = ABS( DBLE( BETA( JC ) ) )
+         SALFAR = ANRM*DBLE( ALPHA( JC ) )
+         SALFAI = ANRM*DIMAG( ALPHA( JC ) )
+         SBETA = BNRM*DBLE( 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*DBLE( ALPHA( JC ) ) )*ANRM
+            SALFAI = ( SCALE*DIMAG( ALPHA( JC ) ) )*ANRM
+            SBETA = ( SCALE*BETA( JC ) )*BNRM
+         END IF
+         ALPHA( JC ) = DCMPLX( SALFAR, SALFAI )
+         BETA( JC ) = SBETA
+   70 CONTINUE
+*
+   80 CONTINUE
+      WORK( 1 ) = LWKOPT
+*
+      RETURN
+*
+*     End of ZGEGV
+*
+      END