|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, 00002 $ ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, 00003 $ LWORK, INFO ) 00004 * 00005 * -- LAPACK routine (version 3.2.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 2009 -- 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 ALPHAI( * ), ALPHAR( * ), BETA( * ), 00016 $ H( LDH, * ), Q( LDQ, * ), T( LDT, * ), 00017 $ WORK( * ), Z( LDZ, * ) 00018 * .. 00019 * 00020 * Purpose 00021 * ======= 00022 * 00023 * DHGEQZ computes the eigenvalues of a real matrix pair (H,T), 00024 * where H is an upper Hessenberg matrix and T is upper triangular, 00025 * using the double-shift QZ method. 00026 * Matrix pairs of this type are produced by the reduction to 00027 * generalized upper Hessenberg form of a real matrix pair (A,B): 00028 * 00029 * A = Q1*H*Z1**T, B = Q1*T*Z1**T, 00030 * 00031 * as computed by DGGHRD. 00032 * 00033 * If JOB='S', then the Hessenberg-triangular pair (H,T) is 00034 * also reduced to generalized Schur form, 00035 * 00036 * H = Q*S*Z**T, T = Q*P*Z**T, 00037 * 00038 * where Q and Z are orthogonal matrices, P is an upper triangular 00039 * matrix, and S is a quasi-triangular matrix with 1-by-1 and 2-by-2 00040 * diagonal blocks. 00041 * 00042 * The 1-by-1 blocks correspond to real eigenvalues of the matrix pair 00043 * (H,T) and the 2-by-2 blocks correspond to complex conjugate pairs of 00044 * eigenvalues. 00045 * 00046 * Additionally, the 2-by-2 upper triangular diagonal blocks of P 00047 * corresponding to 2-by-2 blocks of S are reduced to positive diagonal 00048 * form, i.e., if S(j+1,j) is non-zero, then P(j+1,j) = P(j,j+1) = 0, 00049 * P(j,j) > 0, and P(j+1,j+1) > 0. 00050 * 00051 * Optionally, the orthogonal matrix Q from the generalized Schur 00052 * factorization may be postmultiplied into an input matrix Q1, and the 00053 * orthogonal matrix Z may be postmultiplied into an input matrix Z1. 00054 * If Q1 and Z1 are the orthogonal matrices from DGGHRD that reduced 00055 * the matrix pair (A,B) to generalized upper Hessenberg form, then the 00056 * output matrices Q1*Q and Z1*Z are the orthogonal factors from the 00057 * generalized Schur factorization of (A,B): 00058 * 00059 * A = (Q1*Q)*S*(Z1*Z)**T, B = (Q1*Q)*P*(Z1*Z)**T. 00060 * 00061 * To avoid overflow, eigenvalues of the matrix pair (H,T) (equivalently, 00062 * of (A,B)) are computed as a pair of values (alpha,beta), where alpha is 00063 * complex and beta real. 00064 * If beta is nonzero, lambda = alpha / beta is an eigenvalue of the 00065 * generalized nonsymmetric eigenvalue problem (GNEP) 00066 * A*x = lambda*B*x 00067 * and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the 00068 * alternate form of the GNEP 00069 * mu*A*y = B*y. 00070 * Real eigenvalues can be read directly from the generalized Schur 00071 * form: 00072 * alpha = S(i,i), beta = P(i,i). 00073 * 00074 * Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix 00075 * Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973), 00076 * pp. 241--256. 00077 * 00078 * Arguments 00079 * ========= 00080 * 00081 * JOB (input) CHARACTER*1 00082 * = 'E': Compute eigenvalues only; 00083 * = 'S': Compute eigenvalues and the Schur form. 00084 * 00085 * COMPQ (input) CHARACTER*1 00086 * = 'N': Left Schur vectors (Q) are not computed; 00087 * = 'I': Q is initialized to the unit matrix and the matrix Q 00088 * of left Schur vectors of (H,T) is returned; 00089 * = 'V': Q must contain an orthogonal matrix Q1 on entry and 00090 * the product Q1*Q is returned. 00091 * 00092 * COMPZ (input) CHARACTER*1 00093 * = 'N': Right Schur vectors (Z) are not computed; 00094 * = 'I': Z is initialized to the unit matrix and the matrix Z 00095 * of right Schur vectors of (H,T) is returned; 00096 * = 'V': Z must contain an orthogonal matrix Z1 on entry and 00097 * the product Z1*Z is returned. 00098 * 00099 * N (input) INTEGER 00100 * The order of the matrices H, T, Q, and Z. N >= 0. 00101 * 00102 * ILO (input) INTEGER 00103 * IHI (input) INTEGER 00104 * ILO and IHI mark the rows and columns of H which are in 00105 * Hessenberg form. It is assumed that A is already upper 00106 * triangular in rows and columns 1:ILO-1 and IHI+1:N. 00107 * If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0. 00108 * 00109 * H (input/output) DOUBLE PRECISION array, dimension (LDH, N) 00110 * On entry, the N-by-N upper Hessenberg matrix H. 00111 * On exit, if JOB = 'S', H contains the upper quasi-triangular 00112 * matrix S from the generalized Schur factorization; 00113 * 2-by-2 diagonal blocks (corresponding to complex conjugate 00114 * pairs of eigenvalues) are returned in standard form, with 00115 * H(i,i) = H(i+1,i+1) and H(i+1,i)*H(i,i+1) < 0. 00116 * If JOB = 'E', the diagonal blocks of H match those of S, but 00117 * the rest of H is unspecified. 00118 * 00119 * LDH (input) INTEGER 00120 * The leading dimension of the array H. LDH >= max( 1, N ). 00121 * 00122 * T (input/output) DOUBLE PRECISION array, dimension (LDT, N) 00123 * On entry, the N-by-N upper triangular matrix T. 00124 * On exit, if JOB = 'S', T contains the upper triangular 00125 * matrix P from the generalized Schur factorization; 00126 * 2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks of S 00127 * are reduced to positive diagonal form, i.e., if H(j+1,j) is 00128 * non-zero, then T(j+1,j) = T(j,j+1) = 0, T(j,j) > 0, and 00129 * T(j+1,j+1) > 0. 00130 * If JOB = 'E', the diagonal blocks of T match those of P, but 00131 * the rest of T is unspecified. 00132 * 00133 * LDT (input) INTEGER 00134 * The leading dimension of the array T. LDT >= max( 1, N ). 00135 * 00136 * ALPHAR (output) DOUBLE PRECISION array, dimension (N) 00137 * The real parts of each scalar alpha defining an eigenvalue 00138 * of GNEP. 00139 * 00140 * ALPHAI (output) DOUBLE PRECISION array, dimension (N) 00141 * The imaginary parts of each scalar alpha defining an 00142 * eigenvalue of GNEP. 00143 * If ALPHAI(j) is zero, then the j-th eigenvalue is real; if 00144 * positive, then the j-th and (j+1)-st eigenvalues are a 00145 * complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j). 00146 * 00147 * BETA (output) DOUBLE PRECISION array, dimension (N) 00148 * The scalars beta that define the eigenvalues of GNEP. 00149 * Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and 00150 * beta = BETA(j) represent the j-th eigenvalue of the matrix 00151 * pair (A,B), in one of the forms lambda = alpha/beta or 00152 * mu = beta/alpha. Since either lambda or mu may overflow, 00153 * they should not, in general, be computed. 00154 * 00155 * Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N) 00156 * On entry, if COMPZ = 'V', the orthogonal matrix Q1 used in 00157 * the reduction of (A,B) to generalized Hessenberg form. 00158 * On exit, if COMPZ = 'I', the orthogonal matrix of left Schur 00159 * vectors of (H,T), and if COMPZ = 'V', the orthogonal matrix 00160 * of left Schur vectors of (A,B). 00161 * Not referenced if COMPZ = 'N'. 00162 * 00163 * LDQ (input) INTEGER 00164 * The leading dimension of the array Q. LDQ >= 1. 00165 * If COMPQ='V' or 'I', then LDQ >= N. 00166 * 00167 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N) 00168 * On entry, if COMPZ = 'V', the orthogonal matrix Z1 used in 00169 * the reduction of (A,B) to generalized Hessenberg form. 00170 * On exit, if COMPZ = 'I', the orthogonal matrix of 00171 * right Schur vectors of (H,T), and if COMPZ = 'V', the 00172 * orthogonal matrix of right Schur vectors of (A,B). 00173 * Not referenced if COMPZ = 'N'. 00174 * 00175 * LDZ (input) INTEGER 00176 * The leading dimension of the array Z. LDZ >= 1. 00177 * If COMPZ='V' or 'I', then LDZ >= N. 00178 * 00179 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00180 * On exit, if INFO >= 0, WORK(1) returns the optimal LWORK. 00181 * 00182 * LWORK (input) INTEGER 00183 * The dimension of the array WORK. LWORK >= max(1,N). 00184 * 00185 * If LWORK = -1, then a workspace query is assumed; the routine 00186 * only calculates the optimal size of the WORK array, returns 00187 * this value as the first entry of the WORK array, and no error 00188 * message related to LWORK is issued by XERBLA. 00189 * 00190 * INFO (output) INTEGER 00191 * = 0: successful exit 00192 * < 0: if INFO = -i, the i-th argument had an illegal value 00193 * = 1,...,N: the QZ iteration did not converge. (H,T) is not 00194 * in Schur form, but ALPHAR(i), ALPHAI(i), and 00195 * BETA(i), i=INFO+1,...,N should be correct. 00196 * = N+1,...,2*N: the shift calculation failed. (H,T) is not 00197 * in Schur form, but ALPHAR(i), ALPHAI(i), and 00198 * BETA(i), i=INFO-N+1,...,N should be correct. 00199 * 00200 * Further Details 00201 * =============== 00202 * 00203 * Iteration counters: 00204 * 00205 * JITER -- counts iterations. 00206 * IITER -- counts iterations run since ILAST was last 00207 * changed. This is therefore reset only when a 1-by-1 or 00208 * 2-by-2 block deflates off the bottom. 00209 * 00210 * ===================================================================== 00211 * 00212 * .. Parameters .. 00213 * $ SAFETY = 1.0E+0 ) 00214 DOUBLE PRECISION HALF, ZERO, ONE, SAFETY 00215 PARAMETER ( HALF = 0.5D+0, ZERO = 0.0D+0, ONE = 1.0D+0, 00216 $ SAFETY = 1.0D+2 ) 00217 * .. 00218 * .. Local Scalars .. 00219 LOGICAL ILAZR2, ILAZRO, ILPIVT, ILQ, ILSCHR, ILZ, 00220 $ LQUERY 00221 INTEGER ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST, 00222 $ ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER, 00223 $ JR, MAXIT 00224 DOUBLE PRECISION A11, A12, A1I, A1R, A21, A22, A2I, A2R, AD11, 00225 $ AD11L, AD12, AD12L, AD21, AD21L, AD22, AD22L, 00226 $ AD32L, AN, ANORM, ASCALE, ATOL, B11, B1A, B1I, 00227 $ B1R, B22, B2A, B2I, B2R, BN, BNORM, BSCALE, 00228 $ BTOL, C, C11I, C11R, C12, C21, C22I, C22R, CL, 00229 $ CQ, CR, CZ, ESHIFT, S, S1, S1INV, S2, SAFMAX, 00230 $ SAFMIN, SCALE, SL, SQI, SQR, SR, SZI, SZR, T1, 00231 $ TAU, TEMP, TEMP2, TEMPI, TEMPR, U1, U12, U12L, 00232 $ U2, ULP, VS, W11, W12, W21, W22, WABS, WI, WR, 00233 $ WR2 00234 * .. 00235 * .. Local Arrays .. 00236 DOUBLE PRECISION V( 3 ) 00237 * .. 00238 * .. External Functions .. 00239 LOGICAL LSAME 00240 DOUBLE PRECISION DLAMCH, DLANHS, DLAPY2, DLAPY3 00241 EXTERNAL LSAME, DLAMCH, DLANHS, DLAPY2, DLAPY3 00242 * .. 00243 * .. External Subroutines .. 00244 EXTERNAL DLAG2, DLARFG, DLARTG, DLASET, DLASV2, DROT, 00245 $ XERBLA 00246 * .. 00247 * .. Intrinsic Functions .. 00248 INTRINSIC ABS, DBLE, MAX, MIN, SQRT 00249 * .. 00250 * .. Executable Statements .. 00251 * 00252 * Decode JOB, COMPQ, COMPZ 00253 * 00254 IF( LSAME( JOB, 'E' ) ) THEN 00255 ILSCHR = .FALSE. 00256 ISCHUR = 1 00257 ELSE IF( LSAME( JOB, 'S' ) ) THEN 00258 ILSCHR = .TRUE. 00259 ISCHUR = 2 00260 ELSE 00261 ISCHUR = 0 00262 END IF 00263 * 00264 IF( LSAME( COMPQ, 'N' ) ) THEN 00265 ILQ = .FALSE. 00266 ICOMPQ = 1 00267 ELSE IF( LSAME( COMPQ, 'V' ) ) THEN 00268 ILQ = .TRUE. 00269 ICOMPQ = 2 00270 ELSE IF( LSAME( COMPQ, 'I' ) ) THEN 00271 ILQ = .TRUE. 00272 ICOMPQ = 3 00273 ELSE 00274 ICOMPQ = 0 00275 END IF 00276 * 00277 IF( LSAME( COMPZ, 'N' ) ) THEN 00278 ILZ = .FALSE. 00279 ICOMPZ = 1 00280 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN 00281 ILZ = .TRUE. 00282 ICOMPZ = 2 00283 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN 00284 ILZ = .TRUE. 00285 ICOMPZ = 3 00286 ELSE 00287 ICOMPZ = 0 00288 END IF 00289 * 00290 * Check Argument Values 00291 * 00292 INFO = 0 00293 WORK( 1 ) = MAX( 1, N ) 00294 LQUERY = ( LWORK.EQ.-1 ) 00295 IF( ISCHUR.EQ.0 ) THEN 00296 INFO = -1 00297 ELSE IF( ICOMPQ.EQ.0 ) THEN 00298 INFO = -2 00299 ELSE IF( ICOMPZ.EQ.0 ) THEN 00300 INFO = -3 00301 ELSE IF( N.LT.0 ) THEN 00302 INFO = -4 00303 ELSE IF( ILO.LT.1 ) THEN 00304 INFO = -5 00305 ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN 00306 INFO = -6 00307 ELSE IF( LDH.LT.N ) THEN 00308 INFO = -8 00309 ELSE IF( LDT.LT.N ) THEN 00310 INFO = -10 00311 ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN 00312 INFO = -15 00313 ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN 00314 INFO = -17 00315 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN 00316 INFO = -19 00317 END IF 00318 IF( INFO.NE.0 ) THEN 00319 CALL XERBLA( 'DHGEQZ', -INFO ) 00320 RETURN 00321 ELSE IF( LQUERY ) THEN 00322 RETURN 00323 END IF 00324 * 00325 * Quick return if possible 00326 * 00327 IF( N.LE.0 ) THEN 00328 WORK( 1 ) = DBLE( 1 ) 00329 RETURN 00330 END IF 00331 * 00332 * Initialize Q and Z 00333 * 00334 IF( ICOMPQ.EQ.3 ) 00335 $ CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ ) 00336 IF( ICOMPZ.EQ.3 ) 00337 $ CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ ) 00338 * 00339 * Machine Constants 00340 * 00341 IN = IHI + 1 - ILO 00342 SAFMIN = DLAMCH( 'S' ) 00343 SAFMAX = ONE / SAFMIN 00344 ULP = DLAMCH( 'E' )*DLAMCH( 'B' ) 00345 ANORM = DLANHS( 'F', IN, H( ILO, ILO ), LDH, WORK ) 00346 BNORM = DLANHS( 'F', IN, T( ILO, ILO ), LDT, WORK ) 00347 ATOL = MAX( SAFMIN, ULP*ANORM ) 00348 BTOL = MAX( SAFMIN, ULP*BNORM ) 00349 ASCALE = ONE / MAX( SAFMIN, ANORM ) 00350 BSCALE = ONE / MAX( SAFMIN, BNORM ) 00351 * 00352 * Set Eigenvalues IHI+1:N 00353 * 00354 DO 30 J = IHI + 1, N 00355 IF( T( J, J ).LT.ZERO ) THEN 00356 IF( ILSCHR ) THEN 00357 DO 10 JR = 1, J 00358 H( JR, J ) = -H( JR, J ) 00359 T( JR, J ) = -T( JR, J ) 00360 10 CONTINUE 00361 ELSE 00362 H( J, J ) = -H( J, J ) 00363 T( J, J ) = -T( J, J ) 00364 END IF 00365 IF( ILZ ) THEN 00366 DO 20 JR = 1, N 00367 Z( JR, J ) = -Z( JR, J ) 00368 20 CONTINUE 00369 END IF 00370 END IF 00371 ALPHAR( J ) = H( J, J ) 00372 ALPHAI( J ) = ZERO 00373 BETA( J ) = T( J, J ) 00374 30 CONTINUE 00375 * 00376 * If IHI < ILO, skip QZ steps 00377 * 00378 IF( IHI.LT.ILO ) 00379 $ GO TO 380 00380 * 00381 * MAIN QZ ITERATION LOOP 00382 * 00383 * Initialize dynamic indices 00384 * 00385 * Eigenvalues ILAST+1:N have been found. 00386 * Column operations modify rows IFRSTM:whatever. 00387 * Row operations modify columns whatever:ILASTM. 00388 * 00389 * If only eigenvalues are being computed, then 00390 * IFRSTM is the row of the last splitting row above row ILAST; 00391 * this is always at least ILO. 00392 * IITER counts iterations since the last eigenvalue was found, 00393 * to tell when to use an extraordinary shift. 00394 * MAXIT is the maximum number of QZ sweeps allowed. 00395 * 00396 ILAST = IHI 00397 IF( ILSCHR ) THEN 00398 IFRSTM = 1 00399 ILASTM = N 00400 ELSE 00401 IFRSTM = ILO 00402 ILASTM = IHI 00403 END IF 00404 IITER = 0 00405 ESHIFT = ZERO 00406 MAXIT = 30*( IHI-ILO+1 ) 00407 * 00408 DO 360 JITER = 1, MAXIT 00409 * 00410 * Split the matrix if possible. 00411 * 00412 * Two tests: 00413 * 1: H(j,j-1)=0 or j=ILO 00414 * 2: T(j,j)=0 00415 * 00416 IF( ILAST.EQ.ILO ) THEN 00417 * 00418 * Special case: j=ILAST 00419 * 00420 GO TO 80 00421 ELSE 00422 IF( ABS( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN 00423 H( ILAST, ILAST-1 ) = ZERO 00424 GO TO 80 00425 END IF 00426 END IF 00427 * 00428 IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN 00429 T( ILAST, ILAST ) = ZERO 00430 GO TO 70 00431 END IF 00432 * 00433 * General case: j<ILAST 00434 * 00435 DO 60 J = ILAST - 1, ILO, -1 00436 * 00437 * Test 1: for H(j,j-1)=0 or j=ILO 00438 * 00439 IF( J.EQ.ILO ) THEN 00440 ILAZRO = .TRUE. 00441 ELSE 00442 IF( ABS( H( J, J-1 ) ).LE.ATOL ) THEN 00443 H( J, J-1 ) = ZERO 00444 ILAZRO = .TRUE. 00445 ELSE 00446 ILAZRO = .FALSE. 00447 END IF 00448 END IF 00449 * 00450 * Test 2: for T(j,j)=0 00451 * 00452 IF( ABS( T( J, J ) ).LT.BTOL ) THEN 00453 T( J, J ) = ZERO 00454 * 00455 * Test 1a: Check for 2 consecutive small subdiagonals in A 00456 * 00457 ILAZR2 = .FALSE. 00458 IF( .NOT.ILAZRO ) THEN 00459 TEMP = ABS( H( J, J-1 ) ) 00460 TEMP2 = ABS( H( J, J ) ) 00461 TEMPR = MAX( TEMP, TEMP2 ) 00462 IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN 00463 TEMP = TEMP / TEMPR 00464 TEMP2 = TEMP2 / TEMPR 00465 END IF 00466 IF( TEMP*( ASCALE*ABS( H( J+1, J ) ) ).LE.TEMP2* 00467 $ ( ASCALE*ATOL ) )ILAZR2 = .TRUE. 00468 END IF 00469 * 00470 * If both tests pass (1 & 2), i.e., the leading diagonal 00471 * element of B in the block is zero, split a 1x1 block off 00472 * at the top. (I.e., at the J-th row/column) The leading 00473 * diagonal element of the remainder can also be zero, so 00474 * this may have to be done repeatedly. 00475 * 00476 IF( ILAZRO .OR. ILAZR2 ) THEN 00477 DO 40 JCH = J, ILAST - 1 00478 TEMP = H( JCH, JCH ) 00479 CALL DLARTG( TEMP, H( JCH+1, JCH ), C, S, 00480 $ H( JCH, JCH ) ) 00481 H( JCH+1, JCH ) = ZERO 00482 CALL DROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH, 00483 $ H( JCH+1, JCH+1 ), LDH, C, S ) 00484 CALL DROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT, 00485 $ T( JCH+1, JCH+1 ), LDT, C, S ) 00486 IF( ILQ ) 00487 $ CALL DROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1, 00488 $ C, S ) 00489 IF( ILAZR2 ) 00490 $ H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C 00491 ILAZR2 = .FALSE. 00492 IF( ABS( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN 00493 IF( JCH+1.GE.ILAST ) THEN 00494 GO TO 80 00495 ELSE 00496 IFIRST = JCH + 1 00497 GO TO 110 00498 END IF 00499 END IF 00500 T( JCH+1, JCH+1 ) = ZERO 00501 40 CONTINUE 00502 GO TO 70 00503 ELSE 00504 * 00505 * Only test 2 passed -- chase the zero to T(ILAST,ILAST) 00506 * Then process as in the case T(ILAST,ILAST)=0 00507 * 00508 DO 50 JCH = J, ILAST - 1 00509 TEMP = T( JCH, JCH+1 ) 00510 CALL DLARTG( TEMP, T( JCH+1, JCH+1 ), C, S, 00511 $ T( JCH, JCH+1 ) ) 00512 T( JCH+1, JCH+1 ) = ZERO 00513 IF( JCH.LT.ILASTM-1 ) 00514 $ CALL DROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT, 00515 $ T( JCH+1, JCH+2 ), LDT, C, S ) 00516 CALL DROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH, 00517 $ H( JCH+1, JCH-1 ), LDH, C, S ) 00518 IF( ILQ ) 00519 $ CALL DROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1, 00520 $ C, S ) 00521 TEMP = H( JCH+1, JCH ) 00522 CALL DLARTG( TEMP, H( JCH+1, JCH-1 ), C, S, 00523 $ H( JCH+1, JCH ) ) 00524 H( JCH+1, JCH-1 ) = ZERO 00525 CALL DROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1, 00526 $ H( IFRSTM, JCH-1 ), 1, C, S ) 00527 CALL DROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1, 00528 $ T( IFRSTM, JCH-1 ), 1, C, S ) 00529 IF( ILZ ) 00530 $ CALL DROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1, 00531 $ C, S ) 00532 50 CONTINUE 00533 GO TO 70 00534 END IF 00535 ELSE IF( ILAZRO ) THEN 00536 * 00537 * Only test 1 passed -- work on J:ILAST 00538 * 00539 IFIRST = J 00540 GO TO 110 00541 END IF 00542 * 00543 * Neither test passed -- try next J 00544 * 00545 60 CONTINUE 00546 * 00547 * (Drop-through is "impossible") 00548 * 00549 INFO = N + 1 00550 GO TO 420 00551 * 00552 * T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a 00553 * 1x1 block. 00554 * 00555 70 CONTINUE 00556 TEMP = H( ILAST, ILAST ) 00557 CALL DLARTG( TEMP, H( ILAST, ILAST-1 ), C, S, 00558 $ H( ILAST, ILAST ) ) 00559 H( ILAST, ILAST-1 ) = ZERO 00560 CALL DROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1, 00561 $ H( IFRSTM, ILAST-1 ), 1, C, S ) 00562 CALL DROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1, 00563 $ T( IFRSTM, ILAST-1 ), 1, C, S ) 00564 IF( ILZ ) 00565 $ CALL DROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S ) 00566 * 00567 * H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI, 00568 * and BETA 00569 * 00570 80 CONTINUE 00571 IF( T( ILAST, ILAST ).LT.ZERO ) THEN 00572 IF( ILSCHR ) THEN 00573 DO 90 J = IFRSTM, ILAST 00574 H( J, ILAST ) = -H( J, ILAST ) 00575 T( J, ILAST ) = -T( J, ILAST ) 00576 90 CONTINUE 00577 ELSE 00578 H( ILAST, ILAST ) = -H( ILAST, ILAST ) 00579 T( ILAST, ILAST ) = -T( ILAST, ILAST ) 00580 END IF 00581 IF( ILZ ) THEN 00582 DO 100 J = 1, N 00583 Z( J, ILAST ) = -Z( J, ILAST ) 00584 100 CONTINUE 00585 END IF 00586 END IF 00587 ALPHAR( ILAST ) = H( ILAST, ILAST ) 00588 ALPHAI( ILAST ) = ZERO 00589 BETA( ILAST ) = T( ILAST, ILAST ) 00590 * 00591 * Go to next block -- exit if finished. 00592 * 00593 ILAST = ILAST - 1 00594 IF( ILAST.LT.ILO ) 00595 $ GO TO 380 00596 * 00597 * Reset counters 00598 * 00599 IITER = 0 00600 ESHIFT = ZERO 00601 IF( .NOT.ILSCHR ) THEN 00602 ILASTM = ILAST 00603 IF( IFRSTM.GT.ILAST ) 00604 $ IFRSTM = ILO 00605 END IF 00606 GO TO 350 00607 * 00608 * QZ step 00609 * 00610 * This iteration only involves rows/columns IFIRST:ILAST. We 00611 * assume IFIRST < ILAST, and that the diagonal of B is non-zero. 00612 * 00613 110 CONTINUE 00614 IITER = IITER + 1 00615 IF( .NOT.ILSCHR ) THEN 00616 IFRSTM = IFIRST 00617 END IF 00618 * 00619 * Compute single shifts. 00620 * 00621 * At this point, IFIRST < ILAST, and the diagonal elements of 00622 * T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in 00623 * magnitude) 00624 * 00625 IF( ( IITER / 10 )*10.EQ.IITER ) THEN 00626 * 00627 * Exceptional shift. Chosen for no particularly good reason. 00628 * (Single shift only.) 00629 * 00630 IF( ( DBLE( MAXIT )*SAFMIN )*ABS( H( ILAST-1, ILAST ) ).LT. 00631 $ ABS( T( ILAST-1, ILAST-1 ) ) ) THEN 00632 ESHIFT = ESHIFT + H( ILAST-1, ILAST ) / 00633 $ T( ILAST-1, ILAST-1 ) 00634 ELSE 00635 ESHIFT = ESHIFT + ONE / ( SAFMIN*DBLE( MAXIT ) ) 00636 END IF 00637 S1 = ONE 00638 WR = ESHIFT 00639 * 00640 ELSE 00641 * 00642 * Shifts based on the generalized eigenvalues of the 00643 * bottom-right 2x2 block of A and B. The first eigenvalue 00644 * returned by DLAG2 is the Wilkinson shift (AEP p.512), 00645 * 00646 CALL DLAG2( H( ILAST-1, ILAST-1 ), LDH, 00647 $ T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1, 00648 $ S2, WR, WR2, WI ) 00649 * 00650 TEMP = MAX( S1, SAFMIN*MAX( ONE, ABS( WR ), ABS( WI ) ) ) 00651 IF( WI.NE.ZERO ) 00652 $ GO TO 200 00653 END IF 00654 * 00655 * Fiddle with shift to avoid overflow 00656 * 00657 TEMP = MIN( ASCALE, ONE )*( HALF*SAFMAX ) 00658 IF( S1.GT.TEMP ) THEN 00659 SCALE = TEMP / S1 00660 ELSE 00661 SCALE = ONE 00662 END IF 00663 * 00664 TEMP = MIN( BSCALE, ONE )*( HALF*SAFMAX ) 00665 IF( ABS( WR ).GT.TEMP ) 00666 $ SCALE = MIN( SCALE, TEMP / ABS( WR ) ) 00667 S1 = SCALE*S1 00668 WR = SCALE*WR 00669 * 00670 * Now check for two consecutive small subdiagonals. 00671 * 00672 DO 120 J = ILAST - 1, IFIRST + 1, -1 00673 ISTART = J 00674 TEMP = ABS( S1*H( J, J-1 ) ) 00675 TEMP2 = ABS( S1*H( J, J )-WR*T( J, J ) ) 00676 TEMPR = MAX( TEMP, TEMP2 ) 00677 IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN 00678 TEMP = TEMP / TEMPR 00679 TEMP2 = TEMP2 / TEMPR 00680 END IF 00681 IF( ABS( ( ASCALE*H( J+1, J ) )*TEMP ).LE.( ASCALE*ATOL )* 00682 $ TEMP2 )GO TO 130 00683 120 CONTINUE 00684 * 00685 ISTART = IFIRST 00686 130 CONTINUE 00687 * 00688 * Do an implicit single-shift QZ sweep. 00689 * 00690 * Initial Q 00691 * 00692 TEMP = S1*H( ISTART, ISTART ) - WR*T( ISTART, ISTART ) 00693 TEMP2 = S1*H( ISTART+1, ISTART ) 00694 CALL DLARTG( TEMP, TEMP2, C, S, TEMPR ) 00695 * 00696 * Sweep 00697 * 00698 DO 190 J = ISTART, ILAST - 1 00699 IF( J.GT.ISTART ) THEN 00700 TEMP = H( J, J-1 ) 00701 CALL DLARTG( TEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) ) 00702 H( J+1, J-1 ) = ZERO 00703 END IF 00704 * 00705 DO 140 JC = J, ILASTM 00706 TEMP = C*H( J, JC ) + S*H( J+1, JC ) 00707 H( J+1, JC ) = -S*H( J, JC ) + C*H( J+1, JC ) 00708 H( J, JC ) = TEMP 00709 TEMP2 = C*T( J, JC ) + S*T( J+1, JC ) 00710 T( J+1, JC ) = -S*T( J, JC ) + C*T( J+1, JC ) 00711 T( J, JC ) = TEMP2 00712 140 CONTINUE 00713 IF( ILQ ) THEN 00714 DO 150 JR = 1, N 00715 TEMP = C*Q( JR, J ) + S*Q( JR, J+1 ) 00716 Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 ) 00717 Q( JR, J ) = TEMP 00718 150 CONTINUE 00719 END IF 00720 * 00721 TEMP = T( J+1, J+1 ) 00722 CALL DLARTG( TEMP, T( J+1, J ), C, S, T( J+1, J+1 ) ) 00723 T( J+1, J ) = ZERO 00724 * 00725 DO 160 JR = IFRSTM, MIN( J+2, ILAST ) 00726 TEMP = C*H( JR, J+1 ) + S*H( JR, J ) 00727 H( JR, J ) = -S*H( JR, J+1 ) + C*H( JR, J ) 00728 H( JR, J+1 ) = TEMP 00729 160 CONTINUE 00730 DO 170 JR = IFRSTM, J 00731 TEMP = C*T( JR, J+1 ) + S*T( JR, J ) 00732 T( JR, J ) = -S*T( JR, J+1 ) + C*T( JR, J ) 00733 T( JR, J+1 ) = TEMP 00734 170 CONTINUE 00735 IF( ILZ ) THEN 00736 DO 180 JR = 1, N 00737 TEMP = C*Z( JR, J+1 ) + S*Z( JR, J ) 00738 Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J ) 00739 Z( JR, J+1 ) = TEMP 00740 180 CONTINUE 00741 END IF 00742 190 CONTINUE 00743 * 00744 GO TO 350 00745 * 00746 * Use Francis double-shift 00747 * 00748 * Note: the Francis double-shift should work with real shifts, 00749 * but only if the block is at least 3x3. 00750 * This code may break if this point is reached with 00751 * a 2x2 block with real eigenvalues. 00752 * 00753 200 CONTINUE 00754 IF( IFIRST+1.EQ.ILAST ) THEN 00755 * 00756 * Special case -- 2x2 block with complex eigenvectors 00757 * 00758 * Step 1: Standardize, that is, rotate so that 00759 * 00760 * ( B11 0 ) 00761 * B = ( ) with B11 non-negative. 00762 * ( 0 B22 ) 00763 * 00764 CALL DLASV2( T( ILAST-1, ILAST-1 ), T( ILAST-1, ILAST ), 00765 $ T( ILAST, ILAST ), B22, B11, SR, CR, SL, CL ) 00766 * 00767 IF( B11.LT.ZERO ) THEN 00768 CR = -CR 00769 SR = -SR 00770 B11 = -B11 00771 B22 = -B22 00772 END IF 00773 * 00774 CALL DROT( ILASTM+1-IFIRST, H( ILAST-1, ILAST-1 ), LDH, 00775 $ H( ILAST, ILAST-1 ), LDH, CL, SL ) 00776 CALL DROT( ILAST+1-IFRSTM, H( IFRSTM, ILAST-1 ), 1, 00777 $ H( IFRSTM, ILAST ), 1, CR, SR ) 00778 * 00779 IF( ILAST.LT.ILASTM ) 00780 $ CALL DROT( ILASTM-ILAST, T( ILAST-1, ILAST+1 ), LDT, 00781 $ T( ILAST, ILAST+1 ), LDT, CL, SL ) 00782 IF( IFRSTM.LT.ILAST-1 ) 00783 $ CALL DROT( IFIRST-IFRSTM, T( IFRSTM, ILAST-1 ), 1, 00784 $ T( IFRSTM, ILAST ), 1, CR, SR ) 00785 * 00786 IF( ILQ ) 00787 $ CALL DROT( N, Q( 1, ILAST-1 ), 1, Q( 1, ILAST ), 1, CL, 00788 $ SL ) 00789 IF( ILZ ) 00790 $ CALL DROT( N, Z( 1, ILAST-1 ), 1, Z( 1, ILAST ), 1, CR, 00791 $ SR ) 00792 * 00793 T( ILAST-1, ILAST-1 ) = B11 00794 T( ILAST-1, ILAST ) = ZERO 00795 T( ILAST, ILAST-1 ) = ZERO 00796 T( ILAST, ILAST ) = B22 00797 * 00798 * If B22 is negative, negate column ILAST 00799 * 00800 IF( B22.LT.ZERO ) THEN 00801 DO 210 J = IFRSTM, ILAST 00802 H( J, ILAST ) = -H( J, ILAST ) 00803 T( J, ILAST ) = -T( J, ILAST ) 00804 210 CONTINUE 00805 * 00806 IF( ILZ ) THEN 00807 DO 220 J = 1, N 00808 Z( J, ILAST ) = -Z( J, ILAST ) 00809 220 CONTINUE 00810 END IF 00811 END IF 00812 * 00813 * Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.) 00814 * 00815 * Recompute shift 00816 * 00817 CALL DLAG2( H( ILAST-1, ILAST-1 ), LDH, 00818 $ T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1, 00819 $ TEMP, WR, TEMP2, WI ) 00820 * 00821 * If standardization has perturbed the shift onto real line, 00822 * do another (real single-shift) QR step. 00823 * 00824 IF( WI.EQ.ZERO ) 00825 $ GO TO 350 00826 S1INV = ONE / S1 00827 * 00828 * Do EISPACK (QZVAL) computation of alpha and beta 00829 * 00830 A11 = H( ILAST-1, ILAST-1 ) 00831 A21 = H( ILAST, ILAST-1 ) 00832 A12 = H( ILAST-1, ILAST ) 00833 A22 = H( ILAST, ILAST ) 00834 * 00835 * Compute complex Givens rotation on right 00836 * (Assume some element of C = (sA - wB) > unfl ) 00837 * __ 00838 * (sA - wB) ( CZ -SZ ) 00839 * ( SZ CZ ) 00840 * 00841 C11R = S1*A11 - WR*B11 00842 C11I = -WI*B11 00843 C12 = S1*A12 00844 C21 = S1*A21 00845 C22R = S1*A22 - WR*B22 00846 C22I = -WI*B22 00847 * 00848 IF( ABS( C11R )+ABS( C11I )+ABS( C12 ).GT.ABS( C21 )+ 00849 $ ABS( C22R )+ABS( C22I ) ) THEN 00850 T1 = DLAPY3( C12, C11R, C11I ) 00851 CZ = C12 / T1 00852 SZR = -C11R / T1 00853 SZI = -C11I / T1 00854 ELSE 00855 CZ = DLAPY2( C22R, C22I ) 00856 IF( CZ.LE.SAFMIN ) THEN 00857 CZ = ZERO 00858 SZR = ONE 00859 SZI = ZERO 00860 ELSE 00861 TEMPR = C22R / CZ 00862 TEMPI = C22I / CZ 00863 T1 = DLAPY2( CZ, C21 ) 00864 CZ = CZ / T1 00865 SZR = -C21*TEMPR / T1 00866 SZI = C21*TEMPI / T1 00867 END IF 00868 END IF 00869 * 00870 * Compute Givens rotation on left 00871 * 00872 * ( CQ SQ ) 00873 * ( __ ) A or B 00874 * ( -SQ CQ ) 00875 * 00876 AN = ABS( A11 ) + ABS( A12 ) + ABS( A21 ) + ABS( A22 ) 00877 BN = ABS( B11 ) + ABS( B22 ) 00878 WABS = ABS( WR ) + ABS( WI ) 00879 IF( S1*AN.GT.WABS*BN ) THEN 00880 CQ = CZ*B11 00881 SQR = SZR*B22 00882 SQI = -SZI*B22 00883 ELSE 00884 A1R = CZ*A11 + SZR*A12 00885 A1I = SZI*A12 00886 A2R = CZ*A21 + SZR*A22 00887 A2I = SZI*A22 00888 CQ = DLAPY2( A1R, A1I ) 00889 IF( CQ.LE.SAFMIN ) THEN 00890 CQ = ZERO 00891 SQR = ONE 00892 SQI = ZERO 00893 ELSE 00894 TEMPR = A1R / CQ 00895 TEMPI = A1I / CQ 00896 SQR = TEMPR*A2R + TEMPI*A2I 00897 SQI = TEMPI*A2R - TEMPR*A2I 00898 END IF 00899 END IF 00900 T1 = DLAPY3( CQ, SQR, SQI ) 00901 CQ = CQ / T1 00902 SQR = SQR / T1 00903 SQI = SQI / T1 00904 * 00905 * Compute diagonal elements of QBZ 00906 * 00907 TEMPR = SQR*SZR - SQI*SZI 00908 TEMPI = SQR*SZI + SQI*SZR 00909 B1R = CQ*CZ*B11 + TEMPR*B22 00910 B1I = TEMPI*B22 00911 B1A = DLAPY2( B1R, B1I ) 00912 B2R = CQ*CZ*B22 + TEMPR*B11 00913 B2I = -TEMPI*B11 00914 B2A = DLAPY2( B2R, B2I ) 00915 * 00916 * Normalize so beta > 0, and Im( alpha1 ) > 0 00917 * 00918 BETA( ILAST-1 ) = B1A 00919 BETA( ILAST ) = B2A 00920 ALPHAR( ILAST-1 ) = ( WR*B1A )*S1INV 00921 ALPHAI( ILAST-1 ) = ( WI*B1A )*S1INV 00922 ALPHAR( ILAST ) = ( WR*B2A )*S1INV 00923 ALPHAI( ILAST ) = -( WI*B2A )*S1INV 00924 * 00925 * Step 3: Go to next block -- exit if finished. 00926 * 00927 ILAST = IFIRST - 1 00928 IF( ILAST.LT.ILO ) 00929 $ GO TO 380 00930 * 00931 * Reset counters 00932 * 00933 IITER = 0 00934 ESHIFT = ZERO 00935 IF( .NOT.ILSCHR ) THEN 00936 ILASTM = ILAST 00937 IF( IFRSTM.GT.ILAST ) 00938 $ IFRSTM = ILO 00939 END IF 00940 GO TO 350 00941 ELSE 00942 * 00943 * Usual case: 3x3 or larger block, using Francis implicit 00944 * double-shift 00945 * 00946 * 2 00947 * Eigenvalue equation is w - c w + d = 0, 00948 * 00949 * -1 2 -1 00950 * so compute 1st column of (A B ) - c A B + d 00951 * using the formula in QZIT (from EISPACK) 00952 * 00953 * We assume that the block is at least 3x3 00954 * 00955 AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) / 00956 $ ( BSCALE*T( ILAST-1, ILAST-1 ) ) 00957 AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) / 00958 $ ( BSCALE*T( ILAST-1, ILAST-1 ) ) 00959 AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) / 00960 $ ( BSCALE*T( ILAST, ILAST ) ) 00961 AD22 = ( ASCALE*H( ILAST, ILAST ) ) / 00962 $ ( BSCALE*T( ILAST, ILAST ) ) 00963 U12 = T( ILAST-1, ILAST ) / T( ILAST, ILAST ) 00964 AD11L = ( ASCALE*H( IFIRST, IFIRST ) ) / 00965 $ ( BSCALE*T( IFIRST, IFIRST ) ) 00966 AD21L = ( ASCALE*H( IFIRST+1, IFIRST ) ) / 00967 $ ( BSCALE*T( IFIRST, IFIRST ) ) 00968 AD12L = ( ASCALE*H( IFIRST, IFIRST+1 ) ) / 00969 $ ( BSCALE*T( IFIRST+1, IFIRST+1 ) ) 00970 AD22L = ( ASCALE*H( IFIRST+1, IFIRST+1 ) ) / 00971 $ ( BSCALE*T( IFIRST+1, IFIRST+1 ) ) 00972 AD32L = ( ASCALE*H( IFIRST+2, IFIRST+1 ) ) / 00973 $ ( BSCALE*T( IFIRST+1, IFIRST+1 ) ) 00974 U12L = T( IFIRST, IFIRST+1 ) / T( IFIRST+1, IFIRST+1 ) 00975 * 00976 V( 1 ) = ( AD11-AD11L )*( AD22-AD11L ) - AD12*AD21 + 00977 $ AD21*U12*AD11L + ( AD12L-AD11L*U12L )*AD21L 00978 V( 2 ) = ( ( AD22L-AD11L )-AD21L*U12L-( AD11-AD11L )- 00979 $ ( AD22-AD11L )+AD21*U12 )*AD21L 00980 V( 3 ) = AD32L*AD21L 00981 * 00982 ISTART = IFIRST 00983 * 00984 CALL DLARFG( 3, V( 1 ), V( 2 ), 1, TAU ) 00985 V( 1 ) = ONE 00986 * 00987 * Sweep 00988 * 00989 DO 290 J = ISTART, ILAST - 2 00990 * 00991 * All but last elements: use 3x3 Householder transforms. 00992 * 00993 * Zero (j-1)st column of A 00994 * 00995 IF( J.GT.ISTART ) THEN 00996 V( 1 ) = H( J, J-1 ) 00997 V( 2 ) = H( J+1, J-1 ) 00998 V( 3 ) = H( J+2, J-1 ) 00999 * 01000 CALL DLARFG( 3, H( J, J-1 ), V( 2 ), 1, TAU ) 01001 V( 1 ) = ONE 01002 H( J+1, J-1 ) = ZERO 01003 H( J+2, J-1 ) = ZERO 01004 END IF 01005 * 01006 DO 230 JC = J, ILASTM 01007 TEMP = TAU*( H( J, JC )+V( 2 )*H( J+1, JC )+V( 3 )* 01008 $ H( J+2, JC ) ) 01009 H( J, JC ) = H( J, JC ) - TEMP 01010 H( J+1, JC ) = H( J+1, JC ) - TEMP*V( 2 ) 01011 H( J+2, JC ) = H( J+2, JC ) - TEMP*V( 3 ) 01012 TEMP2 = TAU*( T( J, JC )+V( 2 )*T( J+1, JC )+V( 3 )* 01013 $ T( J+2, JC ) ) 01014 T( J, JC ) = T( J, JC ) - TEMP2 01015 T( J+1, JC ) = T( J+1, JC ) - TEMP2*V( 2 ) 01016 T( J+2, JC ) = T( J+2, JC ) - TEMP2*V( 3 ) 01017 230 CONTINUE 01018 IF( ILQ ) THEN 01019 DO 240 JR = 1, N 01020 TEMP = TAU*( Q( JR, J )+V( 2 )*Q( JR, J+1 )+V( 3 )* 01021 $ Q( JR, J+2 ) ) 01022 Q( JR, J ) = Q( JR, J ) - TEMP 01023 Q( JR, J+1 ) = Q( JR, J+1 ) - TEMP*V( 2 ) 01024 Q( JR, J+2 ) = Q( JR, J+2 ) - TEMP*V( 3 ) 01025 240 CONTINUE 01026 END IF 01027 * 01028 * Zero j-th column of B (see DLAGBC for details) 01029 * 01030 * Swap rows to pivot 01031 * 01032 ILPIVT = .FALSE. 01033 TEMP = MAX( ABS( T( J+1, J+1 ) ), ABS( T( J+1, J+2 ) ) ) 01034 TEMP2 = MAX( ABS( T( J+2, J+1 ) ), ABS( T( J+2, J+2 ) ) ) 01035 IF( MAX( TEMP, TEMP2 ).LT.SAFMIN ) THEN 01036 SCALE = ZERO 01037 U1 = ONE 01038 U2 = ZERO 01039 GO TO 250 01040 ELSE IF( TEMP.GE.TEMP2 ) THEN 01041 W11 = T( J+1, J+1 ) 01042 W21 = T( J+2, J+1 ) 01043 W12 = T( J+1, J+2 ) 01044 W22 = T( J+2, J+2 ) 01045 U1 = T( J+1, J ) 01046 U2 = T( J+2, J ) 01047 ELSE 01048 W21 = T( J+1, J+1 ) 01049 W11 = T( J+2, J+1 ) 01050 W22 = T( J+1, J+2 ) 01051 W12 = T( J+2, J+2 ) 01052 U2 = T( J+1, J ) 01053 U1 = T( J+2, J ) 01054 END IF 01055 * 01056 * Swap columns if nec. 01057 * 01058 IF( ABS( W12 ).GT.ABS( W11 ) ) THEN 01059 ILPIVT = .TRUE. 01060 TEMP = W12 01061 TEMP2 = W22 01062 W12 = W11 01063 W22 = W21 01064 W11 = TEMP 01065 W21 = TEMP2 01066 END IF 01067 * 01068 * LU-factor 01069 * 01070 TEMP = W21 / W11 01071 U2 = U2 - TEMP*U1 01072 W22 = W22 - TEMP*W12 01073 W21 = ZERO 01074 * 01075 * Compute SCALE 01076 * 01077 SCALE = ONE 01078 IF( ABS( W22 ).LT.SAFMIN ) THEN 01079 SCALE = ZERO 01080 U2 = ONE 01081 U1 = -W12 / W11 01082 GO TO 250 01083 END IF 01084 IF( ABS( W22 ).LT.ABS( U2 ) ) 01085 $ SCALE = ABS( W22 / U2 ) 01086 IF( ABS( W11 ).LT.ABS( U1 ) ) 01087 $ SCALE = MIN( SCALE, ABS( W11 / U1 ) ) 01088 * 01089 * Solve 01090 * 01091 U2 = ( SCALE*U2 ) / W22 01092 U1 = ( SCALE*U1-W12*U2 ) / W11 01093 * 01094 250 CONTINUE 01095 IF( ILPIVT ) THEN 01096 TEMP = U2 01097 U2 = U1 01098 U1 = TEMP 01099 END IF 01100 * 01101 * Compute Householder Vector 01102 * 01103 T1 = SQRT( SCALE**2+U1**2+U2**2 ) 01104 TAU = ONE + SCALE / T1 01105 VS = -ONE / ( SCALE+T1 ) 01106 V( 1 ) = ONE 01107 V( 2 ) = VS*U1 01108 V( 3 ) = VS*U2 01109 * 01110 * Apply transformations from the right. 01111 * 01112 DO 260 JR = IFRSTM, MIN( J+3, ILAST ) 01113 TEMP = TAU*( H( JR, J )+V( 2 )*H( JR, J+1 )+V( 3 )* 01114 $ H( JR, J+2 ) ) 01115 H( JR, J ) = H( JR, J ) - TEMP 01116 H( JR, J+1 ) = H( JR, J+1 ) - TEMP*V( 2 ) 01117 H( JR, J+2 ) = H( JR, J+2 ) - TEMP*V( 3 ) 01118 260 CONTINUE 01119 DO 270 JR = IFRSTM, J + 2 01120 TEMP = TAU*( T( JR, J )+V( 2 )*T( JR, J+1 )+V( 3 )* 01121 $ T( JR, J+2 ) ) 01122 T( JR, J ) = T( JR, J ) - TEMP 01123 T( JR, J+1 ) = T( JR, J+1 ) - TEMP*V( 2 ) 01124 T( JR, J+2 ) = T( JR, J+2 ) - TEMP*V( 3 ) 01125 270 CONTINUE 01126 IF( ILZ ) THEN 01127 DO 280 JR = 1, N 01128 TEMP = TAU*( Z( JR, J )+V( 2 )*Z( JR, J+1 )+V( 3 )* 01129 $ Z( JR, J+2 ) ) 01130 Z( JR, J ) = Z( JR, J ) - TEMP 01131 Z( JR, J+1 ) = Z( JR, J+1 ) - TEMP*V( 2 ) 01132 Z( JR, J+2 ) = Z( JR, J+2 ) - TEMP*V( 3 ) 01133 280 CONTINUE 01134 END IF 01135 T( J+1, J ) = ZERO 01136 T( J+2, J ) = ZERO 01137 290 CONTINUE 01138 * 01139 * Last elements: Use Givens rotations 01140 * 01141 * Rotations from the left 01142 * 01143 J = ILAST - 1 01144 TEMP = H( J, J-1 ) 01145 CALL DLARTG( TEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) ) 01146 H( J+1, J-1 ) = ZERO 01147 * 01148 DO 300 JC = J, ILASTM 01149 TEMP = C*H( J, JC ) + S*H( J+1, JC ) 01150 H( J+1, JC ) = -S*H( J, JC ) + C*H( J+1, JC ) 01151 H( J, JC ) = TEMP 01152 TEMP2 = C*T( J, JC ) + S*T( J+1, JC ) 01153 T( J+1, JC ) = -S*T( J, JC ) + C*T( J+1, JC ) 01154 T( J, JC ) = TEMP2 01155 300 CONTINUE 01156 IF( ILQ ) THEN 01157 DO 310 JR = 1, N 01158 TEMP = C*Q( JR, J ) + S*Q( JR, J+1 ) 01159 Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 ) 01160 Q( JR, J ) = TEMP 01161 310 CONTINUE 01162 END IF 01163 * 01164 * Rotations from the right. 01165 * 01166 TEMP = T( J+1, J+1 ) 01167 CALL DLARTG( TEMP, T( J+1, J ), C, S, T( J+1, J+1 ) ) 01168 T( J+1, J ) = ZERO 01169 * 01170 DO 320 JR = IFRSTM, ILAST 01171 TEMP = C*H( JR, J+1 ) + S*H( JR, J ) 01172 H( JR, J ) = -S*H( JR, J+1 ) + C*H( JR, J ) 01173 H( JR, J+1 ) = TEMP 01174 320 CONTINUE 01175 DO 330 JR = IFRSTM, ILAST - 1 01176 TEMP = C*T( JR, J+1 ) + S*T( JR, J ) 01177 T( JR, J ) = -S*T( JR, J+1 ) + C*T( JR, J ) 01178 T( JR, J+1 ) = TEMP 01179 330 CONTINUE 01180 IF( ILZ ) THEN 01181 DO 340 JR = 1, N 01182 TEMP = C*Z( JR, J+1 ) + S*Z( JR, J ) 01183 Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J ) 01184 Z( JR, J+1 ) = TEMP 01185 340 CONTINUE 01186 END IF 01187 * 01188 * End of Double-Shift code 01189 * 01190 END IF 01191 * 01192 GO TO 350 01193 * 01194 * End of iteration loop 01195 * 01196 350 CONTINUE 01197 360 CONTINUE 01198 * 01199 * Drop-through = non-convergence 01200 * 01201 INFO = ILAST 01202 GO TO 420 01203 * 01204 * Successful completion of all QZ steps 01205 * 01206 380 CONTINUE 01207 * 01208 * Set Eigenvalues 1:ILO-1 01209 * 01210 DO 410 J = 1, ILO - 1 01211 IF( T( J, J ).LT.ZERO ) THEN 01212 IF( ILSCHR ) THEN 01213 DO 390 JR = 1, J 01214 H( JR, J ) = -H( JR, J ) 01215 T( JR, J ) = -T( JR, J ) 01216 390 CONTINUE 01217 ELSE 01218 H( J, J ) = -H( J, J ) 01219 T( J, J ) = -T( J, J ) 01220 END IF 01221 IF( ILZ ) THEN 01222 DO 400 JR = 1, N 01223 Z( JR, J ) = -Z( JR, J ) 01224 400 CONTINUE 01225 END IF 01226 END IF 01227 ALPHAR( J ) = H( J, J ) 01228 ALPHAI( J ) = ZERO 01229 BETA( J ) = T( J, J ) 01230 410 CONTINUE 01231 * 01232 * Normal Termination 01233 * 01234 INFO = 0 01235 * 01236 * Exit (other than argument error) -- return optimal workspace size 01237 * 01238 420 CONTINUE 01239 WORK( 1 ) = DBLE( N ) 01240 RETURN 01241 * 01242 * End of DHGEQZ 01243 * 01244 END
1.7.3 
