From 6273f536d15680513e8cddfc4d8baa88ad2c64df Mon Sep 17 00:00:00 2001 From: "philippe.theveny" Date: Tue, 24 Feb 2015 23:50:54 +0000 Subject: Add xGGHD3: blocked Hessenberg reduction, code from Daniel Kressner. Add xGGES3 and xGGEV3: computation of the Schur form, the Schur vectors, and the generalized eigenvalues using the blocked Hessenberg reduction. --- SRC/zgges3.f | 595 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 595 insertions(+) create mode 100644 SRC/zgges3.f (limited to 'SRC/zgges3.f') diff --git a/SRC/zgges3.f b/SRC/zgges3.f new file mode 100644 index 00000000..d4455148 --- /dev/null +++ b/SRC/zgges3.f @@ -0,0 +1,595 @@ +*> \brief ZGGES3 computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices (blocked algorithm) +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZGGES3 + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZGGES3( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, +* $ LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, +* $ WORK, LWORK, RWORK, BWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER JOBVSL, JOBVSR, SORT +* INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM +* .. +* .. Array Arguments .. +* LOGICAL BWORK( * ) +* DOUBLE PRECISION RWORK( * ) +* COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ), +* $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ), +* $ WORK( * ) +* .. +* .. Function Arguments .. +* LOGICAL SELCTG +* EXTERNAL SELCTG +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZGGES3 computes for a pair of N-by-N complex nonsymmetric matrices +*> (A,B), the generalized eigenvalues, the generalized complex Schur +*> form (S, T), and optionally left and/or right Schur vectors (VSL +*> and VSR). This gives the generalized Schur factorization +*> +*> (A,B) = ( (VSL)*S*(VSR)**H, (VSL)*T*(VSR)**H ) +*> +*> where (VSR)**H is the conjugate-transpose of VSR. +*> +*> Optionally, it also orders the eigenvalues so that a selected cluster +*> of eigenvalues appears in the leading diagonal blocks of the upper +*> triangular matrix S and the upper triangular matrix T. The leading +*> columns of VSL and VSR then form an unitary basis for the +*> corresponding left and right eigenspaces (deflating subspaces). +*> +*> (If only the generalized eigenvalues are needed, use the driver +*> ZGGEV instead, which is faster.) +*> +*> A generalized eigenvalue for a pair of matrices (A,B) is a scalar w +*> or a ratio alpha/beta = w, such that A - w*B is singular. It is +*> usually represented as the pair (alpha,beta), as there is a +*> reasonable interpretation for beta=0, and even for both being zero. +*> +*> A pair of matrices (S,T) is in generalized complex Schur form if S +*> and T are upper triangular and, in addition, the diagonal elements +*> of T are non-negative real numbers. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] JOBVSL +*> \verbatim +*> JOBVSL is CHARACTER*1 +*> = 'N': do not compute the left Schur vectors; +*> = 'V': compute the left Schur vectors. +*> \endverbatim +*> +*> \param[in] JOBVSR +*> \verbatim +*> JOBVSR is CHARACTER*1 +*> = 'N': do not compute the right Schur vectors; +*> = 'V': compute the right Schur vectors. +*> \endverbatim +*> +*> \param[in] SORT +*> \verbatim +*> SORT is CHARACTER*1 +*> Specifies whether or not to order the eigenvalues on the +*> diagonal of the generalized Schur form. +*> = 'N': Eigenvalues are not ordered; +*> = 'S': Eigenvalues are ordered (see SELCTG). +*> \endverbatim +*> +*> \param[in] SELCTG +*> \verbatim +*> SELCTG is a LOGICAL FUNCTION of two COMPLEX*16 arguments +*> SELCTG must be declared EXTERNAL in the calling subroutine. +*> If SORT = 'N', SELCTG is not referenced. +*> If SORT = 'S', SELCTG is used to select eigenvalues to sort +*> to the top left of the Schur form. +*> An eigenvalue ALPHA(j)/BETA(j) is selected if +*> SELCTG(ALPHA(j),BETA(j)) is true. +*> +*> Note that a selected complex eigenvalue may no longer satisfy +*> SELCTG(ALPHA(j),BETA(j)) = .TRUE. after ordering, since +*> ordering may change the value of complex eigenvalues +*> (especially if the eigenvalue is ill-conditioned), in this +*> case INFO is set to N+2 (See INFO below). +*> \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 first of the pair of matrices. +*> On exit, A has been overwritten by its generalized Schur +*> form S. +*> \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 second of the pair of matrices. +*> On exit, B has been overwritten by its generalized Schur +*> form T. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of B. LDB >= max(1,N). +*> \endverbatim +*> +*> \param[out] SDIM +*> \verbatim +*> SDIM is INTEGER +*> If SORT = 'N', SDIM = 0. +*> If SORT = 'S', SDIM = number of eigenvalues (after sorting) +*> for which SELCTG is true. +*> \endverbatim +*> +*> \param[out] ALPHA +*> \verbatim +*> ALPHA is COMPLEX*16 array, dimension (N) +*> \endverbatim +*> +*> \param[out] BETA +*> \verbatim +*> BETA is COMPLEX*16 array, dimension (N) +*> On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the +*> generalized eigenvalues. ALPHA(j), j=1,...,N and BETA(j), +*> j=1,...,N are the diagonals of the complex Schur form (A,B) +*> output by ZGGES3. The BETA(j) will be non-negative real. +*> +*> Note: the quotients ALPHA(j)/BETA(j) may easily over- or +*> underflow, and BETA(j) may even be zero. Thus, the user +*> should avoid naively computing the ratio alpha/beta. +*> However, ALPHA will be always less than and usually +*> comparable with norm(A) in magnitude, and BETA always less +*> than and usually comparable with norm(B). +*> \endverbatim +*> +*> \param[out] VSL +*> \verbatim +*> VSL is COMPLEX*16 array, dimension (LDVSL,N) +*> If JOBVSL = 'V', VSL will contain the left Schur vectors. +*> 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', VSR will contain the right Schur vectors. +*> 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. +*> +*> 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] BWORK +*> \verbatim +*> BWORK is LOGICAL array, dimension (N) +*> Not referenced if SORT = '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: =N+1: other than QZ iteration failed in ZHGEQZ +*> =N+2: after reordering, roundoff changed values of +*> some complex eigenvalues so that leading +*> eigenvalues in the Generalized Schur form no +*> longer satisfy SELCTG=.TRUE. This could also +*> be caused due to scaling. +*> =N+3: reordering failed in ZTGSEN. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date January 2015 +* +*> \ingroup complex16GEeigen +* +* ===================================================================== + SUBROUTINE ZGGES3( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, + $ LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, + $ WORK, LWORK, RWORK, BWORK, INFO ) +* +* -- LAPACK driver routine (version 3.6.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* January 2015 +* +* .. Scalar Arguments .. + CHARACTER JOBVSL, JOBVSR, SORT + INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM +* .. +* .. Array Arguments .. + LOGICAL BWORK( * ) + DOUBLE PRECISION RWORK( * ) + COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ), + $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ), + $ WORK( * ) +* .. +* .. Function Arguments .. + LOGICAL SELCTG + EXTERNAL SELCTG +* .. +* +* ===================================================================== +* +* .. 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 CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL, + $ LQUERY, WANTST + INTEGER I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, + $ ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK, LWKOPT + DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL, + $ PVSR, SMLNUM +* .. +* .. Local Arrays .. + INTEGER IDUM( 1 ) + DOUBLE PRECISION DIF( 2 ) +* .. +* .. External Subroutines .. + EXTERNAL DLABAD, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHD3, + $ ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGSEN, ZUNGQR, + $ ZUNMQR +* .. +* .. External Functions .. + LOGICAL LSAME + DOUBLE PRECISION DLAMCH, ZLANGE + EXTERNAL LSAME, DLAMCH, ZLANGE +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, SQRT +* .. +* .. 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 +* + WANTST = LSAME( SORT, 'S' ) +* +* Test the input arguments +* + INFO = 0 + LQUERY = ( LWORK.EQ.-1 ) + IF( IJOBVL.LE.0 ) THEN + INFO = -1 + ELSE IF( IJOBVR.LE.0 ) THEN + INFO = -2 + ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN + INFO = -3 + ELSE IF( N.LT.0 ) THEN + INFO = -5 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -7 + ELSE IF( LDB.LT.MAX( 1, N ) ) THEN + INFO = -9 + ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN + INFO = -14 + ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN + INFO = -16 + ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN + INFO = -18 + END IF +* +* Compute workspace +* + IF( INFO.EQ.0 ) THEN + CALL ZGEQRF( N, N, B, LDB, WORK, WORK, -1, IERR ) + LWKOPT = MAX( 1, N + INT ( WORK( 1 ) ) ) + CALL ZUNMQR( 'L', 'C', N, N, N, B, LDB, WORK, A, LDA, WORK, + $ -1, IERR ) + LWKOPT = MAX( LWKOPT, N + INT ( WORK( 1 ) ) ) + IF( ILVSL ) THEN + CALL ZUNGQR( N, N, N, VSL, LDVSL, WORK, WORK, -1, IERR ) + LWKOPT = MAX( LWKOPT, N + INT ( WORK( 1 ) ) ) + END IF + CALL ZGGHD3( JOBVSL, JOBVSR, N, 1, N, A, LDA, B, LDB, VSL, + $ LDVSL, VSR, LDVSR, WORK, -1, IERR ) + LWKOPT = MAX( LWKOPT, N + INT ( WORK( 1 ) ) ) + CALL ZHGEQZ( 'S', JOBVSL, JOBVSR, N, 1, N, A, LDA, B, LDB, + $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK, + $ -1, RWORK, IERR ) + LWKOPT = MAX( LWKOPT, INT ( WORK( 1 ) ) ) + IF( WANTST ) THEN + CALL ZTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, + $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, SDIM, + $ PVSL, PVSR, DIF, WORK, -1, IDUM, 1, IERR ) + LWKOPT = MAX( LWKOPT, INT ( WORK( 1 ) ) ) + END IF + WORK( 1 ) = DCMPLX( WKOPT ) + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZGGES3 ', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) THEN + SDIM = 0 + RETURN + END IF +* +* Get machine constants +* + EPS = DLAMCH( 'P' ) + SMLNUM = DLAMCH( 'S' ) + BIGNUM = ONE / SMLNUM + CALL DLABAD( SMLNUM, BIGNUM ) + SMLNUM = SQRT( SMLNUM ) / 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 ) + $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR ) +* +* 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 ) + $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR ) +* +* Permute the matrix to make it more nearly triangular +* + ILEFT = 1 + IRIGHT = N + 1 + IRWRK = IRIGHT + N + CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ), + $ RWORK( IRIGHT ), RWORK( IRWRK ), IERR ) +* +* Reduce B to triangular form (QR decomposition of B) +* + IROWS = IHI + 1 - ILO + ICOLS = N + 1 - ILO + ITAU = 1 + IWRK = ITAU + IROWS + CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ), + $ WORK( IWRK ), LWORK+1-IWRK, IERR ) +* +* Apply the orthogonal transformation to matrix A +* + CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB, + $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ), + $ LWORK+1-IWRK, IERR ) +* +* Initialize VSL +* + IF( ILVSL ) THEN + CALL ZLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL ) + IF( IROWS.GT.1 ) THEN + CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB, + $ VSL( ILO+1, ILO ), LDVSL ) + END IF + CALL ZUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL, + $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR ) + END IF +* +* Initialize VSR +* + IF( ILVSR ) + $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR ) +* +* Reduce to generalized Hessenberg form +* + CALL ZGGHD3( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL, + $ LDVSL, VSR, LDVSR, WORK( IWRK ), LWORK+1-IWRK, IERR ) +* + SDIM = 0 +* +* Perform QZ algorithm, computing Schur vectors if desired +* + IWRK = ITAU + CALL ZHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, + $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWRK ), + $ LWORK+1-IWRK, RWORK( IRWRK ), IERR ) + IF( IERR.NE.0 ) THEN + IF( IERR.GT.0 .AND. IERR.LE.N ) THEN + INFO = IERR + ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN + INFO = IERR - N + ELSE + INFO = N + 1 + END IF + GO TO 30 + END IF +* +* Sort eigenvalues ALPHA/BETA if desired +* + IF( WANTST ) THEN +* +* Undo scaling on eigenvalues before selecting +* + IF( ILASCL ) + $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, 1, ALPHA, N, IERR ) + IF( ILBSCL ) + $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, 1, BETA, N, IERR ) +* +* Select eigenvalues +* + DO 10 I = 1, N + BWORK( I ) = SELCTG( ALPHA( I ), BETA( I ) ) + 10 CONTINUE +* + CALL ZTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, ALPHA, + $ BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PVSL, PVSR, + $ DIF, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1, IERR ) + IF( IERR.EQ.1 ) + $ INFO = N + 3 +* + END IF +* +* Apply back-permutation to VSL and VSR +* + IF( ILVSL ) + $ CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ), + $ RWORK( IRIGHT ), N, VSL, LDVSL, IERR ) + IF( ILVSR ) + $ CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ), + $ RWORK( IRIGHT ), N, VSR, LDVSR, IERR ) +* +* Undo scaling +* + IF( ILASCL ) THEN + CALL ZLASCL( 'U', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR ) + CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR ) + END IF +* + IF( ILBSCL ) THEN + CALL ZLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR ) + CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR ) + END IF +* + IF( WANTST ) THEN +* +* Check if reordering is correct +* + LASTSL = .TRUE. + SDIM = 0 + DO 20 I = 1, N + CURSL = SELCTG( ALPHA( I ), BETA( I ) ) + IF( CURSL ) + $ SDIM = SDIM + 1 + IF( CURSL .AND. .NOT.LASTSL ) + $ INFO = N + 2 + LASTSL = CURSL + 20 CONTINUE +* + END IF +* + 30 CONTINUE +* + WORK( 1 ) = DCMPLX( LWKOPT ) +* + RETURN +* +* End of ZGGES3 +* + END -- cgit v1.2.3