diff options
Diffstat (limited to 'src/lib/lapack/zhgeqz.f')
| -rw-r--r-- | src/lib/lapack/zhgeqz.f | 759 |
1 files changed, 0 insertions, 759 deletions
diff --git a/src/lib/lapack/zhgeqz.f b/src/lib/lapack/zhgeqz.f deleted file mode 100644 index 6a9403bd..00000000 --- a/src/lib/lapack/zhgeqz.f +++ /dev/null @@ -1,759 +0,0 @@ - SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, - $ ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, - $ RWORK, INFO ) -* -* -- LAPACK routine (version 3.1) -- -* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. -* November 2006 -* -* .. 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, * ) -* .. -* -* Purpose -* ======= -* -* 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. -* -* Arguments -* ========= -* -* JOB (input) CHARACTER*1 -* = 'E': Compute eigenvalues only; -* = 'S': Computer eigenvalues and the Schur form. -* -* COMPQ (input) 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. -* -* COMPZ (input) 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. -* -* N (input) INTEGER -* The order of the matrices H, T, Q, and Z. N >= 0. -* -* ILO (input) INTEGER -* IHI (input) 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. -* -* H (input/output) 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. -* -* LDH (input) INTEGER -* The leading dimension of the array H. LDH >= max( 1, N ). -* -* T (input/output) 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. -* -* LDT (input) INTEGER -* The leading dimension of the array T. LDT >= max( 1, N ). -* -* ALPHA (output) 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. -* -* BETA (output) 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. -* -* Q (input/output) COMPLEX*16 array, dimension (LDQ, N) -* On entry, if COMPZ = 'V', the unitary matrix Q1 used in the -* reduction of (A,B) to generalized Hessenberg form. -* On exit, if COMPZ = 'I', the unitary matrix of left Schur -* vectors of (H,T), and if COMPZ = 'V', the unitary matrix of -* left Schur vectors of (A,B). -* Not referenced if COMPZ = 'N'. -* -* LDQ (input) INTEGER -* The leading dimension of the array Q. LDQ >= 1. -* If COMPQ='V' or 'I', then LDQ >= N. -* -* Z (input/output) 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'. -* -* LDZ (input) INTEGER -* The leading dimension of the array Z. LDZ >= 1. -* If COMPZ='V' or 'I', then LDZ >= N. -* -* WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) -* On exit, if INFO >= 0, WORK(1) returns the optimal LWORK. -* -* LWORK (input) 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. -* -* RWORK (workspace) DOUBLE PRECISION array, dimension (N) -* -* INFO (output) 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. -* -* Further Details -* =============== -* -* We assume that complex ABS works as long as its value is less than -* overflow. -* -* ===================================================================== -* -* .. 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 - H( J, J ) = H( J, J )*SIGNBC - 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<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( ABS1( H( J, J-1 ) ).LE.ATOL ) 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( ABS1( H( J, J-1 ) )*( ASCALE*ABS1( H( J+1, - $ J ) ) ).LE.ABS1( 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 ZLARTG( CTEMP, H( JCH+1, JCH ), C, S, - $ H( JCH, JCH ) ) - H( JCH+1, JCH ) = CZERO - CALL ZROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH, - $ H( JCH+1, JCH+1 ), LDH, C, S ) - CALL ZROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT, - $ T( JCH+1, JCH+1 ), LDT, C, S ) - IF( ILQ ) - $ CALL ZROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1, - $ C, DCONJG( S ) ) - IF( ILAZR2 ) - $ H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C - ILAZR2 = .FALSE. - IF( ABS1( 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 ZLARTG( CTEMP, T( JCH+1, JCH+1 ), C, S, - $ T( JCH, JCH+1 ) ) - T( JCH+1, JCH+1 ) = CZERO - IF( JCH.LT.ILASTM-1 ) - $ CALL ZROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT, - $ T( JCH+1, JCH+2 ), LDT, C, S ) - CALL ZROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH, - $ H( JCH+1, JCH-1 ), LDH, C, S ) - IF( ILQ ) - $ CALL ZROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1, - $ C, DCONJG( S ) ) - CTEMP = H( JCH+1, JCH ) - CALL ZLARTG( CTEMP, H( JCH+1, JCH-1 ), C, S, - $ H( JCH+1, JCH ) ) - H( JCH+1, JCH-1 ) = CZERO - CALL ZROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1, - $ H( IFRSTM, JCH-1 ), 1, C, S ) - CALL ZROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1, - $ T( IFRSTM, JCH-1 ), 1, C, S ) - IF( ILZ ) - $ CALL ZROT( 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 ZLARTG( CTEMP, H( ILAST, ILAST-1 ), C, S, - $ H( ILAST, ILAST ) ) - H( ILAST, ILAST-1 ) = CZERO - CALL ZROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1, - $ H( IFRSTM, ILAST-1 ), 1, C, S ) - CALL ZROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1, - $ T( IFRSTM, ILAST-1 ), 1, C, S ) - IF( ILZ ) - $ CALL ZROT( 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 = DCONJG( T( ILAST, ILAST ) / ABSB ) - T( ILAST, ILAST ) = ABSB - IF( ILSCHR ) THEN - CALL ZSCAL( ILAST-IFRSTM, SIGNBC, T( IFRSTM, ILAST ), 1 ) - CALL ZSCAL( ILAST+1-IFRSTM, SIGNBC, H( IFRSTM, ILAST ), - $ 1 ) - ELSE - H( ILAST, ILAST ) = H( ILAST, ILAST )*SIGNBC - END IF - IF( ILZ ) - $ CALL ZSCAL( 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 -* - T1 = HALF*( AD11+ABI22 ) - RTDISC = SQRT( T1**2+AD12*AD21-AD11*AD22 ) - TEMP = DBLE( T1-ABI22 )*DBLE( RTDISC ) + - $ DIMAG( T1-ABI22 )*DIMAG( RTDISC ) - IF( TEMP.LE.ZERO ) THEN - SHIFT = T1 + RTDISC - ELSE - SHIFT = T1 - RTDISC - END IF - ELSE -* -* Exceptional shift. Chosen for no particularly good reason. -* - ESHIFT = ESHIFT + DCONJG( ( ASCALE*H( ILAST-1, ILAST ) ) / - $ ( BSCALE*T( ILAST-1, ILAST-1 ) ) ) - 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 = ABS1( CTEMP ) - TEMP2 = ASCALE*ABS1( 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( ABS1( 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 ZLARTG( CTEMP, CTEMP2, C, S, CTEMP3 ) -* -* Sweep -* - DO 150 J = ISTART, ILAST - 1 - IF( J.GT.ISTART ) THEN - CTEMP = H( J, J-1 ) - CALL ZLARTG( 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 ) = -DCONJG( 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 ) = -DCONJG( 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 ) + DCONJG( 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 ZLARTG( 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 ) = -DCONJG( 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 ) = -DCONJG( 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 ) = -DCONJG( 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 = 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 - H( J, J ) = H( J, J )*SIGNBC - 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 ) - 200 CONTINUE -* -* Normal Termination -* - INFO = 0 -* -* Exit (other than argument error) -- return optimal workspace size -* - 210 CONTINUE - WORK( 1 ) = DCMPLX( N ) - RETURN -* -* End of ZHGEQZ -* - END |