diff --git a/lapack-netlib/SRC/cgegs.f b/lapack-netlib/SRC/cgegs.f new file mode 100644 index 000000000..b6adf9111 --- /dev/null +++ b/lapack-netlib/SRC/cgegs.f @@ -0,0 +1,528 @@ +*> \brief CGEGS computes the eigenvalues, Schur form, and, optionally, the left and or/right Schur vectors of a complex matrix pair (A,B) +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CGEGS + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \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 diff --git a/lapack-netlib/SRC/cgegv.f b/lapack-netlib/SRC/cgegv.f new file mode 100644 index 000000000..d2b254255 --- /dev/null +++ b/lapack-netlib/SRC/cgegv.f @@ -0,0 +1,703 @@ +*> \brief CGEGV computes the eigenvalues and, optionally, the left and/or right eigenvectors of a complex matrix pair (A,B). +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CGEGV + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \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 diff --git a/lapack-netlib/SRC/dgegs.f b/lapack-netlib/SRC/dgegs.f new file mode 100644 index 000000000..02e9fdcb2 --- /dev/null +++ b/lapack-netlib/SRC/dgegs.f @@ -0,0 +1,538 @@ +*> \brief DGEGS computes the eigenvalues, real Schur form, and, optionally, the left and/or right Schur vectors of a real matrix pair (A,B) +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DGEGS + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \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 diff --git a/lapack-netlib/SRC/dgegv.f b/lapack-netlib/SRC/dgegv.f new file mode 100644 index 000000000..0b5c48922 --- /dev/null +++ b/lapack-netlib/SRC/dgegv.f @@ -0,0 +1,766 @@ +*> \brief DGEGV computes the eigenvalues and, optionally, the left and/or right eigenvectors of a real matrix pair (A,B). +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DGEGV + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \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 diff --git a/lapack-netlib/SRC/sgegs.f b/lapack-netlib/SRC/sgegs.f new file mode 100644 index 000000000..11ecc67ac --- /dev/null +++ b/lapack-netlib/SRC/sgegs.f @@ -0,0 +1,538 @@ +*> \brief SGEGS computes the eigenvalues, real Schur form, and, optionally, the left and/or right Schur vectors of a real matrix pair (A,B) +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download SGEGS + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \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 diff --git a/lapack-netlib/SRC/sgegv.f b/lapack-netlib/SRC/sgegv.f new file mode 100644 index 000000000..97556e371 --- /dev/null +++ b/lapack-netlib/SRC/sgegv.f @@ -0,0 +1,766 @@ +*> \brief SGEGV computes the eigenvalues and, optionally, the left and/or right eigenvectors of a real matrix pair (A,B). +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download SGEGV + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \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 diff --git a/lapack-netlib/SRC/zgegs.f b/lapack-netlib/SRC/zgegs.f new file mode 100644 index 000000000..23f8d43d1 --- /dev/null +++ b/lapack-netlib/SRC/zgegs.f @@ -0,0 +1,528 @@ +*> \brief ZGEGS computes the eigenvalues, Schur form, and, optionally, the left and or/right Schur vectors of a complex matrix pair (A,B) +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZGEGS + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \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 diff --git a/lapack-netlib/SRC/zgegv.f b/lapack-netlib/SRC/zgegv.f new file mode 100644 index 000000000..542d3f4ff --- /dev/null +++ b/lapack-netlib/SRC/zgegv.f @@ -0,0 +1,703 @@ +*> \brief ZGEGV computes the eigenvalues and, optionally, the left and/or right eigenvectors of a complex matrix pair (A,B). +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZGEGV + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \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