|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE ZHGEQZ( 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.2) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * November 2006 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 DOUBLE PRECISION RWORK( * ) 00016 COMPLEX*16 ALPHA( * ), BETA( * ), H( LDH, * ), 00017 $ Q( LDQ, * ), T( LDT, * ), WORK( * ), 00018 $ Z( LDZ, * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * ZHGEQZ 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 ZGGHRD. 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 ZGGHRD 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*16 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*16 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*16 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*16 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*16 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*16 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*16 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) DOUBLE PRECISION 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*16 CZERO, CONE 00191 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00192 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00193 DOUBLE PRECISION ZERO, ONE 00194 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00195 DOUBLE PRECISION HALF 00196 PARAMETER ( HALF = 0.5D+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 DOUBLE PRECISION ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL, 00204 $ C, SAFMIN, TEMP, TEMP2, TEMPR, ULP 00205 COMPLEX*16 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 DOUBLE PRECISION DLAMCH, ZLANHS 00212 EXTERNAL LSAME, DLAMCH, ZLANHS 00213 * .. 00214 * .. External Subroutines .. 00215 EXTERNAL XERBLA, ZLARTG, ZLASET, ZROT, ZSCAL 00216 * .. 00217 * .. Intrinsic Functions .. 00218 INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN, 00219 $ SQRT 00220 * .. 00221 * .. Statement Functions .. 00222 DOUBLE PRECISION ABS1 00223 * .. 00224 * .. Statement Function definitions .. 00225 ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) ) 00226 * .. 00227 * .. Executable Statements .. 00228 * 00229 * Decode JOB, COMPQ, COMPZ 00230 * 00231 IF( LSAME( JOB, 'E' ) ) THEN 00232 ILSCHR = .FALSE. 00233 ISCHUR = 1 00234 ELSE IF( LSAME( JOB, 'S' ) ) THEN 00235 ILSCHR = .TRUE. 00236 ISCHUR = 2 00237 ELSE 00238 ISCHUR = 0 00239 END IF 00240 * 00241 IF( LSAME( COMPQ, 'N' ) ) THEN 00242 ILQ = .FALSE. 00243 ICOMPQ = 1 00244 ELSE IF( LSAME( COMPQ, 'V' ) ) THEN 00245 ILQ = .TRUE. 00246 ICOMPQ = 2 00247 ELSE IF( LSAME( COMPQ, 'I' ) ) THEN 00248 ILQ = .TRUE. 00249 ICOMPQ = 3 00250 ELSE 00251 ICOMPQ = 0 00252 END IF 00253 * 00254 IF( LSAME( COMPZ, 'N' ) ) THEN 00255 ILZ = .FALSE. 00256 ICOMPZ = 1 00257 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN 00258 ILZ = .TRUE. 00259 ICOMPZ = 2 00260 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN 00261 ILZ = .TRUE. 00262 ICOMPZ = 3 00263 ELSE 00264 ICOMPZ = 0 00265 END IF 00266 * 00267 * Check Argument Values 00268 * 00269 INFO = 0 00270 WORK( 1 ) = MAX( 1, N ) 00271 LQUERY = ( LWORK.EQ.-1 ) 00272 IF( ISCHUR.EQ.0 ) THEN 00273 INFO = -1 00274 ELSE IF( ICOMPQ.EQ.0 ) THEN 00275 INFO = -2 00276 ELSE IF( ICOMPZ.EQ.0 ) THEN 00277 INFO = -3 00278 ELSE IF( N.LT.0 ) THEN 00279 INFO = -4 00280 ELSE IF( ILO.LT.1 ) THEN 00281 INFO = -5 00282 ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN 00283 INFO = -6 00284 ELSE IF( LDH.LT.N ) THEN 00285 INFO = -8 00286 ELSE IF( LDT.LT.N ) THEN 00287 INFO = -10 00288 ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN 00289 INFO = -14 00290 ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN 00291 INFO = -16 00292 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN 00293 INFO = -18 00294 END IF 00295 IF( INFO.NE.0 ) THEN 00296 CALL XERBLA( 'ZHGEQZ', -INFO ) 00297 RETURN 00298 ELSE IF( LQUERY ) THEN 00299 RETURN 00300 END IF 00301 * 00302 * Quick return if possible 00303 * 00304 * WORK( 1 ) = CMPLX( 1 ) 00305 IF( N.LE.0 ) THEN 00306 WORK( 1 ) = DCMPLX( 1 ) 00307 RETURN 00308 END IF 00309 * 00310 * Initialize Q and Z 00311 * 00312 IF( ICOMPQ.EQ.3 ) 00313 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ ) 00314 IF( ICOMPZ.EQ.3 ) 00315 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ ) 00316 * 00317 * Machine Constants 00318 * 00319 IN = IHI + 1 - ILO 00320 SAFMIN = DLAMCH( 'S' ) 00321 ULP = DLAMCH( 'E' )*DLAMCH( 'B' ) 00322 ANORM = ZLANHS( 'F', IN, H( ILO, ILO ), LDH, RWORK ) 00323 BNORM = ZLANHS( 'F', IN, T( ILO, ILO ), LDT, RWORK ) 00324 ATOL = MAX( SAFMIN, ULP*ANORM ) 00325 BTOL = MAX( SAFMIN, ULP*BNORM ) 00326 ASCALE = ONE / MAX( SAFMIN, ANORM ) 00327 BSCALE = ONE / MAX( SAFMIN, BNORM ) 00328 * 00329 * 00330 * Set Eigenvalues IHI+1:N 00331 * 00332 DO 10 J = IHI + 1, N 00333 ABSB = ABS( T( J, J ) ) 00334 IF( ABSB.GT.SAFMIN ) THEN 00335 SIGNBC = DCONJG( T( J, J ) / ABSB ) 00336 T( J, J ) = ABSB 00337 IF( ILSCHR ) THEN 00338 CALL ZSCAL( J-1, SIGNBC, T( 1, J ), 1 ) 00339 CALL ZSCAL( J, SIGNBC, H( 1, J ), 1 ) 00340 ELSE 00341 H( J, J ) = H( J, J )*SIGNBC 00342 END IF 00343 IF( ILZ ) 00344 $ CALL ZSCAL( N, SIGNBC, Z( 1, J ), 1 ) 00345 ELSE 00346 T( J, J ) = CZERO 00347 END IF 00348 ALPHA( J ) = H( J, J ) 00349 BETA( J ) = T( J, J ) 00350 10 CONTINUE 00351 * 00352 * If IHI < ILO, skip QZ steps 00353 * 00354 IF( IHI.LT.ILO ) 00355 $ GO TO 190 00356 * 00357 * MAIN QZ ITERATION LOOP 00358 * 00359 * Initialize dynamic indices 00360 * 00361 * Eigenvalues ILAST+1:N have been found. 00362 * Column operations modify rows IFRSTM:whatever 00363 * Row operations modify columns whatever:ILASTM 00364 * 00365 * If only eigenvalues are being computed, then 00366 * IFRSTM is the row of the last splitting row above row ILAST; 00367 * this is always at least ILO. 00368 * IITER counts iterations since the last eigenvalue was found, 00369 * to tell when to use an extraordinary shift. 00370 * MAXIT is the maximum number of QZ sweeps allowed. 00371 * 00372 ILAST = IHI 00373 IF( ILSCHR ) THEN 00374 IFRSTM = 1 00375 ILASTM = N 00376 ELSE 00377 IFRSTM = ILO 00378 ILASTM = IHI 00379 END IF 00380 IITER = 0 00381 ESHIFT = CZERO 00382 MAXIT = 30*( IHI-ILO+1 ) 00383 * 00384 DO 170 JITER = 1, MAXIT 00385 * 00386 * Check for too many iterations. 00387 * 00388 IF( JITER.GT.MAXIT ) 00389 $ GO TO 180 00390 * 00391 * Split the matrix if possible. 00392 * 00393 * Two tests: 00394 * 1: H(j,j-1)=0 or j=ILO 00395 * 2: T(j,j)=0 00396 * 00397 * Special case: j=ILAST 00398 * 00399 IF( ILAST.EQ.ILO ) THEN 00400 GO TO 60 00401 ELSE 00402 IF( ABS1( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN 00403 H( ILAST, ILAST-1 ) = CZERO 00404 GO TO 60 00405 END IF 00406 END IF 00407 * 00408 IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN 00409 T( ILAST, ILAST ) = CZERO 00410 GO TO 50 00411 END IF 00412 * 00413 * General case: j<ILAST 00414 * 00415 DO 40 J = ILAST - 1, ILO, -1 00416 * 00417 * Test 1: for H(j,j-1)=0 or j=ILO 00418 * 00419 IF( J.EQ.ILO ) THEN 00420 ILAZRO = .TRUE. 00421 ELSE 00422 IF( ABS1( H( J, J-1 ) ).LE.ATOL ) THEN 00423 H( J, J-1 ) = CZERO 00424 ILAZRO = .TRUE. 00425 ELSE 00426 ILAZRO = .FALSE. 00427 END IF 00428 END IF 00429 * 00430 * Test 2: for T(j,j)=0 00431 * 00432 IF( ABS( T( J, J ) ).LT.BTOL ) THEN 00433 T( J, J ) = CZERO 00434 * 00435 * Test 1a: Check for 2 consecutive small subdiagonals in A 00436 * 00437 ILAZR2 = .FALSE. 00438 IF( .NOT.ILAZRO ) THEN 00439 IF( ABS1( H( J, J-1 ) )*( ASCALE*ABS1( H( J+1, 00440 $ J ) ) ).LE.ABS1( H( J, J ) )*( ASCALE*ATOL ) ) 00441 $ ILAZR2 = .TRUE. 00442 END IF 00443 * 00444 * If both tests pass (1 & 2), i.e., the leading diagonal 00445 * element of B in the block is zero, split a 1x1 block off 00446 * at the top. (I.e., at the J-th row/column) The leading 00447 * diagonal element of the remainder can also be zero, so 00448 * this may have to be done repeatedly. 00449 * 00450 IF( ILAZRO .OR. ILAZR2 ) THEN 00451 DO 20 JCH = J, ILAST - 1 00452 CTEMP = H( JCH, JCH ) 00453 CALL ZLARTG( CTEMP, H( JCH+1, JCH ), C, S, 00454 $ H( JCH, JCH ) ) 00455 H( JCH+1, JCH ) = CZERO 00456 CALL ZROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH, 00457 $ H( JCH+1, JCH+1 ), LDH, C, S ) 00458 CALL ZROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT, 00459 $ T( JCH+1, JCH+1 ), LDT, C, S ) 00460 IF( ILQ ) 00461 $ CALL ZROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1, 00462 $ C, DCONJG( S ) ) 00463 IF( ILAZR2 ) 00464 $ H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C 00465 ILAZR2 = .FALSE. 00466 IF( ABS1( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN 00467 IF( JCH+1.GE.ILAST ) THEN 00468 GO TO 60 00469 ELSE 00470 IFIRST = JCH + 1 00471 GO TO 70 00472 END IF 00473 END IF 00474 T( JCH+1, JCH+1 ) = CZERO 00475 20 CONTINUE 00476 GO TO 50 00477 ELSE 00478 * 00479 * Only test 2 passed -- chase the zero to T(ILAST,ILAST) 00480 * Then process as in the case T(ILAST,ILAST)=0 00481 * 00482 DO 30 JCH = J, ILAST - 1 00483 CTEMP = T( JCH, JCH+1 ) 00484 CALL ZLARTG( CTEMP, T( JCH+1, JCH+1 ), C, S, 00485 $ T( JCH, JCH+1 ) ) 00486 T( JCH+1, JCH+1 ) = CZERO 00487 IF( JCH.LT.ILASTM-1 ) 00488 $ CALL ZROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT, 00489 $ T( JCH+1, JCH+2 ), LDT, C, S ) 00490 CALL ZROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH, 00491 $ H( JCH+1, JCH-1 ), LDH, C, S ) 00492 IF( ILQ ) 00493 $ CALL ZROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1, 00494 $ C, DCONJG( S ) ) 00495 CTEMP = H( JCH+1, JCH ) 00496 CALL ZLARTG( CTEMP, H( JCH+1, JCH-1 ), C, S, 00497 $ H( JCH+1, JCH ) ) 00498 H( JCH+1, JCH-1 ) = CZERO 00499 CALL ZROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1, 00500 $ H( IFRSTM, JCH-1 ), 1, C, S ) 00501 CALL ZROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1, 00502 $ T( IFRSTM, JCH-1 ), 1, C, S ) 00503 IF( ILZ ) 00504 $ CALL ZROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1, 00505 $ C, S ) 00506 30 CONTINUE 00507 GO TO 50 00508 END IF 00509 ELSE IF( ILAZRO ) THEN 00510 * 00511 * Only test 1 passed -- work on J:ILAST 00512 * 00513 IFIRST = J 00514 GO TO 70 00515 END IF 00516 * 00517 * Neither test passed -- try next J 00518 * 00519 40 CONTINUE 00520 * 00521 * (Drop-through is "impossible") 00522 * 00523 INFO = 2*N + 1 00524 GO TO 210 00525 * 00526 * T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a 00527 * 1x1 block. 00528 * 00529 50 CONTINUE 00530 CTEMP = H( ILAST, ILAST ) 00531 CALL ZLARTG( CTEMP, H( ILAST, ILAST-1 ), C, S, 00532 $ H( ILAST, ILAST ) ) 00533 H( ILAST, ILAST-1 ) = CZERO 00534 CALL ZROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1, 00535 $ H( IFRSTM, ILAST-1 ), 1, C, S ) 00536 CALL ZROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1, 00537 $ T( IFRSTM, ILAST-1 ), 1, C, S ) 00538 IF( ILZ ) 00539 $ CALL ZROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S ) 00540 * 00541 * H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHA and BETA 00542 * 00543 60 CONTINUE 00544 ABSB = ABS( T( ILAST, ILAST ) ) 00545 IF( ABSB.GT.SAFMIN ) THEN 00546 SIGNBC = DCONJG( T( ILAST, ILAST ) / ABSB ) 00547 T( ILAST, ILAST ) = ABSB 00548 IF( ILSCHR ) THEN 00549 CALL ZSCAL( ILAST-IFRSTM, SIGNBC, T( IFRSTM, ILAST ), 1 ) 00550 CALL ZSCAL( ILAST+1-IFRSTM, SIGNBC, H( IFRSTM, ILAST ), 00551 $ 1 ) 00552 ELSE 00553 H( ILAST, ILAST ) = H( ILAST, ILAST )*SIGNBC 00554 END IF 00555 IF( ILZ ) 00556 $ CALL ZSCAL( N, SIGNBC, Z( 1, ILAST ), 1 ) 00557 ELSE 00558 T( ILAST, ILAST ) = CZERO 00559 END IF 00560 ALPHA( ILAST ) = H( ILAST, ILAST ) 00561 BETA( ILAST ) = T( ILAST, ILAST ) 00562 * 00563 * Go to next block -- exit if finished. 00564 * 00565 ILAST = ILAST - 1 00566 IF( ILAST.LT.ILO ) 00567 $ GO TO 190 00568 * 00569 * Reset counters 00570 * 00571 IITER = 0 00572 ESHIFT = CZERO 00573 IF( .NOT.ILSCHR ) THEN 00574 ILASTM = ILAST 00575 IF( IFRSTM.GT.ILAST ) 00576 $ IFRSTM = ILO 00577 END IF 00578 GO TO 160 00579 * 00580 * QZ step 00581 * 00582 * This iteration only involves rows/columns IFIRST:ILAST. We 00583 * assume IFIRST < ILAST, and that the diagonal of B is non-zero. 00584 * 00585 70 CONTINUE 00586 IITER = IITER + 1 00587 IF( .NOT.ILSCHR ) THEN 00588 IFRSTM = IFIRST 00589 END IF 00590 * 00591 * Compute the Shift. 00592 * 00593 * At this point, IFIRST < ILAST, and the diagonal elements of 00594 * T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in 00595 * magnitude) 00596 * 00597 IF( ( IITER / 10 )*10.NE.IITER ) THEN 00598 * 00599 * The Wilkinson shift (AEP p.512), i.e., the eigenvalue of 00600 * the bottom-right 2x2 block of A inv(B) which is nearest to 00601 * the bottom-right element. 00602 * 00603 * We factor B as U*D, where U has unit diagonals, and 00604 * compute (A*inv(D))*inv(U). 00605 * 00606 U12 = ( BSCALE*T( ILAST-1, ILAST ) ) / 00607 $ ( BSCALE*T( ILAST, ILAST ) ) 00608 AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) / 00609 $ ( BSCALE*T( ILAST-1, ILAST-1 ) ) 00610 AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) / 00611 $ ( BSCALE*T( ILAST-1, ILAST-1 ) ) 00612 AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) / 00613 $ ( BSCALE*T( ILAST, ILAST ) ) 00614 AD22 = ( ASCALE*H( ILAST, ILAST ) ) / 00615 $ ( BSCALE*T( ILAST, ILAST ) ) 00616 ABI22 = AD22 - U12*AD21 00617 * 00618 T1 = HALF*( AD11+ABI22 ) 00619 RTDISC = SQRT( T1**2+AD12*AD21-AD11*AD22 ) 00620 TEMP = DBLE( T1-ABI22 )*DBLE( RTDISC ) + 00621 $ DIMAG( T1-ABI22 )*DIMAG( RTDISC ) 00622 IF( TEMP.LE.ZERO ) THEN 00623 SHIFT = T1 + RTDISC 00624 ELSE 00625 SHIFT = T1 - RTDISC 00626 END IF 00627 ELSE 00628 * 00629 * Exceptional shift. Chosen for no particularly good reason. 00630 * 00631 ESHIFT = ESHIFT + DCONJG( ( ASCALE*H( ILAST-1, ILAST ) ) / 00632 $ ( BSCALE*T( ILAST-1, ILAST-1 ) ) ) 00633 SHIFT = ESHIFT 00634 END IF 00635 * 00636 * Now check for two consecutive small subdiagonals. 00637 * 00638 DO 80 J = ILAST - 1, IFIRST + 1, -1 00639 ISTART = J 00640 CTEMP = ASCALE*H( J, J ) - SHIFT*( BSCALE*T( J, J ) ) 00641 TEMP = ABS1( CTEMP ) 00642 TEMP2 = ASCALE*ABS1( H( J+1, J ) ) 00643 TEMPR = MAX( TEMP, TEMP2 ) 00644 IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN 00645 TEMP = TEMP / TEMPR 00646 TEMP2 = TEMP2 / TEMPR 00647 END IF 00648 IF( ABS1( H( J, J-1 ) )*TEMP2.LE.TEMP*ATOL ) 00649 $ GO TO 90 00650 80 CONTINUE 00651 * 00652 ISTART = IFIRST 00653 CTEMP = ASCALE*H( IFIRST, IFIRST ) - 00654 $ SHIFT*( BSCALE*T( IFIRST, IFIRST ) ) 00655 90 CONTINUE 00656 * 00657 * Do an implicit-shift QZ sweep. 00658 * 00659 * Initial Q 00660 * 00661 CTEMP2 = ASCALE*H( ISTART+1, ISTART ) 00662 CALL ZLARTG( CTEMP, CTEMP2, C, S, CTEMP3 ) 00663 * 00664 * Sweep 00665 * 00666 DO 150 J = ISTART, ILAST - 1 00667 IF( J.GT.ISTART ) THEN 00668 CTEMP = H( J, J-1 ) 00669 CALL ZLARTG( CTEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) ) 00670 H( J+1, J-1 ) = CZERO 00671 END IF 00672 * 00673 DO 100 JC = J, ILASTM 00674 CTEMP = C*H( J, JC ) + S*H( J+1, JC ) 00675 H( J+1, JC ) = -DCONJG( S )*H( J, JC ) + C*H( J+1, JC ) 00676 H( J, JC ) = CTEMP 00677 CTEMP2 = C*T( J, JC ) + S*T( J+1, JC ) 00678 T( J+1, JC ) = -DCONJG( S )*T( J, JC ) + C*T( J+1, JC ) 00679 T( J, JC ) = CTEMP2 00680 100 CONTINUE 00681 IF( ILQ ) THEN 00682 DO 110 JR = 1, N 00683 CTEMP = C*Q( JR, J ) + DCONJG( S )*Q( JR, J+1 ) 00684 Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 ) 00685 Q( JR, J ) = CTEMP 00686 110 CONTINUE 00687 END IF 00688 * 00689 CTEMP = T( J+1, J+1 ) 00690 CALL ZLARTG( CTEMP, T( J+1, J ), C, S, T( J+1, J+1 ) ) 00691 T( J+1, J ) = CZERO 00692 * 00693 DO 120 JR = IFRSTM, MIN( J+2, ILAST ) 00694 CTEMP = C*H( JR, J+1 ) + S*H( JR, J ) 00695 H( JR, J ) = -DCONJG( S )*H( JR, J+1 ) + C*H( JR, J ) 00696 H( JR, J+1 ) = CTEMP 00697 120 CONTINUE 00698 DO 130 JR = IFRSTM, J 00699 CTEMP = C*T( JR, J+1 ) + S*T( JR, J ) 00700 T( JR, J ) = -DCONJG( S )*T( JR, J+1 ) + C*T( JR, J ) 00701 T( JR, J+1 ) = CTEMP 00702 130 CONTINUE 00703 IF( ILZ ) THEN 00704 DO 140 JR = 1, N 00705 CTEMP = C*Z( JR, J+1 ) + S*Z( JR, J ) 00706 Z( JR, J ) = -DCONJG( S )*Z( JR, J+1 ) + C*Z( JR, J ) 00707 Z( JR, J+1 ) = CTEMP 00708 140 CONTINUE 00709 END IF 00710 150 CONTINUE 00711 * 00712 160 CONTINUE 00713 * 00714 170 CONTINUE 00715 * 00716 * Drop-through = non-convergence 00717 * 00718 180 CONTINUE 00719 INFO = ILAST 00720 GO TO 210 00721 * 00722 * Successful completion of all QZ steps 00723 * 00724 190 CONTINUE 00725 * 00726 * Set Eigenvalues 1:ILO-1 00727 * 00728 DO 200 J = 1, ILO - 1 00729 ABSB = ABS( T( J, J ) ) 00730 IF( ABSB.GT.SAFMIN ) THEN 00731 SIGNBC = DCONJG( T( J, J ) / ABSB ) 00732 T( J, J ) = ABSB 00733 IF( ILSCHR ) THEN 00734 CALL ZSCAL( J-1, SIGNBC, T( 1, J ), 1 ) 00735 CALL ZSCAL( J, SIGNBC, H( 1, J ), 1 ) 00736 ELSE 00737 H( J, J ) = H( J, J )*SIGNBC 00738 END IF 00739 IF( ILZ ) 00740 $ CALL ZSCAL( N, SIGNBC, Z( 1, J ), 1 ) 00741 ELSE 00742 T( J, J ) = CZERO 00743 END IF 00744 ALPHA( J ) = H( J, J ) 00745 BETA( J ) = T( J, J ) 00746 200 CONTINUE 00747 * 00748 * Normal Termination 00749 * 00750 INFO = 0 00751 * 00752 * Exit (other than argument error) -- return optimal workspace size 00753 * 00754 210 CONTINUE 00755 WORK( 1 ) = DCMPLX( N ) 00756 RETURN 00757 * 00758 * End of ZHGEQZ 00759 * 00760 END
1.7.3 
