phynexis-deps/lapack
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
Files
9aed94de6d78beddb218baaebce809193b8160fc
lapack/SRC/chgeqz.f
T
NAKATA Maho ac82605ea9 refactor: rename ABS1 statement function to CABS1 for consistency 
Several LAPACK, BLAS, and CBLAS source files defined a local statement
function named ABS1 for the complex 1-norm approximation:

    ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
    ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )

The majority of the codebase already uses CABS1 for this identical
purpose. This commit renames ABS1 to CABS1 in all remaining files
(definition line, declaration line, and all call sites within the
same file) to make the naming consistent across the repository.

A small number of fixed-form lines required continuation-line splits
to stay within the 72-column limit after the rename.

No numerical change. Statement functions are file-local in Fortran,
so there is no ABI or interface impact.

This is a preparatory cleanup before inlining these statement
functions (see issue #1200).
2026-03-28 14:46:42 +09:00

900 lines
29 KiB
FortranFixed
Raw Blame History

*> \brief \b CHGEQZ
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> Download CHGEQZ + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chgeqz.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chgeqz.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chgeqz.f">
*> [TXT]</a>
*
* Definition:
* ===========
*
* SUBROUTINE CHGEQZ( 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 ..
* REAL RWORK( * )
* COMPLEX ALPHA( * ), BETA( * ), H( LDH, * ),
* $ Q( LDQ, * ), T( LDT, * ), WORK( * ),
* $ Z( LDZ, * )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> CHGEQZ 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 CGGHRD.
*>
*> 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 CGGHRD 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 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 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 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 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 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 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 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 REAL 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.
*
*> \ingroup hgeqz
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> We assume that complex ABS works as long as its value is less than
*> overflow.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T,
$ LDT,
$ ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK,
$ RWORK, INFO )
IMPLICIT NONE
*
* -- LAPACK computational 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 COMPQ, COMPZ, JOB
INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
* ..
* .. Array Arguments ..
REAL RWORK( * )
COMPLEX ALPHA( * ), BETA( * ), H( LDH, * ),
$ Q( LDQ, * ), T( LDT, * ), WORK( * ),
$ Z( LDZ, * )
* ..
*
* =====================================================================
*
* .. Parameters ..
COMPLEX CZERO, CONE
PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
$ CONE = ( 1.0E+0, 0.0E+0 ) )
REAL ZERO, ONE
PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
REAL HALF
PARAMETER ( HALF = 0.5E+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
REAL ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL,
$ C, SAFMIN, TEMP, TEMP2, TEMPR, ULP
COMPLEX ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2,
$ CTEMP3, ESHIFT, S, SHIFT, SIGNBC,
$ U12, X, ABI12, Y
* ..
* .. External Functions ..
COMPLEX CLADIV
LOGICAL LSAME
REAL CLANHS, SLAMCH
EXTERNAL CLADIV, LSAME, CLANHS, SLAMCH
* ..
* .. External Subroutines ..
EXTERNAL CLARTG, CLASET, CROT, CSCAL, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL, SQRT
* ..
* .. Statement Functions ..
REAL CABS1
* ..
* .. Statement Function definitions ..
CABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( 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
ILSCHR = .TRUE.
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
ILQ = .TRUE.
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
ILZ = .TRUE.
ICOMPZ = 0
END IF
*
* Check Argument Values
*
INFO = 0
WORK( 1 ) = CMPLX( 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( 'CHGEQZ', -INFO )
RETURN
ELSE IF( LQUERY ) THEN
RETURN
END IF
*
* Quick return if possible
*
* WORK( 1 ) = CMPLX( 1 )
IF( N.LE.0 ) THEN
WORK( 1 ) = CMPLX( 1 )
RETURN
END IF
*
* Initialize Q and Z
*
IF( ICOMPQ.EQ.3 )
$ CALL CLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
IF( ICOMPZ.EQ.3 )
$ CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
*
* Machine Constants
*
IN = IHI + 1 - ILO
SAFMIN = SLAMCH( 'S' )
ULP = SLAMCH( 'E' )*SLAMCH( 'B' )
ANORM = CLANHS( 'F', IN, H( ILO, ILO ), LDH, RWORK )
BNORM = CLANHS( '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 = CONJG( T( J, J ) / ABSB )
T( J, J ) = ABSB
IF( ILSCHR ) THEN
CALL CSCAL( J-1, SIGNBC, T( 1, J ), 1 )
CALL CSCAL( J, SIGNBC, H( 1, J ), 1 )
ELSE
CALL CSCAL( 1, SIGNBC, H( J, J ), 1 )
END IF
IF( ILZ )
$ CALL CSCAL( 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( CABS1( H( ILAST, ILAST-1 ) ).LE.MAX( SAFMIN,
$ ULP*( CABS1( H( ILAST, ILAST ) )
$ + CABS1( H( ILAST-1, ILAST-1 ) ) ) ) ) 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<ILAST
*
DO 40 J = ILAST - 1, ILO, -1
*
* Test 1: for H(j,j-1)=0 or j=ILO
*
IF( J.EQ.ILO ) THEN
ILAZRO = .TRUE.
ELSE
IF( CABS1( H( J, J-1 ) ).LE.MAX( SAFMIN, ULP*(
$ CABS1( H( J, J ) ) + CABS1( H( J-1, J-1 ) )
$ ) ) ) THEN
H( J, J-1 ) = CZERO
ILAZRO = .TRUE.
ELSE
ILAZRO = .FALSE.
END IF
END IF
*
* Test 2: for T(j,j)=0
*
IF( ABS( T( J, J ) ).LT.BTOL ) THEN
T( J, J ) = CZERO
*
* Test 1a: Check for 2 consecutive small subdiagonals in A
*
ILAZR2 = .FALSE.
IF( .NOT.ILAZRO ) THEN
IF( CABS1( H( J, J-1 ) )*( ASCALE*CABS1( H( J+1,
$ J ) ) ).LE.CABS1( H( J, J ) )*( ASCALE*ATOL ) )
$ ILAZR2 = .TRUE.
END IF
*
* If both tests pass (1 & 2), i.e., the leading diagonal
* element of B in the block is zero, split a 1x1 block off
* at the top. (I.e., at the J-th row/column) The leading
* diagonal element of the remainder can also be zero, so
* this may have to be done repeatedly.
*
IF( ILAZRO .OR. ILAZR2 ) THEN
DO 20 JCH = J, ILAST - 1
CTEMP = H( JCH, JCH )
CALL CLARTG( CTEMP, H( JCH+1, JCH ), C, S,
$ H( JCH, JCH ) )
H( JCH+1, JCH ) = CZERO
CALL CROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH,
$ H( JCH+1, JCH+1 ), LDH, C, S )
CALL CROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT,
$ T( JCH+1, JCH+1 ), LDT, C, S )
IF( ILQ )
$ CALL CROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ),
$ 1,
$ C, CONJG( S ) )
IF( ILAZR2 )
$ H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C
ILAZR2 = .FALSE.
IF( CABS1( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN
IF( JCH+1.GE.ILAST ) THEN
GO TO 60
ELSE
IFIRST = JCH + 1
GO TO 70
END IF
END IF
T( JCH+1, JCH+1 ) = CZERO
20 CONTINUE
GO TO 50
ELSE
*
* Only test 2 passed -- chase the zero to T(ILAST,ILAST)
* Then process as in the case T(ILAST,ILAST)=0
*
DO 30 JCH = J, ILAST - 1
CTEMP = T( JCH, JCH+1 )
CALL CLARTG( CTEMP, T( JCH+1, JCH+1 ), C, S,
$ T( JCH, JCH+1 ) )
T( JCH+1, JCH+1 ) = CZERO
IF( JCH.LT.ILASTM-1 )
$ CALL CROT( ILASTM-JCH-1, T( JCH, JCH+2 ),
$ LDT,
$ T( JCH+1, JCH+2 ), LDT, C, S )
CALL CROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH,
$ H( JCH+1, JCH-1 ), LDH, C, S )
IF( ILQ )
$ CALL CROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ),
$ 1,
$ C, CONJG( S ) )
CTEMP = H( JCH+1, JCH )
CALL CLARTG( CTEMP, H( JCH+1, JCH-1 ), C, S,
$ H( JCH+1, JCH ) )
H( JCH+1, JCH-1 ) = CZERO
CALL CROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1,
$ H( IFRSTM, JCH-1 ), 1, C, S )
CALL CROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1,
$ T( IFRSTM, JCH-1 ), 1, C, S )
IF( ILZ )
$ CALL CROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ),
$ 1,
$ C, S )
30 CONTINUE
GO TO 50
END IF
ELSE IF( ILAZRO ) THEN
*
* Only test 1 passed -- work on J:ILAST
*
IFIRST = J
GO TO 70
END IF
*
* Neither test passed -- try next J
*
40 CONTINUE
*
* (Drop-through is "impossible")
*
INFO = 2*N + 1
GO TO 210
*
* T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
* 1x1 block.
*
50 CONTINUE
CTEMP = H( ILAST, ILAST )
CALL CLARTG( CTEMP, H( ILAST, ILAST-1 ), C, S,
$ H( ILAST, ILAST ) )
H( ILAST, ILAST-1 ) = CZERO
CALL CROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1,
$ H( IFRSTM, ILAST-1 ), 1, C, S )
CALL CROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1,
$ T( IFRSTM, ILAST-1 ), 1, C, S )
IF( ILZ )
$ CALL CROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C,
$ S )
*
* H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHA and BETA
*
60 CONTINUE
ABSB = ABS( T( ILAST, ILAST ) )
IF( ABSB.GT.SAFMIN ) THEN
SIGNBC = CONJG( T( ILAST, ILAST ) / ABSB )
T( ILAST, ILAST ) = ABSB
IF( ILSCHR ) THEN
CALL CSCAL( ILAST-IFRSTM, SIGNBC, T( IFRSTM, ILAST ),
$ 1 )
CALL CSCAL( ILAST+1-IFRSTM, SIGNBC, H( IFRSTM,
$ ILAST ),
$ 1 )
ELSE
CALL CSCAL( 1, SIGNBC, H( ILAST, ILAST ), 1 )
END IF
IF( ILZ )
$ CALL CSCAL( N, SIGNBC, Z( 1, ILAST ), 1 )
ELSE
T( ILAST, ILAST ) = CZERO
END IF
ALPHA( ILAST ) = H( ILAST, ILAST )
BETA( ILAST ) = T( ILAST, ILAST )
*
* Go to next block -- exit if finished.
*
ILAST = ILAST - 1
IF( ILAST.LT.ILO )
$ GO TO 190
*
* Reset counters
*
IITER = 0
ESHIFT = CZERO
IF( .NOT.ILSCHR ) THEN
ILASTM = ILAST
IF( IFRSTM.GT.ILAST )
$ IFRSTM = ILO
END IF
GO TO 160
*
* QZ step
*
* This iteration only involves rows/columns IFIRST:ILAST. We
* assume IFIRST < ILAST, and that the diagonal of B is non-zero.
*
70 CONTINUE
IITER = IITER + 1
IF( .NOT.ILSCHR ) THEN
IFRSTM = IFIRST
END IF
*
* Compute the Shift.
*
* At this point, IFIRST < ILAST, and the diagonal elements of
* T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
* magnitude)
*
IF( ( IITER / 10 )*10.NE.IITER ) THEN
*
* The Wilkinson shift (AEP p.512), i.e., the eigenvalue of
* the bottom-right 2x2 block of A inv(B) which is nearest to
* the bottom-right element.
*
* We factor B as U*D, where U has unit diagonals, and
* compute (A*inv(D))*inv(U).
*
U12 = ( BSCALE*T( ILAST-1, ILAST ) ) /
$ ( BSCALE*T( ILAST, ILAST ) )
AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) /
$ ( BSCALE*T( ILAST-1, ILAST-1 ) )
AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) /
$ ( BSCALE*T( ILAST-1, ILAST-1 ) )
AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) /
$ ( BSCALE*T( ILAST, ILAST ) )
AD22 = ( ASCALE*H( ILAST, ILAST ) ) /
$ ( BSCALE*T( ILAST, ILAST ) )
ABI22 = AD22 - U12*AD21
ABI12 = AD12 - U12*AD11
*
SHIFT = ABI22
CTEMP = SQRT( ABI12 )*SQRT( AD21 )
TEMP = CABS1( CTEMP )
IF( CTEMP.NE.ZERO ) THEN
X = HALF*( AD11-SHIFT )
TEMP2 = CABS1( X )
TEMP = MAX( TEMP, CABS1( X ) )
Y = TEMP*SQRT( ( X / TEMP )**2+( CTEMP / TEMP )**2 )
IF( TEMP2.GT.ZERO ) THEN
IF( REAL( X / TEMP2 )*REAL( Y )+
$ AIMAG( X / TEMP2 )*AIMAG( Y ).LT.ZERO )Y = -Y
END IF
SHIFT = SHIFT - CTEMP*CLADIV( CTEMP, ( X+Y ) )
END IF
ELSE
*
* Exceptional shift. Chosen for no particularly good reason.
*
IF( ( IITER / 20 )*20.EQ.IITER .AND.
$ BSCALE*CABS1(T( ILAST, ILAST )).GT.SAFMIN ) THEN
ESHIFT = ESHIFT + ( ASCALE*H( ILAST,
$ ILAST ) )/( BSCALE*T( ILAST, ILAST ) )
ELSE
ESHIFT = ESHIFT + ( ASCALE*H( ILAST,
$ ILAST-1 ) )/( BSCALE*T( ILAST-1, ILAST-1 ) )
END IF
SHIFT = ESHIFT
END IF
*
* Now check for two consecutive small subdiagonals.
*
DO 80 J = ILAST - 1, IFIRST + 1, -1
ISTART = J
CTEMP = ASCALE*H( J, J ) - SHIFT*( BSCALE*T( J, J ) )
TEMP = CABS1( CTEMP )
TEMP2 = ASCALE*CABS1( H( J+1, J ) )
TEMPR = MAX( TEMP, TEMP2 )
IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
TEMP = TEMP / TEMPR
TEMP2 = TEMP2 / TEMPR
END IF
IF( CABS1( H( J, J-1 ) )*TEMP2.LE.TEMP*ATOL )
$ GO TO 90
80 CONTINUE
*
ISTART = IFIRST
CTEMP = ASCALE*H( IFIRST, IFIRST ) -
$ SHIFT*( BSCALE*T( IFIRST, IFIRST ) )
90 CONTINUE
*
* Do an implicit-shift QZ sweep.
*
* Initial Q
*
CTEMP2 = ASCALE*H( ISTART+1, ISTART )
CALL CLARTG( CTEMP, CTEMP2, C, S, CTEMP3 )
*
* Sweep
*
DO 150 J = ISTART, ILAST - 1
IF( J.GT.ISTART ) THEN
CTEMP = H( J, J-1 )
CALL CLARTG( CTEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) )
H( J+1, J-1 ) = CZERO
END IF
*
DO 100 JC = J, ILASTM
CTEMP = C*H( J, JC ) + S*H( J+1, JC )
H( J+1, JC ) = -CONJG( S )*H( J, JC ) + C*H( J+1, JC )
H( J, JC ) = CTEMP
CTEMP2 = C*T( J, JC ) + S*T( J+1, JC )
T( J+1, JC ) = -CONJG( S )*T( J, JC ) + C*T( J+1, JC )
T( J, JC ) = CTEMP2
100 CONTINUE
IF( ILQ ) THEN
DO 110 JR = 1, N
CTEMP = C*Q( JR, J ) + CONJG( S )*Q( JR, J+1 )
Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
Q( JR, J ) = CTEMP
110 CONTINUE
END IF
*
CTEMP = T( J+1, J+1 )
CALL CLARTG( CTEMP, T( J+1, J ), C, S, T( J+1, J+1 ) )
T( J+1, J ) = CZERO
*
DO 120 JR = IFRSTM, MIN( J+2, ILAST )
CTEMP = C*H( JR, J+1 ) + S*H( JR, J )
H( JR, J ) = -CONJG( S )*H( JR, J+1 ) + C*H( JR, J )
H( JR, J+1 ) = CTEMP
120 CONTINUE
DO 130 JR = IFRSTM, J
CTEMP = C*T( JR, J+1 ) + S*T( JR, J )
T( JR, J ) = -CONJG( S )*T( JR, J+1 ) + C*T( JR, J )
T( JR, J+1 ) = CTEMP
130 CONTINUE
IF( ILZ ) THEN
DO 140 JR = 1, N
CTEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
Z( JR, J ) = -CONJG( S )*Z( JR, J+1 ) + C*Z( JR, J )
Z( JR, J+1 ) = CTEMP
140 CONTINUE
END IF
150 CONTINUE
*
160 CONTINUE
*
170 CONTINUE
*
* Drop-through = non-convergence
*
180 CONTINUE
INFO = ILAST
GO TO 210
*
* Successful completion of all QZ steps
*
190 CONTINUE
*
* Set Eigenvalues 1:ILO-1
*
DO 200 J = 1, ILO - 1
ABSB = ABS( T( J, J ) )
IF( ABSB.GT.SAFMIN ) THEN
SIGNBC = CONJG( T( J, J ) / ABSB )
T( J, J ) = ABSB
IF( ILSCHR ) THEN
CALL CSCAL( J-1, SIGNBC, T( 1, J ), 1 )
CALL CSCAL( J, SIGNBC, H( 1, J ), 1 )
ELSE
CALL CSCAL( 1, SIGNBC, H( J, J ), 1 )
END IF
IF( ILZ )
$ CALL CSCAL( N, SIGNBC, Z( 1, J ), 1 )
ELSE
T( J, J ) = CZERO
END IF
ALPHA( J ) = H( J, J )
BETA( J ) = T( J, J )
200 CONTINUE
*
* Normal Termination
*
INFO = 0
*
* Exit (other than argument error) -- return optimal workspace size
*
210 CONTINUE
WORK( 1 ) = CMPLX( N )
RETURN
*
* End of CHGEQZ
*
END
Reference in New Issue View Git Blame Copy Permalink