 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, 00002 $ ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, 00003 $ RWORK, INFO ) 00004 * 00005 * -- LAPACK routine (version 3.3.1) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * -- April 2011 -- 00009 * 00010 * .. Scalar Arguments .. 00011 CHARACTER COMPQ, COMPZ, JOB 00012 INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N 00013 * .. 00014 * .. Array Arguments .. 00015 REAL RWORK( * ) 00016 COMPLEX ALPHA( * ), BETA( * ), H( LDH, * ), 00017 $ Q( LDQ, * ), T( LDT, * ), WORK( * ), 00018 $ Z( LDZ, * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * CHGEQZ computes the eigenvalues of a complex matrix pair (H,T), 00025 * where H is an upper Hessenberg matrix and T is upper triangular, 00026 * using the single-shift QZ method. 00027 * Matrix pairs of this type are produced by the reduction to 00028 * generalized upper Hessenberg form of a complex matrix pair (A,B): 00029 * 00030 * A = Q1*H*Z1**H, B = Q1*T*Z1**H, 00031 * 00032 * as computed by CGGHRD. 00033 * 00034 * If JOB='S', then the Hessenberg-triangular pair (H,T) is 00035 * also reduced to generalized Schur form, 00036 * 00037 * H = Q*S*Z**H, T = Q*P*Z**H, 00038 * 00039 * where Q and Z are unitary matrices and S and P are upper triangular. 00040 * 00041 * Optionally, the unitary matrix Q from the generalized Schur 00042 * factorization may be postmultiplied into an input matrix Q1, and the 00043 * unitary matrix Z may be postmultiplied into an input matrix Z1. 00044 * If Q1 and Z1 are the unitary matrices from CGGHRD that reduced 00045 * the matrix pair (A,B) to generalized Hessenberg form, then the output 00046 * matrices Q1*Q and Z1*Z are the unitary factors from the generalized 00047 * Schur factorization of (A,B): 00048 * 00049 * A = (Q1*Q)*S*(Z1*Z)**H, B = (Q1*Q)*P*(Z1*Z)**H. 00050 * 00051 * To avoid overflow, eigenvalues of the matrix pair (H,T) 00052 * (equivalently, of (A,B)) are computed as a pair of complex values 00053 * (alpha,beta). If beta is nonzero, lambda = alpha / beta is an 00054 * eigenvalue of the generalized nonsymmetric eigenvalue problem (GNEP) 00055 * A*x = lambda*B*x 00056 * and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the 00057 * alternate form of the GNEP 00058 * mu*A*y = B*y. 00059 * The values of alpha and beta for the i-th eigenvalue can be read 00060 * directly from the generalized Schur form: alpha = S(i,i), 00061 * beta = P(i,i). 00062 * 00063 * Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix 00064 * Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973), 00065 * pp. 241--256. 00066 * 00067 * Arguments 00068 * ========= 00069 * 00070 * JOB (input) CHARACTER*1 00071 * = 'E': Compute eigenvalues only; 00072 * = 'S': Computer eigenvalues and the Schur form. 00073 * 00074 * COMPQ (input) CHARACTER*1 00075 * = 'N': Left Schur vectors (Q) are not computed; 00076 * = 'I': Q is initialized to the unit matrix and the matrix Q 00077 * of left Schur vectors of (H,T) is returned; 00078 * = 'V': Q must contain a unitary matrix Q1 on entry and 00079 * the product Q1*Q is returned. 00080 * 00081 * COMPZ (input) CHARACTER*1 00082 * = 'N': Right Schur vectors (Z) are not computed; 00083 * = 'I': Q is initialized to the unit matrix and the matrix Z 00084 * of right Schur vectors of (H,T) is returned; 00085 * = 'V': Z must contain a unitary matrix Z1 on entry and 00086 * the product Z1*Z is returned. 00087 * 00088 * N (input) INTEGER 00089 * The order of the matrices H, T, Q, and Z. N >= 0. 00090 * 00091 * ILO (input) INTEGER 00092 * IHI (input) INTEGER 00093 * ILO and IHI mark the rows and columns of H which are in 00094 * Hessenberg form. It is assumed that A is already upper 00095 * triangular in rows and columns 1:ILO-1 and IHI+1:N. 00096 * If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0. 00097 * 00098 * H (input/output) COMPLEX array, dimension (LDH, N) 00099 * On entry, the N-by-N upper Hessenberg matrix H. 00100 * On exit, if JOB = 'S', H contains the upper triangular 00101 * matrix S from the generalized Schur factorization. 00102 * If JOB = 'E', the diagonal of H matches that of S, but 00103 * the rest of H is unspecified. 00104 * 00105 * LDH (input) INTEGER 00106 * The leading dimension of the array H. LDH >= max( 1, N ). 00107 * 00108 * T (input/output) COMPLEX array, dimension (LDT, N) 00109 * On entry, the N-by-N upper triangular matrix T. 00110 * On exit, if JOB = 'S', T contains the upper triangular 00111 * matrix P from the generalized Schur factorization. 00112 * If JOB = 'E', the diagonal of T matches that of P, but 00113 * the rest of T is unspecified. 00114 * 00115 * LDT (input) INTEGER 00116 * The leading dimension of the array T. LDT >= max( 1, N ). 00117 * 00118 * ALPHA (output) COMPLEX array, dimension (N) 00119 * The complex scalars alpha that define the eigenvalues of 00120 * GNEP. ALPHA(i) = S(i,i) in the generalized Schur 00121 * factorization. 00122 * 00123 * BETA (output) COMPLEX array, dimension (N) 00124 * The real non-negative scalars beta that define the 00125 * eigenvalues of GNEP. BETA(i) = P(i,i) in the generalized 00126 * Schur factorization. 00127 * 00128 * Together, the quantities alpha = ALPHA(j) and beta = BETA(j) 00129 * represent the j-th eigenvalue of the matrix pair (A,B), in 00130 * one of the forms lambda = alpha/beta or mu = beta/alpha. 00131 * Since either lambda or mu may overflow, they should not, 00132 * in general, be computed. 00133 * 00134 * Q (input/output) COMPLEX array, dimension (LDQ, N) 00135 * On entry, if COMPZ = 'V', the unitary matrix Q1 used in the 00136 * reduction of (A,B) to generalized Hessenberg form. 00137 * On exit, if COMPZ = 'I', the unitary matrix of left Schur 00138 * vectors of (H,T), and if COMPZ = 'V', the unitary matrix of 00139 * left Schur vectors of (A,B). 00140 * Not referenced if COMPZ = 'N'. 00141 * 00142 * LDQ (input) INTEGER 00143 * The leading dimension of the array Q. LDQ >= 1. 00144 * If COMPQ='V' or 'I', then LDQ >= N. 00145 * 00146 * Z (input/output) COMPLEX array, dimension (LDZ, N) 00147 * On entry, if COMPZ = 'V', the unitary matrix Z1 used in the 00148 * reduction of (A,B) to generalized Hessenberg form. 00149 * On exit, if COMPZ = 'I', the unitary matrix of right Schur 00150 * vectors of (H,T), and if COMPZ = 'V', the unitary matrix of 00151 * right Schur vectors of (A,B). 00152 * Not referenced if COMPZ = 'N'. 00153 * 00154 * LDZ (input) INTEGER 00155 * The leading dimension of the array Z. LDZ >= 1. 00156 * If COMPZ='V' or 'I', then LDZ >= N. 00157 * 00158 * WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) 00159 * On exit, if INFO >= 0, WORK(1) returns the optimal LWORK. 00160 * 00161 * LWORK (input) INTEGER 00162 * The dimension of the array WORK. LWORK >= max(1,N). 00163 * 00164 * If LWORK = -1, then a workspace query is assumed; the routine 00165 * only calculates the optimal size of the WORK array, returns 00166 * this value as the first entry of the WORK array, and no error 00167 * message related to LWORK is issued by XERBLA. 00168 * 00169 * RWORK (workspace) REAL array, dimension (N) 00170 * 00171 * INFO (output) INTEGER 00172 * = 0: successful exit 00173 * < 0: if INFO = -i, the i-th argument had an illegal value 00174 * = 1,...,N: the QZ iteration did not converge. (H,T) is not 00175 * in Schur form, but ALPHA(i) and BETA(i), 00176 * i=INFO+1,...,N should be correct. 00177 * = N+1,...,2*N: the shift calculation failed. (H,T) is not 00178 * in Schur form, but ALPHA(i) and BETA(i), 00179 * i=INFO-N+1,...,N should be correct. 00180 * 00181 * Further Details 00182 * =============== 00183 * 00184 * We assume that complex ABS works as long as its value is less than 00185 * overflow. 00186 * 00187 * ===================================================================== 00188 * 00189 * .. Parameters .. 00190 COMPLEX CZERO, CONE 00191 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00192 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00193 REAL ZERO, ONE 00194 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00195 REAL HALF 00196 PARAMETER ( HALF = 0.5E+0 ) 00197 * .. 00198 * .. Local Scalars .. 00199 LOGICAL ILAZR2, ILAZRO, ILQ, ILSCHR, ILZ, LQUERY 00200 INTEGER ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST, 00201 $ ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER, 00202 $ JR, MAXIT 00203 REAL ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL, 00204 $ C, SAFMIN, TEMP, TEMP2, TEMPR, ULP 00205 COMPLEX ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2, 00206 $ CTEMP3, ESHIFT, RTDISC, S, SHIFT, SIGNBC, T1, 00207 $ U12, X 00208 * .. 00209 * .. External Functions .. 00210 LOGICAL LSAME 00211 REAL CLANHS, SLAMCH 00212 EXTERNAL LSAME, CLANHS, SLAMCH 00213 * .. 00214 * .. External Subroutines .. 00215 EXTERNAL CLARTG, CLASET, CROT, CSCAL, XERBLA 00216 * .. 00217 * .. Intrinsic Functions .. 00218 INTRINSIC ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL, SQRT 00219 * .. 00220 * .. Statement Functions .. 00221 REAL ABS1 00222 * .. 00223 * .. Statement Function definitions .. 00224 ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) ) 00225 * .. 00226 * .. Executable Statements .. 00227 * 00228 * Decode JOB, COMPQ, COMPZ 00229 * 00230 IF( LSAME( JOB, 'E' ) ) THEN 00231 ILSCHR = .FALSE. 00232 ISCHUR = 1 00233 ELSE IF( LSAME( JOB, 'S' ) ) THEN 00234 ILSCHR = .TRUE. 00235 ISCHUR = 2 00236 ELSE 00237 ISCHUR = 0 00238 END IF 00239 * 00240 IF( LSAME( COMPQ, 'N' ) ) THEN 00241 ILQ = .FALSE. 00242 ICOMPQ = 1 00243 ELSE IF( LSAME( COMPQ, 'V' ) ) THEN 00244 ILQ = .TRUE. 00245 ICOMPQ = 2 00246 ELSE IF( LSAME( COMPQ, 'I' ) ) THEN 00247 ILQ = .TRUE. 00248 ICOMPQ = 3 00249 ELSE 00250 ICOMPQ = 0 00251 END IF 00252 * 00253 IF( LSAME( COMPZ, 'N' ) ) THEN 00254 ILZ = .FALSE. 00255 ICOMPZ = 1 00256 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN 00257 ILZ = .TRUE. 00258 ICOMPZ = 2 00259 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN 00260 ILZ = .TRUE. 00261 ICOMPZ = 3 00262 ELSE 00263 ICOMPZ = 0 00264 END IF 00265 * 00266 * Check Argument Values 00267 * 00268 INFO = 0 00269 WORK( 1 ) = MAX( 1, N ) 00270 LQUERY = ( LWORK.EQ.-1 ) 00271 IF( ISCHUR.EQ.0 ) THEN 00272 INFO = -1 00273 ELSE IF( ICOMPQ.EQ.0 ) THEN 00274 INFO = -2 00275 ELSE IF( ICOMPZ.EQ.0 ) THEN 00276 INFO = -3 00277 ELSE IF( N.LT.0 ) THEN 00278 INFO = -4 00279 ELSE IF( ILO.LT.1 ) THEN 00280 INFO = -5 00281 ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN 00282 INFO = -6 00283 ELSE IF( LDH.LT.N ) THEN 00284 INFO = -8 00285 ELSE IF( LDT.LT.N ) THEN 00286 INFO = -10 00287 ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN 00288 INFO = -14 00289 ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN 00290 INFO = -16 00291 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN 00292 INFO = -18 00293 END IF 00294 IF( INFO.NE.0 ) THEN 00295 CALL XERBLA( 'CHGEQZ', -INFO ) 00296 RETURN 00297 ELSE IF( LQUERY ) THEN 00298 RETURN 00299 END IF 00300 * 00301 * Quick return if possible 00302 * 00303 * WORK( 1 ) = CMPLX( 1 ) 00304 IF( N.LE.0 ) THEN 00305 WORK( 1 ) = CMPLX( 1 ) 00306 RETURN 00307 END IF 00308 * 00309 * Initialize Q and Z 00310 * 00311 IF( ICOMPQ.EQ.3 ) 00312 $ CALL CLASET( 'Full', N, N, CZERO, CONE, Q, LDQ ) 00313 IF( ICOMPZ.EQ.3 ) 00314 $ CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDZ ) 00315 * 00316 * Machine Constants 00317 * 00318 IN = IHI + 1 - ILO 00319 SAFMIN = SLAMCH( 'S' ) 00320 ULP = SLAMCH( 'E' )*SLAMCH( 'B' ) 00321 ANORM = CLANHS( 'F', IN, H( ILO, ILO ), LDH, RWORK ) 00322 BNORM = CLANHS( 'F', IN, T( ILO, ILO ), LDT, RWORK ) 00323 ATOL = MAX( SAFMIN, ULP*ANORM ) 00324 BTOL = MAX( SAFMIN, ULP*BNORM ) 00325 ASCALE = ONE / MAX( SAFMIN, ANORM ) 00326 BSCALE = ONE / MAX( SAFMIN, BNORM ) 00327 * 00328 * 00329 * Set Eigenvalues IHI+1:N 00330 * 00331 DO 10 J = IHI + 1, N 00332 ABSB = ABS( T( J, J ) ) 00333 IF( ABSB.GT.SAFMIN ) THEN 00334 SIGNBC = CONJG( T( J, J ) / ABSB ) 00335 T( J, J ) = ABSB 00336 IF( ILSCHR ) THEN 00337 CALL CSCAL( J-1, SIGNBC, T( 1, J ), 1 ) 00338 CALL CSCAL( J, SIGNBC, H( 1, J ), 1 ) 00339 ELSE 00340 H( J, J ) = H( J, J )*SIGNBC 00341 END IF 00342 IF( ILZ ) 00343 $ CALL CSCAL( N, SIGNBC, Z( 1, J ), 1 ) 00344 ELSE 00345 T( J, J ) = CZERO 00346 END IF 00347 ALPHA( J ) = H( J, J ) 00348 BETA( J ) = T( J, J ) 00349 10 CONTINUE 00350 * 00351 * If IHI < ILO, skip QZ steps 00352 * 00353 IF( IHI.LT.ILO ) 00354 $ GO TO 190 00355 * 00356 * MAIN QZ ITERATION LOOP 00357 * 00358 * Initialize dynamic indices 00359 * 00360 * Eigenvalues ILAST+1:N have been found. 00361 * Column operations modify rows IFRSTM:whatever 00362 * Row operations modify columns whatever:ILASTM 00363 * 00364 * If only eigenvalues are being computed, then 00365 * IFRSTM is the row of the last splitting row above row ILAST; 00366 * this is always at least ILO. 00367 * IITER counts iterations since the last eigenvalue was found, 00368 * to tell when to use an extraordinary shift. 00369 * MAXIT is the maximum number of QZ sweeps allowed. 00370 * 00371 ILAST = IHI 00372 IF( ILSCHR ) THEN 00373 IFRSTM = 1 00374 ILASTM = N 00375 ELSE 00376 IFRSTM = ILO 00377 ILASTM = IHI 00378 END IF 00379 IITER = 0 00380 ESHIFT = CZERO 00381 MAXIT = 30*( IHI-ILO+1 ) 00382 * 00383 DO 170 JITER = 1, MAXIT 00384 * 00385 * Check for too many iterations. 00386 * 00387 IF( JITER.GT.MAXIT ) 00388 $ GO TO 180 00389 * 00390 * Split the matrix if possible. 00391 * 00392 * Two tests: 00393 * 1: H(j,j-1)=0 or j=ILO 00394 * 2: T(j,j)=0 00395 * 00396 * Special case: j=ILAST 00397 * 00398 IF( ILAST.EQ.ILO ) THEN 00399 GO TO 60 00400 ELSE 00401 IF( ABS1( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN 00402 H( ILAST, ILAST-1 ) = CZERO 00403 GO TO 60 00404 END IF 00405 END IF 00406 * 00407 IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN 00408 T( ILAST, ILAST ) = CZERO 00409 GO TO 50 00410 END IF 00411 * 00412 * General case: j<ILAST 00413 * 00414 DO 40 J = ILAST - 1, ILO, -1 00415 * 00416 * Test 1: for H(j,j-1)=0 or j=ILO 00417 * 00418 IF( J.EQ.ILO ) THEN 00419 ILAZRO = .TRUE. 00420 ELSE 00421 IF( ABS1( H( J, J-1 ) ).LE.ATOL ) THEN 00422 H( J, J-1 ) = CZERO 00423 ILAZRO = .TRUE. 00424 ELSE 00425 ILAZRO = .FALSE. 00426 END IF 00427 END IF 00428 * 00429 * Test 2: for T(j,j)=0 00430 * 00431 IF( ABS( T( J, J ) ).LT.BTOL ) THEN 00432 T( J, J ) = CZERO 00433 * 00434 * Test 1a: Check for 2 consecutive small subdiagonals in A 00435 * 00436 ILAZR2 = .FALSE. 00437 IF( .NOT.ILAZRO ) THEN 00438 IF( ABS1( H( J, J-1 ) )*( ASCALE*ABS1( H( J+1, 00439 $ J ) ) ).LE.ABS1( H( J, J ) )*( ASCALE*ATOL ) ) 00440 $ ILAZR2 = .TRUE. 00441 END IF 00442 * 00443 * If both tests pass (1 & 2), i.e., the leading diagonal 00444 * element of B in the block is zero, split a 1x1 block off 00445 * at the top. (I.e., at the J-th row/column) The leading 00446 * diagonal element of the remainder can also be zero, so 00447 * this may have to be done repeatedly. 00448 * 00449 IF( ILAZRO .OR. ILAZR2 ) THEN 00450 DO 20 JCH = J, ILAST - 1 00451 CTEMP = H( JCH, JCH ) 00452 CALL CLARTG( CTEMP, H( JCH+1, JCH ), C, S, 00453 $ H( JCH, JCH ) ) 00454 H( JCH+1, JCH ) = CZERO 00455 CALL CROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH, 00456 $ H( JCH+1, JCH+1 ), LDH, C, S ) 00457 CALL CROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT, 00458 $ T( JCH+1, JCH+1 ), LDT, C, S ) 00459 IF( ILQ ) 00460 $ CALL CROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1, 00461 $ C, CONJG( S ) ) 00462 IF( ILAZR2 ) 00463 $ H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C 00464 ILAZR2 = .FALSE. 00465 IF( ABS1( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN 00466 IF( JCH+1.GE.ILAST ) THEN 00467 GO TO 60 00468 ELSE 00469 IFIRST = JCH + 1 00470 GO TO 70 00471 END IF 00472 END IF 00473 T( JCH+1, JCH+1 ) = CZERO 00474 20 CONTINUE 00475 GO TO 50 00476 ELSE 00477 * 00478 * Only test 2 passed -- chase the zero to T(ILAST,ILAST) 00479 * Then process as in the case T(ILAST,ILAST)=0 00480 * 00481 DO 30 JCH = J, ILAST - 1 00482 CTEMP = T( JCH, JCH+1 ) 00483 CALL CLARTG( CTEMP, T( JCH+1, JCH+1 ), C, S, 00484 $ T( JCH, JCH+1 ) ) 00485 T( JCH+1, JCH+1 ) = CZERO 00486 IF( JCH.LT.ILASTM-1 ) 00487 $ CALL CROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT, 00488 $ T( JCH+1, JCH+2 ), LDT, C, S ) 00489 CALL CROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH, 00490 $ H( JCH+1, JCH-1 ), LDH, C, S ) 00491 IF( ILQ ) 00492 $ CALL CROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1, 00493 $ C, CONJG( S ) ) 00494 CTEMP = H( JCH+1, JCH ) 00495 CALL CLARTG( CTEMP, H( JCH+1, JCH-1 ), C, S, 00496 $ H( JCH+1, JCH ) ) 00497 H( JCH+1, JCH-1 ) = CZERO 00498 CALL CROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1, 00499 $ H( IFRSTM, JCH-1 ), 1, C, S ) 00500 CALL CROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1, 00501 $ T( IFRSTM, JCH-1 ), 1, C, S ) 00502 IF( ILZ ) 00503 $ CALL CROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1, 00504 $ C, S ) 00505 30 CONTINUE 00506 GO TO 50 00507 END IF 00508 ELSE IF( ILAZRO ) THEN 00509 * 00510 * Only test 1 passed -- work on J:ILAST 00511 * 00512 IFIRST = J 00513 GO TO 70 00514 END IF 00515 * 00516 * Neither test passed -- try next J 00517 * 00518 40 CONTINUE 00519 * 00520 * (Drop-through is "impossible") 00521 * 00522 INFO = 2*N + 1 00523 GO TO 210 00524 * 00525 * T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a 00526 * 1x1 block. 00527 * 00528 50 CONTINUE 00529 CTEMP = H( ILAST, ILAST ) 00530 CALL CLARTG( CTEMP, H( ILAST, ILAST-1 ), C, S, 00531 $ H( ILAST, ILAST ) ) 00532 H( ILAST, ILAST-1 ) = CZERO 00533 CALL CROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1, 00534 $ H( IFRSTM, ILAST-1 ), 1, C, S ) 00535 CALL CROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1, 00536 $ T( IFRSTM, ILAST-1 ), 1, C, S ) 00537 IF( ILZ ) 00538 $ CALL CROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S ) 00539 * 00540 * H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHA and BETA 00541 * 00542 60 CONTINUE 00543 ABSB = ABS( T( ILAST, ILAST ) ) 00544 IF( ABSB.GT.SAFMIN ) THEN 00545 SIGNBC = CONJG( T( ILAST, ILAST ) / ABSB ) 00546 T( ILAST, ILAST ) = ABSB 00547 IF( ILSCHR ) THEN 00548 CALL CSCAL( ILAST-IFRSTM, SIGNBC, T( IFRSTM, ILAST ), 1 ) 00549 CALL CSCAL( ILAST+1-IFRSTM, SIGNBC, H( IFRSTM, ILAST ), 00550 $ 1 ) 00551 ELSE 00552 H( ILAST, ILAST ) = H( ILAST, ILAST )*SIGNBC 00553 END IF 00554 IF( ILZ ) 00555 $ CALL CSCAL( N, SIGNBC, Z( 1, ILAST ), 1 ) 00556 ELSE 00557 T( ILAST, ILAST ) = CZERO 00558 END IF 00559 ALPHA( ILAST ) = H( ILAST, ILAST ) 00560 BETA( ILAST ) = T( ILAST, ILAST ) 00561 * 00562 * Go to next block -- exit if finished. 00563 * 00564 ILAST = ILAST - 1 00565 IF( ILAST.LT.ILO ) 00566 $ GO TO 190 00567 * 00568 * Reset counters 00569 * 00570 IITER = 0 00571 ESHIFT = CZERO 00572 IF( .NOT.ILSCHR ) THEN 00573 ILASTM = ILAST 00574 IF( IFRSTM.GT.ILAST ) 00575 $ IFRSTM = ILO 00576 END IF 00577 GO TO 160 00578 * 00579 * QZ step 00580 * 00581 * This iteration only involves rows/columns IFIRST:ILAST. We 00582 * assume IFIRST < ILAST, and that the diagonal of B is non-zero. 00583 * 00584 70 CONTINUE 00585 IITER = IITER + 1 00586 IF( .NOT.ILSCHR ) THEN 00587 IFRSTM = IFIRST 00588 END IF 00589 * 00590 * Compute the Shift. 00591 * 00592 * At this point, IFIRST < ILAST, and the diagonal elements of 00593 * T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in 00594 * magnitude) 00595 * 00596 IF( ( IITER / 10 )*10.NE.IITER ) THEN 00597 * 00598 * The Wilkinson shift (AEP p.512), i.e., the eigenvalue of 00599 * the bottom-right 2x2 block of A inv(B) which is nearest to 00600 * the bottom-right element. 00601 * 00602 * We factor B as U*D, where U has unit diagonals, and 00603 * compute (A*inv(D))*inv(U). 00604 * 00605 U12 = ( BSCALE*T( ILAST-1, ILAST ) ) / 00606 $ ( BSCALE*T( ILAST, ILAST ) ) 00607 AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) / 00608 $ ( BSCALE*T( ILAST-1, ILAST-1 ) ) 00609 AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) / 00610 $ ( BSCALE*T( ILAST-1, ILAST-1 ) ) 00611 AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) / 00612 $ ( BSCALE*T( ILAST, ILAST ) ) 00613 AD22 = ( ASCALE*H( ILAST, ILAST ) ) / 00614 $ ( BSCALE*T( ILAST, ILAST ) ) 00615 ABI22 = AD22 - U12*AD21 00616 * 00617 T1 = HALF*( AD11+ABI22 ) 00618 RTDISC = SQRT( T1**2+AD12*AD21-AD11*AD22 ) 00619 TEMP = REAL( T1-ABI22 )*REAL( RTDISC ) + 00620 $ AIMAG( T1-ABI22 )*AIMAG( RTDISC ) 00621 IF( TEMP.LE.ZERO ) THEN 00622 SHIFT = T1 + RTDISC 00623 ELSE 00624 SHIFT = T1 - RTDISC 00625 END IF 00626 ELSE 00627 * 00628 * Exceptional shift. Chosen for no particularly good reason. 00629 * 00630 ESHIFT = ESHIFT + CONJG( ( ASCALE*H( ILAST-1, ILAST ) ) / 00631 $ ( BSCALE*T( ILAST-1, ILAST-1 ) ) ) 00632 SHIFT = ESHIFT 00633 END IF 00634 * 00635 * Now check for two consecutive small subdiagonals. 00636 * 00637 DO 80 J = ILAST - 1, IFIRST + 1, -1 00638 ISTART = J 00639 CTEMP = ASCALE*H( J, J ) - SHIFT*( BSCALE*T( J, J ) ) 00640 TEMP = ABS1( CTEMP ) 00641 TEMP2 = ASCALE*ABS1( H( J+1, J ) ) 00642 TEMPR = MAX( TEMP, TEMP2 ) 00643 IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN 00644 TEMP = TEMP / TEMPR 00645 TEMP2 = TEMP2 / TEMPR 00646 END IF 00647 IF( ABS1( H( J, J-1 ) )*TEMP2.LE.TEMP*ATOL ) 00648 $ GO TO 90 00649 80 CONTINUE 00650 * 00651 ISTART = IFIRST 00652 CTEMP = ASCALE*H( IFIRST, IFIRST ) - 00653 $ SHIFT*( BSCALE*T( IFIRST, IFIRST ) ) 00654 90 CONTINUE 00655 * 00656 * Do an implicit-shift QZ sweep. 00657 * 00658 * Initial Q 00659 * 00660 CTEMP2 = ASCALE*H( ISTART+1, ISTART ) 00661 CALL CLARTG( CTEMP, CTEMP2, C, S, CTEMP3 ) 00662 * 00663 * Sweep 00664 * 00665 DO 150 J = ISTART, ILAST - 1 00666 IF( J.GT.ISTART ) THEN 00667 CTEMP = H( J, J-1 ) 00668 CALL CLARTG( CTEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) ) 00669 H( J+1, J-1 ) = CZERO 00670 END IF 00671 * 00672 DO 100 JC = J, ILASTM 00673 CTEMP = C*H( J, JC ) + S*H( J+1, JC ) 00674 H( J+1, JC ) = -CONJG( S )*H( J, JC ) + C*H( J+1, JC ) 00675 H( J, JC ) = CTEMP 00676 CTEMP2 = C*T( J, JC ) + S*T( J+1, JC ) 00677 T( J+1, JC ) = -CONJG( S )*T( J, JC ) + C*T( J+1, JC ) 00678 T( J, JC ) = CTEMP2 00679 100 CONTINUE 00680 IF( ILQ ) THEN 00681 DO 110 JR = 1, N 00682 CTEMP = C*Q( JR, J ) + CONJG( S )*Q( JR, J+1 ) 00683 Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 ) 00684 Q( JR, J ) = CTEMP 00685 110 CONTINUE 00686 END IF 00687 * 00688 CTEMP = T( J+1, J+1 ) 00689 CALL CLARTG( CTEMP, T( J+1, J ), C, S, T( J+1, J+1 ) ) 00690 T( J+1, J ) = CZERO 00691 * 00692 DO 120 JR = IFRSTM, MIN( J+2, ILAST ) 00693 CTEMP = C*H( JR, J+1 ) + S*H( JR, J ) 00694 H( JR, J ) = -CONJG( S )*H( JR, J+1 ) + C*H( JR, J ) 00695 H( JR, J+1 ) = CTEMP 00696 120 CONTINUE 00697 DO 130 JR = IFRSTM, J 00698 CTEMP = C*T( JR, J+1 ) + S*T( JR, J ) 00699 T( JR, J ) = -CONJG( S )*T( JR, J+1 ) + C*T( JR, J ) 00700 T( JR, J+1 ) = CTEMP 00701 130 CONTINUE 00702 IF( ILZ ) THEN 00703 DO 140 JR = 1, N 00704 CTEMP = C*Z( JR, J+1 ) + S*Z( JR, J ) 00705 Z( JR, J ) = -CONJG( S )*Z( JR, J+1 ) + C*Z( JR, J ) 00706 Z( JR, J+1 ) = CTEMP 00707 140 CONTINUE 00708 END IF 00709 150 CONTINUE 00710 * 00711 160 CONTINUE 00712 * 00713 170 CONTINUE 00714 * 00715 * Drop-through = non-convergence 00716 * 00717 180 CONTINUE 00718 INFO = ILAST 00719 GO TO 210 00720 * 00721 * Successful completion of all QZ steps 00722 * 00723 190 CONTINUE 00724 * 00725 * Set Eigenvalues 1:ILO-1 00726 * 00727 DO 200 J = 1, ILO - 1 00728 ABSB = ABS( T( J, J ) ) 00729 IF( ABSB.GT.SAFMIN ) THEN 00730 SIGNBC = CONJG( T( J, J ) / ABSB ) 00731 T( J, J ) = ABSB 00732 IF( ILSCHR ) THEN 00733 CALL CSCAL( J-1, SIGNBC, T( 1, J ), 1 ) 00734 CALL CSCAL( J, SIGNBC, H( 1, J ), 1 ) 00735 ELSE 00736 H( J, J ) = H( J, J )*SIGNBC 00737 END IF 00738 IF( ILZ ) 00739 $ CALL CSCAL( N, SIGNBC, Z( 1, J ), 1 ) 00740 ELSE 00741 T( J, J ) = CZERO 00742 END IF 00743 ALPHA( J ) = H( J, J ) 00744 BETA( J ) = T( J, J ) 00745 200 CONTINUE 00746 * 00747 * Normal Termination 00748 * 00749 INFO = 0 00750 * 00751 * Exit (other than argument error) -- return optimal workspace size 00752 * 00753 210 CONTINUE 00754 WORK( 1 ) = CMPLX( N ) 00755 RETURN 00756 * 00757 * End of CHGEQZ 00758 * 00759 END
1.7.3 
