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