*> \brief \b ZHGEQZ * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * *> \htmlonly *> Download ZHGEQZ + dependencies *> *> [TGZ] *> *> [ZIP] *> *> [TXT] *> \endhtmlonly * * Definition: * =========== * * SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, * ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, * RWORK, INFO ) * * .. Scalar Arguments .. * CHARACTER COMPQ, COMPZ, JOB * INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N * .. * .. Array Arguments .. * DOUBLE PRECISION RWORK( * ) * COMPLEX*16 ALPHA( * ), BETA( * ), H( LDH, * ), * \$ Q( LDQ, * ), T( LDT, * ), WORK( * ), * \$ Z( LDZ, * ) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> ZHGEQZ computes the eigenvalues of a complex matrix pair (H,T), *> where H is an upper Hessenberg matrix and T is upper triangular, *> using the single-shift QZ method. *> Matrix pairs of this type are produced by the reduction to *> generalized upper Hessenberg form of a complex matrix pair (A,B): *> *> A = Q1*H*Z1**H, B = Q1*T*Z1**H, *> *> as computed by ZGGHRD. *> *> If JOB='S', then the Hessenberg-triangular pair (H,T) is *> also reduced to generalized Schur form, *> *> H = Q*S*Z**H, T = Q*P*Z**H, *> *> where Q and Z are unitary matrices and S and P are upper triangular. *> *> Optionally, the unitary matrix Q from the generalized Schur *> factorization may be postmultiplied into an input matrix Q1, and the *> unitary matrix Z may be postmultiplied into an input matrix Z1. *> If Q1 and Z1 are the unitary matrices from ZGGHRD that reduced *> the matrix pair (A,B) to generalized Hessenberg form, then the output *> matrices Q1*Q and Z1*Z are the unitary factors from the generalized *> Schur factorization of (A,B): *> *> A = (Q1*Q)*S*(Z1*Z)**H, B = (Q1*Q)*P*(Z1*Z)**H. *> *> To avoid overflow, eigenvalues of the matrix pair (H,T) *> (equivalently, of (A,B)) are computed as a pair of complex values *> (alpha,beta). If beta is nonzero, lambda = alpha / beta is an *> eigenvalue of the generalized nonsymmetric eigenvalue problem (GNEP) *> A*x = lambda*B*x *> and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the *> alternate form of the GNEP *> mu*A*y = B*y. *> The values of alpha and beta for the i-th eigenvalue can be read *> directly from the generalized Schur form: alpha = S(i,i), *> beta = P(i,i). *> *> Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix *> Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973), *> pp. 241--256. *> \endverbatim * * Arguments: * ========== * *> \param[in] JOB *> \verbatim *> JOB is CHARACTER*1 *> = 'E': Compute eigenvalues only; *> = 'S': Computer eigenvalues and the Schur form. *> \endverbatim *> *> \param[in] COMPQ *> \verbatim *> COMPQ is CHARACTER*1 *> = 'N': Left Schur vectors (Q) are not computed; *> = 'I': Q is initialized to the unit matrix and the matrix Q *> of left Schur vectors of (H,T) is returned; *> = 'V': Q must contain a unitary matrix Q1 on entry and *> the product Q1*Q is returned. *> \endverbatim *> *> \param[in] COMPZ *> \verbatim *> COMPZ is CHARACTER*1 *> = 'N': Right Schur vectors (Z) are not computed; *> = 'I': Q is initialized to the unit matrix and the matrix Z *> of right Schur vectors of (H,T) is returned; *> = 'V': Z must contain a unitary matrix Z1 on entry and *> the product Z1*Z is returned. *> \endverbatim *> *> \param[in] N *> \verbatim *> N is INTEGER *> The order of the matrices H, T, Q, and Z. N >= 0. *> \endverbatim *> *> \param[in] ILO *> \verbatim *> ILO is INTEGER *> \endverbatim *> *> \param[in] IHI *> \verbatim *> IHI is INTEGER *> ILO and IHI mark the rows and columns of H which are in *> Hessenberg form. It is assumed that A is already upper *> triangular in rows and columns 1:ILO-1 and IHI+1:N. *> If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0. *> \endverbatim *> *> \param[in,out] H *> \verbatim *> H is COMPLEX*16 array, dimension (LDH, N) *> On entry, the N-by-N upper Hessenberg matrix H. *> On exit, if JOB = 'S', H contains the upper triangular *> matrix S from the generalized Schur factorization. *> If JOB = 'E', the diagonal of H matches that of S, but *> the rest of H is unspecified. *> \endverbatim *> *> \param[in] LDH *> \verbatim *> LDH is INTEGER *> The leading dimension of the array H. LDH >= max( 1, N ). *> \endverbatim *> *> \param[in,out] T *> \verbatim *> T is COMPLEX*16 array, dimension (LDT, N) *> On entry, the N-by-N upper triangular matrix T. *> On exit, if JOB = 'S', T contains the upper triangular *> matrix P from the generalized Schur factorization. *> If JOB = 'E', the diagonal of T matches that of P, but *> the rest of T is unspecified. *> \endverbatim *> *> \param[in] LDT *> \verbatim *> LDT is INTEGER *> The leading dimension of the array T. LDT >= 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(i) = S(i,i) in the generalized Schur *> factorization. *> \endverbatim *> *> \param[out] BETA *> \verbatim *> BETA is COMPLEX*16 array, dimension (N) *> The real non-negative scalars beta that define the *> eigenvalues of GNEP. BETA(i) = P(i,i) in the generalized *> Schur factorization. *> *> 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[in,out] Q *> \verbatim *> Q is COMPLEX*16 array, dimension (LDQ, N) *> On entry, if COMPQ = 'V', the unitary matrix Q1 used in the *> reduction of (A,B) to generalized Hessenberg form. *> On exit, if COMPQ = 'I', the unitary matrix of left Schur *> vectors of (H,T), and if COMPQ = 'V', the unitary matrix of *> left Schur vectors of (A,B). *> Not referenced if COMPQ = 'N'. *> \endverbatim *> *> \param[in] LDQ *> \verbatim *> LDQ is INTEGER *> The leading dimension of the array Q. LDQ >= 1. *> If COMPQ='V' or 'I', then LDQ >= N. *> \endverbatim *> *> \param[in,out] Z *> \verbatim *> Z is COMPLEX*16 array, dimension (LDZ, N) *> On entry, if COMPZ = 'V', the unitary matrix Z1 used in the *> reduction of (A,B) to generalized Hessenberg form. *> On exit, if COMPZ = 'I', the unitary matrix of right Schur *> vectors of (H,T), and if COMPZ = 'V', the unitary matrix of *> right Schur vectors of (A,B). *> Not referenced if COMPZ = 'N'. *> \endverbatim *> *> \param[in] LDZ *> \verbatim *> LDZ is INTEGER *> The leading dimension of the array Z. LDZ >= 1. *> If COMPZ='V' or 'I', then LDZ >= 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,N). *> *> 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 (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 did not converge. (H,T) is not *> in Schur form, but ALPHA(i) and BETA(i), *> i=INFO+1,...,N should be correct. *> = N+1,...,2*N: the shift calculation failed. (H,T) is not *> in Schur form, but ALPHA(i) and BETA(i), *> i=INFO-N+1,...,N should be correct. *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date April 2012 * *> \ingroup complex16GEcomputational * *> \par Further Details: * ===================== *> *> \verbatim *> *> We assume that complex ABS works as long as its value is less than *> overflow. *> \endverbatim *> * ===================================================================== SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, \$ ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, \$ RWORK, INFO ) * * -- LAPACK computational routine (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * April 2012 * * .. Scalar Arguments .. CHARACTER COMPQ, COMPZ, JOB INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N * .. * .. Array Arguments .. DOUBLE PRECISION RWORK( * ) COMPLEX*16 ALPHA( * ), BETA( * ), H( LDH, * ), \$ Q( LDQ, * ), T( LDT, * ), WORK( * ), \$ Z( LDZ, * ) * .. * * ===================================================================== * * .. Parameters .. COMPLEX*16 CZERO, CONE PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), \$ CONE = ( 1.0D+0, 0.0D+0 ) ) DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) DOUBLE PRECISION HALF PARAMETER ( HALF = 0.5D+0 ) * .. * .. Local Scalars .. LOGICAL ILAZR2, ILAZRO, ILQ, ILSCHR, ILZ, LQUERY INTEGER ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST, \$ ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER, \$ JR, MAXIT DOUBLE PRECISION ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL, \$ C, SAFMIN, TEMP, TEMP2, TEMPR, ULP COMPLEX*16 ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2, \$ CTEMP3, ESHIFT, RTDISC, S, SHIFT, SIGNBC, T1, \$ U12, X * .. * .. External Functions .. LOGICAL LSAME DOUBLE PRECISION DLAMCH, ZLANHS EXTERNAL LSAME, DLAMCH, ZLANHS * .. * .. External Subroutines .. EXTERNAL XERBLA, ZLARTG, ZLASET, ZROT, ZSCAL * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN, \$ SQRT * .. * .. Statement Functions .. DOUBLE PRECISION ABS1 * .. * .. Statement Function definitions .. ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) ) * .. * .. Executable Statements .. * * Decode JOB, COMPQ, COMPZ * IF( LSAME( JOB, 'E' ) ) THEN ILSCHR = .FALSE. ISCHUR = 1 ELSE IF( LSAME( JOB, 'S' ) ) THEN ILSCHR = .TRUE. ISCHUR = 2 ELSE ISCHUR = 0 END IF * IF( LSAME( COMPQ, 'N' ) ) THEN ILQ = .FALSE. ICOMPQ = 1 ELSE IF( LSAME( COMPQ, 'V' ) ) THEN ILQ = .TRUE. ICOMPQ = 2 ELSE IF( LSAME( COMPQ, 'I' ) ) THEN ILQ = .TRUE. ICOMPQ = 3 ELSE ICOMPQ = 0 END IF * IF( LSAME( COMPZ, 'N' ) ) THEN ILZ = .FALSE. ICOMPZ = 1 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN ILZ = .TRUE. ICOMPZ = 2 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN ILZ = .TRUE. ICOMPZ = 3 ELSE ICOMPZ = 0 END IF * * Check Argument Values * INFO = 0 WORK( 1 ) = MAX( 1, N ) LQUERY = ( LWORK.EQ.-1 ) IF( ISCHUR.EQ.0 ) THEN INFO = -1 ELSE IF( ICOMPQ.EQ.0 ) THEN INFO = -2 ELSE IF( ICOMPZ.EQ.0 ) THEN INFO = -3 ELSE IF( N.LT.0 ) THEN INFO = -4 ELSE IF( ILO.LT.1 ) THEN INFO = -5 ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN INFO = -6 ELSE IF( LDH.LT.N ) THEN INFO = -8 ELSE IF( LDT.LT.N ) THEN INFO = -10 ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN INFO = -14 ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN INFO = -16 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN INFO = -18 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'ZHGEQZ', -INFO ) RETURN ELSE IF( LQUERY ) THEN RETURN END IF * * Quick return if possible * * WORK( 1 ) = CMPLX( 1 ) IF( N.LE.0 ) THEN WORK( 1 ) = DCMPLX( 1 ) RETURN END IF * * Initialize Q and Z * IF( ICOMPQ.EQ.3 ) \$ CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ ) IF( ICOMPZ.EQ.3 ) \$ CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ ) * * Machine Constants * IN = IHI + 1 - ILO SAFMIN = DLAMCH( 'S' ) ULP = DLAMCH( 'E' )*DLAMCH( 'B' ) ANORM = ZLANHS( 'F', IN, H( ILO, ILO ), LDH, RWORK ) BNORM = ZLANHS( 'F', IN, T( ILO, ILO ), LDT, RWORK ) ATOL = MAX( SAFMIN, ULP*ANORM ) BTOL = MAX( SAFMIN, ULP*BNORM ) ASCALE = ONE / MAX( SAFMIN, ANORM ) BSCALE = ONE / MAX( SAFMIN, BNORM ) * * * Set Eigenvalues IHI+1:N * DO 10 J = IHI + 1, N ABSB = ABS( T( J, J ) ) IF( ABSB.GT.SAFMIN ) THEN SIGNBC = DCONJG( T( J, J ) / ABSB ) T( J, J ) = ABSB IF( ILSCHR ) THEN CALL ZSCAL( J-1, SIGNBC, T( 1, J ), 1 ) CALL ZSCAL( J, SIGNBC, H( 1, J ), 1 ) ELSE CALL ZSCAL( 1, SIGNBC, H( J, J ), 1 ) END IF IF( ILZ ) \$ CALL ZSCAL( N, SIGNBC, Z( 1, J ), 1 ) ELSE T( J, J ) = CZERO END IF ALPHA( J ) = H( J, J ) BETA( J ) = T( J, J ) 10 CONTINUE * * If IHI < ILO, skip QZ steps * IF( IHI.LT.ILO ) \$ GO TO 190 * * MAIN QZ ITERATION LOOP * * Initialize dynamic indices * * Eigenvalues ILAST+1:N have been found. * Column operations modify rows IFRSTM:whatever * Row operations modify columns whatever:ILASTM * * If only eigenvalues are being computed, then * IFRSTM is the row of the last splitting row above row ILAST; * this is always at least ILO. * IITER counts iterations since the last eigenvalue was found, * to tell when to use an extraordinary shift. * MAXIT is the maximum number of QZ sweeps allowed. * ILAST = IHI IF( ILSCHR ) THEN IFRSTM = 1 ILASTM = N ELSE IFRSTM = ILO ILASTM = IHI END IF IITER = 0 ESHIFT = CZERO MAXIT = 30*( IHI-ILO+1 ) * DO 170 JITER = 1, MAXIT * * Check for too many iterations. * IF( JITER.GT.MAXIT ) \$ GO TO 180 * * Split the matrix if possible. * * Two tests: * 1: H(j,j-1)=0 or j=ILO * 2: T(j,j)=0 * * Special case: j=ILAST * IF( ILAST.EQ.ILO ) THEN GO TO 60 ELSE IF( ABS1( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN H( ILAST, ILAST-1 ) = CZERO GO TO 60 END IF END IF * IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN T( ILAST, ILAST ) = CZERO GO TO 50 END IF * * General case: j