1 SUBROUTINE SHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, 2 $ ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, 3 $ LWORK, INFO ) 4* 5* -- LAPACK routine (instrumented to count operations, version 3.0) -- 6* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 7* Courant Institute, Argonne National Lab, and Rice University 8* June 30, 1999 9* 10* .. Scalar Arguments .. 11 CHARACTER COMPQ, COMPZ, JOB 12 INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, LWORK, N 13* .. 14* .. Array Arguments .. 15 REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ), 16 $ B( LDB, * ), BETA( * ), Q( LDQ, * ), WORK( * ), 17 $ Z( LDZ, * ) 18* .. 19* ---------------------- Begin Timing Code ------------------------- 20* Common block to return operation count and iteration count 21* ITCNT is initialized to 0, OPS is only incremented 22* OPST is used to accumulate small contributions to OPS 23* to avoid roundoff error 24* .. Common blocks .. 25 COMMON / LATIME / OPS, ITCNT 26* .. 27* .. Scalars in Common .. 28 REAL ITCNT, OPS 29* .. 30* ----------------------- End Timing Code -------------------------- 31* 32* Purpose 33* ======= 34* 35* SHGEQZ implements a single-/double-shift version of the QZ method for 36* finding the generalized eigenvalues 37* 38* w(j)=(ALPHAR(j) + i*ALPHAI(j))/BETAR(j) of the equation 39* 40* det( A - w(i) B ) = 0 41* 42* In addition, the pair A,B may be reduced to generalized Schur form: 43* B is upper triangular, and A is block upper triangular, where the 44* diagonal blocks are either 1-by-1 or 2-by-2, the 2-by-2 blocks having 45* complex generalized eigenvalues (see the description of the argument 46* JOB.) 47* 48* If JOB='S', then the pair (A,B) is simultaneously reduced to Schur 49* form by applying one orthogonal tranformation (usually called Q) on 50* the left and another (usually called Z) on the right. The 2-by-2 51* upper-triangular diagonal blocks of B corresponding to 2-by-2 blocks 52* of A will be reduced to positive diagonal matrices. (I.e., 53* if A(j+1,j) is non-zero, then B(j+1,j)=B(j,j+1)=0 and B(j,j) and 54* B(j+1,j+1) will be positive.) 55* 56* If JOB='E', then at each iteration, the same transformations 57* are computed, but they are only applied to those parts of A and B 58* which are needed to compute ALPHAR, ALPHAI, and BETAR. 59* 60* If JOB='S' and COMPQ and COMPZ are 'V' or 'I', then the orthogonal 61* transformations used to reduce (A,B) are accumulated into the arrays 62* Q and Z s.t.: 63* 64* Q(in) A(in) Z(in)* = Q(out) A(out) Z(out)* 65* Q(in) B(in) Z(in)* = Q(out) B(out) Z(out)* 66* 67* Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix 68* Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973), 69* pp. 241--256. 70* 71* Arguments 72* ========= 73* 74* JOB (input) CHARACTER*1 75* = 'E': compute only ALPHAR, ALPHAI, and BETA. A and B will 76* not necessarily be put into generalized Schur form. 77* = 'S': put A and B into generalized Schur form, as well 78* as computing ALPHAR, ALPHAI, and BETA. 79* 80* COMPQ (input) CHARACTER*1 81* = 'N': do not modify Q. 82* = 'V': multiply the array Q on the right by the transpose of 83* the orthogonal tranformation that is applied to the 84* left side of A and B to reduce them to Schur form. 85* = 'I': like COMPQ='V', except that Q will be initialized to 86* the identity first. 87* 88* COMPZ (input) CHARACTER*1 89* = 'N': do not modify Z. 90* = 'V': multiply the array Z on the right by the orthogonal 91* tranformation that is applied to the right side of 92* A and B to reduce them to Schur form. 93* = 'I': like COMPZ='V', except that Z will be initialized to 94* the identity first. 95* 96* N (input) INTEGER 97* The order of the matrices A, B, Q, and Z. N >= 0. 98* 99* ILO (input) INTEGER 100* IHI (input) INTEGER 101* It is assumed that A is already upper triangular in rows and 102* columns 1:ILO-1 and IHI+1:N. 103* 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. 104* 105* A (input/output) REAL array, dimension (LDA, N) 106* On entry, the N-by-N upper Hessenberg matrix A. Elements 107* below the subdiagonal must be zero. 108* If JOB='S', then on exit A and B will have been 109* simultaneously reduced to generalized Schur form. 110* If JOB='E', then on exit A will have been destroyed. 111* The diagonal blocks will be correct, but the off-diagonal 112* portion will be meaningless. 113* 114* LDA (input) INTEGER 115* The leading dimension of the array A. LDA >= max( 1, N ). 116* 117* B (input/output) REAL array, dimension (LDB, N) 118* On entry, the N-by-N upper triangular matrix B. Elements 119* below the diagonal must be zero. 2-by-2 blocks in B 120* corresponding to 2-by-2 blocks in A will be reduced to 121* positive diagonal form. (I.e., if A(j+1,j) is non-zero, 122* then B(j+1,j)=B(j,j+1)=0 and B(j,j) and B(j+1,j+1) will be 123* positive.) 124* If JOB='S', then on exit A and B will have been 125* simultaneously reduced to Schur form. 126* If JOB='E', then on exit B will have been destroyed. 127* Elements corresponding to diagonal blocks of A will be 128* correct, but the off-diagonal portion will be meaningless. 129* 130* LDB (input) INTEGER 131* The leading dimension of the array B. LDB >= max( 1, N ). 132* 133* ALPHAR (output) REAL array, dimension (N) 134* ALPHAR(1:N) will be set to real parts of the diagonal 135* elements of A that would result from reducing A and B to 136* Schur form and then further reducing them both to triangular 137* form using unitary transformations s.t. the diagonal of B 138* was non-negative real. Thus, if A(j,j) is in a 1-by-1 block 139* (i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=A(j,j). 140* Note that the (real or complex) values 141* (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the 142* generalized eigenvalues of the matrix pencil A - wB. 143* 144* ALPHAI (output) REAL array, dimension (N) 145* ALPHAI(1:N) will be set to imaginary parts of the diagonal 146* elements of A that would result from reducing A and B to 147* Schur form and then further reducing them both to triangular 148* form using unitary transformations s.t. the diagonal of B 149* was non-negative real. Thus, if A(j,j) is in a 1-by-1 block 150* (i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=0. 151* Note that the (real or complex) values 152* (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the 153* generalized eigenvalues of the matrix pencil A - wB. 154* 155* BETA (output) REAL array, dimension (N) 156* BETA(1:N) will be set to the (real) diagonal elements of B 157* that would result from reducing A and B to Schur form and 158* then further reducing them both to triangular form using 159* unitary transformations s.t. the diagonal of B was 160* non-negative real. Thus, if A(j,j) is in a 1-by-1 block 161* (i.e., A(j+1,j)=A(j,j+1)=0), then BETA(j)=B(j,j). 162* Note that the (real or complex) values 163* (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the 164* generalized eigenvalues of the matrix pencil A - wB. 165* (Note that BETA(1:N) will always be non-negative, and no 166* BETAI is necessary.) 167* 168* Q (input/output) REAL array, dimension (LDQ, N) 169* If COMPQ='N', then Q will not be referenced. 170* If COMPQ='V' or 'I', then the transpose of the orthogonal 171* transformations which are applied to A and B on the left 172* will be applied to the array Q on the right. 173* 174* LDQ (input) INTEGER 175* The leading dimension of the array Q. LDQ >= 1. 176* If COMPQ='V' or 'I', then LDQ >= N. 177* 178* Z (input/output) REAL array, dimension (LDZ, N) 179* If COMPZ='N', then Z will not be referenced. 180* If COMPZ='V' or 'I', then the orthogonal transformations 181* which are applied to A and B on the right will be applied 182* to the array Z on the right. 183* 184* LDZ (input) INTEGER 185* The leading dimension of the array Z. LDZ >= 1. 186* If COMPZ='V' or 'I', then LDZ >= N. 187* 188* WORK (workspace/output) REAL array, dimension (LWORK) 189* On exit, if INFO >= 0, WORK(1) returns the optimal LWORK. 190* 191* LWORK (input) INTEGER 192* The dimension of the array WORK. LWORK >= max(1,N). 193* 194* If LWORK = -1, then a workspace query is assumed; the routine 195* only calculates the optimal size of the WORK array, returns 196* this value as the first entry of the WORK array, and no error 197* message related to LWORK is issued by XERBLA. 198* 199* INFO (output) INTEGER 200* = 0: successful exit 201* < 0: if INFO = -i, the i-th argument had an illegal value 202* = 1,...,N: the QZ iteration did not converge. (A,B) is not 203* in Schur form, but ALPHAR(i), ALPHAI(i), and 204* BETA(i), i=INFO+1,...,N should be correct. 205* = N+1,...,2*N: the shift calculation failed. (A,B) is not 206* in Schur form, but ALPHAR(i), ALPHAI(i), and 207* BETA(i), i=INFO-N+1,...,N should be correct. 208* > 2*N: various "impossible" errors. 209* 210* Further Details 211* =============== 212* 213* Iteration counters: 214* 215* JITER -- counts iterations. 216* IITER -- counts iterations run since ILAST was last 217* changed. This is therefore reset only when a 1-by-1 or 218* 2-by-2 block deflates off the bottom. 219* 220* ===================================================================== 221* 222* .. Parameters .. 223 REAL HALF, ZERO, ONE, SAFETY 224 PARAMETER ( HALF = 0.5E+0, ZERO = 0.0E+0, ONE = 1.0E+0, 225 $ SAFETY = 1.0E+2 ) 226* .. 227* .. Local Scalars .. 228 LOGICAL ILAZR2, ILAZRO, ILPIVT, ILQ, ILSCHR, ILZ, 229 $ LQUERY 230 INTEGER ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST, 231 $ ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER, 232 $ JR, MAXIT, NQ, NZ 233 REAL A11, A12, A1I, A1R, A21, A22, A2I, A2R, AD11, 234 $ AD11L, AD12, AD12L, AD21, AD21L, AD22, AD22L, 235 $ AD32L, AN, ANORM, ASCALE, ATOL, B11, B1A, B1I, 236 $ B1R, B22, B2A, B2I, B2R, BN, BNORM, BSCALE, 237 $ BTOL, C, C11I, C11R, C12, C21, C22I, C22R, CL, 238 $ CQ, CR, CZ, ESHIFT, OPST, S, S1, S1INV, S2, 239 $ SAFMAX, SAFMIN, SCALE, SL, SQI, SQR, SR, SZI, 240 $ SZR, T, TAU, TEMP, TEMP2, TEMPI, TEMPR, U1, 241 $ U12, U12L, U2, ULP, VS, W11, W12, W21, W22, 242 $ WABS, WI, WR, WR2 243* .. 244* .. Local Arrays .. 245 REAL V( 3 ) 246* .. 247* .. External Functions .. 248 LOGICAL LSAME 249 REAL SLAMCH, SLANHS, SLAPY2, SLAPY3 250 EXTERNAL LSAME, SLAMCH, SLANHS, SLAPY2, SLAPY3 251* .. 252* .. External Subroutines .. 253 EXTERNAL SLAG2, SLARFG, SLARTG, SLASET, SLASV2, SROT, 254 $ XERBLA 255* .. 256* .. Intrinsic Functions .. 257 INTRINSIC ABS, MAX, MIN, REAL, SQRT 258* .. 259* .. Executable Statements .. 260* 261* Decode JOB, COMPQ, COMPZ 262* 263 IF( LSAME( JOB, 'E' ) ) THEN 264 ILSCHR = .FALSE. 265 ISCHUR = 1 266 ELSE IF( LSAME( JOB, 'S' ) ) THEN 267 ILSCHR = .TRUE. 268 ISCHUR = 2 269 ELSE 270 ISCHUR = 0 271 END IF 272* 273 IF( LSAME( COMPQ, 'N' ) ) THEN 274 ILQ = .FALSE. 275 ICOMPQ = 1 276 NQ = 0 277 ELSE IF( LSAME( COMPQ, 'V' ) ) THEN 278 ILQ = .TRUE. 279 ICOMPQ = 2 280 NQ = N 281 ELSE IF( LSAME( COMPQ, 'I' ) ) THEN 282 ILQ = .TRUE. 283 ICOMPQ = 3 284 NQ = N 285 ELSE 286 ICOMPQ = 0 287 END IF 288* 289 IF( LSAME( COMPZ, 'N' ) ) THEN 290 ILZ = .FALSE. 291 ICOMPZ = 1 292 NZ = 0 293 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN 294 ILZ = .TRUE. 295 ICOMPZ = 2 296 NZ = N 297 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN 298 ILZ = .TRUE. 299 ICOMPZ = 3 300 NZ = N 301 ELSE 302 ICOMPZ = 0 303 END IF 304* 305* Check Argument Values 306* 307 INFO = 0 308 WORK( 1 ) = MAX( 1, N ) 309 LQUERY = ( LWORK.EQ.-1 ) 310 IF( ISCHUR.EQ.0 ) THEN 311 INFO = -1 312 ELSE IF( ICOMPQ.EQ.0 ) THEN 313 INFO = -2 314 ELSE IF( ICOMPZ.EQ.0 ) THEN 315 INFO = -3 316 ELSE IF( N.LT.0 ) THEN 317 INFO = -4 318 ELSE IF( ILO.LT.1 ) THEN 319 INFO = -5 320 ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN 321 INFO = -6 322 ELSE IF( LDA.LT.N ) THEN 323 INFO = -8 324 ELSE IF( LDB.LT.N ) THEN 325 INFO = -10 326 ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN 327 INFO = -15 328 ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN 329 INFO = -17 330 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN 331 INFO = -19 332 END IF 333 IF( INFO.NE.0 ) THEN 334 CALL XERBLA( 'SHGEQZ', -INFO ) 335 RETURN 336 ELSE IF( LQUERY ) THEN 337 RETURN 338 END IF 339* 340* Quick return if possible 341* 342 IF( N.LE.0 ) THEN 343 WORK( 1 ) = REAL( 1 ) 344* --------------------- Begin Timing Code ----------------------- 345 ITCNT = ZERO 346* ---------------------- End Timing Code ------------------------ 347 RETURN 348 END IF 349* 350* Initialize Q and Z 351* 352 IF( ICOMPQ.EQ.3 ) 353 $ CALL SLASET( 'Full', N, N, ZERO, ONE, Q, LDQ ) 354 IF( ICOMPZ.EQ.3 ) 355 $ CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDZ ) 356* 357* Machine Constants 358* 359 IN = IHI + 1 - ILO 360 SAFMIN = SLAMCH( 'S' ) 361 SAFMAX = ONE / SAFMIN 362 ULP = SLAMCH( 'E' )*SLAMCH( 'B' ) 363 ANORM = SLANHS( 'F', IN, A( ILO, ILO ), LDA, WORK ) 364 BNORM = SLANHS( 'F', IN, B( ILO, ILO ), LDB, WORK ) 365 ATOL = MAX( SAFMIN, ULP*ANORM ) 366 BTOL = MAX( SAFMIN, ULP*BNORM ) 367 ASCALE = ONE / MAX( SAFMIN, ANORM ) 368 BSCALE = ONE / MAX( SAFMIN, BNORM ) 369* 370* Set Eigenvalues IHI+1:N 371* 372 DO 30 J = IHI + 1, N 373 IF( B( J, J ).LT.ZERO ) THEN 374 IF( ILSCHR ) THEN 375 DO 10 JR = 1, J 376 A( JR, J ) = -A( JR, J ) 377 B( JR, J ) = -B( JR, J ) 378 10 CONTINUE 379 ELSE 380 A( J, J ) = -A( J, J ) 381 B( J, J ) = -B( J, J ) 382 END IF 383 IF( ILZ ) THEN 384 DO 20 JR = 1, N 385 Z( JR, J ) = -Z( JR, J ) 386 20 CONTINUE 387 END IF 388 END IF 389 ALPHAR( J ) = A( J, J ) 390 ALPHAI( J ) = ZERO 391 BETA( J ) = B( J, J ) 392 30 CONTINUE 393* 394* ---------------------- Begin Timing Code ------------------------- 395* Count ops for norms, etc. 396 OPST = ZERO 397 OPS = OPS + REAL( 2*N**2+6*N ) 398* ----------------------- End Timing Code -------------------------- 399* 400* 401* If IHI < ILO, skip QZ steps 402* 403 IF( IHI.LT.ILO ) 404 $ GO TO 380 405* 406* MAIN QZ ITERATION LOOP 407* 408* Initialize dynamic indices 409* 410* Eigenvalues ILAST+1:N have been found. 411* Column operations modify rows IFRSTM:whatever. 412* Row operations modify columns whatever:ILASTM. 413* 414* If only eigenvalues are being computed, then 415* IFRSTM is the row of the last splitting row above row ILAST; 416* this is always at least ILO. 417* IITER counts iterations since the last eigenvalue was found, 418* to tell when to use an extraordinary shift. 419* MAXIT is the maximum number of QZ sweeps allowed. 420* 421 ILAST = IHI 422 IF( ILSCHR ) THEN 423 IFRSTM = 1 424 ILASTM = N 425 ELSE 426 IFRSTM = ILO 427 ILASTM = IHI 428 END IF 429 IITER = 0 430 ESHIFT = ZERO 431 MAXIT = 30*( IHI-ILO+1 ) 432* 433 DO 360 JITER = 1, MAXIT 434* 435* Split the matrix if possible. 436* 437* Two tests: 438* 1: A(j,j-1)=0 or j=ILO 439* 2: B(j,j)=0 440* 441 IF( ILAST.EQ.ILO ) THEN 442* 443* Special case: j=ILAST 444* 445 GO TO 80 446 ELSE 447 IF( ABS( A( ILAST, ILAST-1 ) ).LE.ATOL ) THEN 448 A( ILAST, ILAST-1 ) = ZERO 449 GO TO 80 450 END IF 451 END IF 452* 453 IF( ABS( B( ILAST, ILAST ) ).LE.BTOL ) THEN 454 B( ILAST, ILAST ) = ZERO 455 GO TO 70 456 END IF 457* 458* General case: j<ILAST 459* 460 DO 60 J = ILAST - 1, ILO, -1 461* 462* Test 1: for A(j,j-1)=0 or j=ILO 463* 464 IF( J.EQ.ILO ) THEN 465 ILAZRO = .TRUE. 466 ELSE 467 IF( ABS( A( J, J-1 ) ).LE.ATOL ) THEN 468 A( J, J-1 ) = ZERO 469 ILAZRO = .TRUE. 470 ELSE 471 ILAZRO = .FALSE. 472 END IF 473 END IF 474* 475* Test 2: for B(j,j)=0 476* 477 IF( ABS( B( J, J ) ).LT.BTOL ) THEN 478 B( J, J ) = ZERO 479* 480* Test 1a: Check for 2 consecutive small subdiagonals in A 481* 482 ILAZR2 = .FALSE. 483 IF( .NOT.ILAZRO ) THEN 484 TEMP = ABS( A( J, J-1 ) ) 485 TEMP2 = ABS( A( J, J ) ) 486 TEMPR = MAX( TEMP, TEMP2 ) 487 IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN 488 TEMP = TEMP / TEMPR 489 TEMP2 = TEMP2 / TEMPR 490 END IF 491 IF( TEMP*( ASCALE*ABS( A( J+1, J ) ) ).LE.TEMP2* 492 $ ( ASCALE*ATOL ) )ILAZR2 = .TRUE. 493 END IF 494* 495* If both tests pass (1 & 2), i.e., the leading diagonal 496* element of B in the block is zero, split a 1x1 block off 497* at the top. (I.e., at the J-th row/column) The leading 498* diagonal element of the remainder can also be zero, so 499* this may have to be done repeatedly. 500* 501 IF( ILAZRO .OR. ILAZR2 ) THEN 502 DO 40 JCH = J, ILAST - 1 503 TEMP = A( JCH, JCH ) 504 CALL SLARTG( TEMP, A( JCH+1, JCH ), C, S, 505 $ A( JCH, JCH ) ) 506 A( JCH+1, JCH ) = ZERO 507 CALL SROT( ILASTM-JCH, A( JCH, JCH+1 ), LDA, 508 $ A( JCH+1, JCH+1 ), LDA, C, S ) 509 CALL SROT( ILASTM-JCH, B( JCH, JCH+1 ), LDB, 510 $ B( JCH+1, JCH+1 ), LDB, C, S ) 511 IF( ILQ ) 512 $ CALL SROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1, 513 $ C, S ) 514 IF( ILAZR2 ) 515 $ A( JCH, JCH-1 ) = A( JCH, JCH-1 )*C 516 ILAZR2 = .FALSE. 517* 518* --------------- Begin Timing Code ----------------- 519 OPST = OPST + REAL( 7+12*( ILASTM-JCH )+6*NQ ) 520* ---------------- End Timing Code ------------------ 521* 522 IF( ABS( B( JCH+1, JCH+1 ) ).GE.BTOL ) THEN 523 IF( JCH+1.GE.ILAST ) THEN 524 GO TO 80 525 ELSE 526 IFIRST = JCH + 1 527 GO TO 110 528 END IF 529 END IF 530 B( JCH+1, JCH+1 ) = ZERO 531 40 CONTINUE 532 GO TO 70 533 ELSE 534* 535* Only test 2 passed -- chase the zero to B(ILAST,ILAST) 536* Then process as in the case B(ILAST,ILAST)=0 537* 538 DO 50 JCH = J, ILAST - 1 539 TEMP = B( JCH, JCH+1 ) 540 CALL SLARTG( TEMP, B( JCH+1, JCH+1 ), C, S, 541 $ B( JCH, JCH+1 ) ) 542 B( JCH+1, JCH+1 ) = ZERO 543 IF( JCH.LT.ILASTM-1 ) 544 $ CALL SROT( ILASTM-JCH-1, B( JCH, JCH+2 ), LDB, 545 $ B( JCH+1, JCH+2 ), LDB, C, S ) 546 CALL SROT( ILASTM-JCH+2, A( JCH, JCH-1 ), LDA, 547 $ A( JCH+1, JCH-1 ), LDA, C, S ) 548 IF( ILQ ) 549 $ CALL SROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1, 550 $ C, S ) 551 TEMP = A( JCH+1, JCH ) 552 CALL SLARTG( TEMP, A( JCH+1, JCH-1 ), C, S, 553 $ A( JCH+1, JCH ) ) 554 A( JCH+1, JCH-1 ) = ZERO 555 CALL SROT( JCH+1-IFRSTM, A( IFRSTM, JCH ), 1, 556 $ A( IFRSTM, JCH-1 ), 1, C, S ) 557 CALL SROT( JCH-IFRSTM, B( IFRSTM, JCH ), 1, 558 $ B( IFRSTM, JCH-1 ), 1, C, S ) 559 IF( ILZ ) 560 $ CALL SROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1, 561 $ C, S ) 562 50 CONTINUE 563* 564* ---------------- Begin Timing Code ------------------- 565 OPST = OPST + REAL( 26+12*( ILASTM-IFRSTM )+6* 566 $ ( NQ+NZ ) )*REAL( ILAST-J ) 567* ----------------- End Timing Code -------------------- 568* 569 GO TO 70 570 END IF 571 ELSE IF( ILAZRO ) THEN 572* 573* Only test 1 passed -- work on J:ILAST 574* 575 IFIRST = J 576 GO TO 110 577 END IF 578* 579* Neither test passed -- try next J 580* 581 60 CONTINUE 582* 583* (Drop-through is "impossible") 584* 585 INFO = N + 1 586 GO TO 420 587* 588* B(ILAST,ILAST)=0 -- clear A(ILAST,ILAST-1) to split off a 589* 1x1 block. 590* 591 70 CONTINUE 592 TEMP = A( ILAST, ILAST ) 593 CALL SLARTG( TEMP, A( ILAST, ILAST-1 ), C, S, 594 $ A( ILAST, ILAST ) ) 595 A( ILAST, ILAST-1 ) = ZERO 596 CALL SROT( ILAST-IFRSTM, A( IFRSTM, ILAST ), 1, 597 $ A( IFRSTM, ILAST-1 ), 1, C, S ) 598 CALL SROT( ILAST-IFRSTM, B( IFRSTM, ILAST ), 1, 599 $ B( IFRSTM, ILAST-1 ), 1, C, S ) 600 IF( ILZ ) 601 $ CALL SROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S ) 602* 603* --------------------- Begin Timing Code ----------------------- 604 OPST = OPST + REAL( 7+12*( ILAST-IFRSTM )+6*NZ ) 605* ---------------------- End Timing Code ------------------------ 606* 607* 608* A(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI, 609* and BETA 610* 611 80 CONTINUE 612 IF( B( ILAST, ILAST ).LT.ZERO ) THEN 613 IF( ILSCHR ) THEN 614 DO 90 J = IFRSTM, ILAST 615 A( J, ILAST ) = -A( J, ILAST ) 616 B( J, ILAST ) = -B( J, ILAST ) 617 90 CONTINUE 618 ELSE 619 A( ILAST, ILAST ) = -A( ILAST, ILAST ) 620 B( ILAST, ILAST ) = -B( ILAST, ILAST ) 621 END IF 622 IF( ILZ ) THEN 623 DO 100 J = 1, N 624 Z( J, ILAST ) = -Z( J, ILAST ) 625 100 CONTINUE 626 END IF 627 END IF 628 ALPHAR( ILAST ) = A( ILAST, ILAST ) 629 ALPHAI( ILAST ) = ZERO 630 BETA( ILAST ) = B( ILAST, ILAST ) 631* 632* Go to next block -- exit if finished. 633* 634 ILAST = ILAST - 1 635 IF( ILAST.LT.ILO ) 636 $ GO TO 380 637* 638* Reset counters 639* 640 IITER = 0 641 ESHIFT = ZERO 642 IF( .NOT.ILSCHR ) THEN 643 ILASTM = ILAST 644 IF( IFRSTM.GT.ILAST ) 645 $ IFRSTM = ILO 646 END IF 647 GO TO 350 648* 649* QZ step 650* 651* This iteration only involves rows/columns IFIRST:ILAST. We 652* assume IFIRST < ILAST, and that the diagonal of B is non-zero. 653* 654 110 CONTINUE 655 IITER = IITER + 1 656 IF( .NOT.ILSCHR ) THEN 657 IFRSTM = IFIRST 658 END IF 659* 660* Compute single shifts. 661* 662* At this point, IFIRST < ILAST, and the diagonal elements of 663* B(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in 664* magnitude) 665* 666 IF( ( IITER / 10 )*10.EQ.IITER ) THEN 667* 668* Exceptional shift. Chosen for no particularly good reason. 669* (Single shift only.) 670* 671 IF( ( REAL( MAXIT )*SAFMIN )*ABS( A( ILAST-1, ILAST ) ).LT. 672 $ ABS( B( ILAST-1, ILAST-1 ) ) ) THEN 673 ESHIFT = ESHIFT + A( ILAST-1, ILAST ) / 674 $ B( ILAST-1, ILAST-1 ) 675 ELSE 676 ESHIFT = ESHIFT + ONE / ( SAFMIN*REAL( MAXIT ) ) 677 END IF 678 S1 = ONE 679 WR = ESHIFT 680* 681* ------------------- Begin Timing Code ---------------------- 682 OPST = OPST + REAL( 4 ) 683* -------------------- End Timing Code ----------------------- 684* 685 ELSE 686* 687* Shifts based on the generalized eigenvalues of the 688* bottom-right 2x2 block of A and B. The first eigenvalue 689* returned by SLAG2 is the Wilkinson shift (AEP p.512), 690* 691 CALL SLAG2( A( ILAST-1, ILAST-1 ), LDA, 692 $ B( ILAST-1, ILAST-1 ), LDB, SAFMIN*SAFETY, S1, 693 $ S2, WR, WR2, WI ) 694* 695 TEMP = MAX( S1, SAFMIN*MAX( ONE, ABS( WR ), ABS( WI ) ) ) 696* 697* ------------------- Begin Timing Code ---------------------- 698 OPST = OPST + REAL( 57 ) 699* -------------------- End Timing Code ----------------------- 700* 701 IF( WI.NE.ZERO ) 702 $ GO TO 200 703 END IF 704* 705* Fiddle with shift to avoid overflow 706* 707 TEMP = MIN( ASCALE, ONE )*( HALF*SAFMAX ) 708 IF( S1.GT.TEMP ) THEN 709 SCALE = TEMP / S1 710 ELSE 711 SCALE = ONE 712 END IF 713* 714 TEMP = MIN( BSCALE, ONE )*( HALF*SAFMAX ) 715 IF( ABS( WR ).GT.TEMP ) 716 $ SCALE = MIN( SCALE, TEMP / ABS( WR ) ) 717 S1 = SCALE*S1 718 WR = SCALE*WR 719* 720* Now check for two consecutive small subdiagonals. 721* 722 DO 120 J = ILAST - 1, IFIRST + 1, -1 723 ISTART = J 724 TEMP = ABS( S1*A( J, J-1 ) ) 725 TEMP2 = ABS( S1*A( J, J )-WR*B( J, J ) ) 726 TEMPR = MAX( TEMP, TEMP2 ) 727 IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN 728 TEMP = TEMP / TEMPR 729 TEMP2 = TEMP2 / TEMPR 730 END IF 731 IF( ABS( ( ASCALE*A( J+1, J ) )*TEMP ).LE.( ASCALE*ATOL )* 732 $ TEMP2 )GO TO 130 733 120 CONTINUE 734* 735 ISTART = IFIRST 736 130 CONTINUE 737* 738* Do an implicit single-shift QZ sweep. 739* 740* Initial Q 741* 742 TEMP = S1*A( ISTART, ISTART ) - WR*B( ISTART, ISTART ) 743 TEMP2 = S1*A( ISTART+1, ISTART ) 744 CALL SLARTG( TEMP, TEMP2, C, S, TEMPR ) 745* 746* Sweep 747* 748 DO 190 J = ISTART, ILAST - 1 749 IF( J.GT.ISTART ) THEN 750 TEMP = A( J, J-1 ) 751 CALL SLARTG( TEMP, A( J+1, J-1 ), C, S, A( J, J-1 ) ) 752 A( J+1, J-1 ) = ZERO 753 END IF 754* 755 DO 140 JC = J, ILASTM 756 TEMP = C*A( J, JC ) + S*A( J+1, JC ) 757 A( J+1, JC ) = -S*A( J, JC ) + C*A( J+1, JC ) 758 A( J, JC ) = TEMP 759 TEMP2 = C*B( J, JC ) + S*B( J+1, JC ) 760 B( J+1, JC ) = -S*B( J, JC ) + C*B( J+1, JC ) 761 B( J, JC ) = TEMP2 762 140 CONTINUE 763 IF( ILQ ) THEN 764 DO 150 JR = 1, N 765 TEMP = C*Q( JR, J ) + S*Q( JR, J+1 ) 766 Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 ) 767 Q( JR, J ) = TEMP 768 150 CONTINUE 769 END IF 770* 771 TEMP = B( J+1, J+1 ) 772 CALL SLARTG( TEMP, B( J+1, J ), C, S, B( J+1, J+1 ) ) 773 B( J+1, J ) = ZERO 774* 775 DO 160 JR = IFRSTM, MIN( J+2, ILAST ) 776 TEMP = C*A( JR, J+1 ) + S*A( JR, J ) 777 A( JR, J ) = -S*A( JR, J+1 ) + C*A( JR, J ) 778 A( JR, J+1 ) = TEMP 779 160 CONTINUE 780 DO 170 JR = IFRSTM, J 781 TEMP = C*B( JR, J+1 ) + S*B( JR, J ) 782 B( JR, J ) = -S*B( JR, J+1 ) + C*B( JR, J ) 783 B( JR, J+1 ) = TEMP 784 170 CONTINUE 785 IF( ILZ ) THEN 786 DO 180 JR = 1, N 787 TEMP = C*Z( JR, J+1 ) + S*Z( JR, J ) 788 Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J ) 789 Z( JR, J+1 ) = TEMP 790 180 CONTINUE 791 END IF 792 190 CONTINUE 793* 794* --------------------- Begin Timing Code ----------------------- 795 OPST = OPST + REAL( 6+( ILAST-ISTART )* 796 $ ( 8+14+36+12*( ILASTM-IFRSTM )+6*( NQ+NZ ) ) ) 797* ---------------------- End Timing Code ------------------------ 798* 799 GO TO 350 800* 801* Use Francis double-shift 802* 803* Note: the Francis double-shift should work with real shifts, 804* but only if the block is at least 3x3. 805* This code may break if this point is reached with 806* a 2x2 block with real eigenvalues. 807* 808 200 CONTINUE 809 IF( IFIRST+1.EQ.ILAST ) THEN 810* 811* Special case -- 2x2 block with complex eigenvectors 812* 813* Step 1: Standardize, that is, rotate so that 814* 815* ( B11 0 ) 816* B = ( ) with B11 non-negative. 817* ( 0 B22 ) 818* 819 CALL SLASV2( B( ILAST-1, ILAST-1 ), B( ILAST-1, ILAST ), 820 $ B( ILAST, ILAST ), B22, B11, SR, CR, SL, CL ) 821* 822 IF( B11.LT.ZERO ) THEN 823 CR = -CR 824 SR = -SR 825 B11 = -B11 826 B22 = -B22 827 END IF 828* 829 CALL SROT( ILASTM+1-IFIRST, A( ILAST-1, ILAST-1 ), LDA, 830 $ A( ILAST, ILAST-1 ), LDA, CL, SL ) 831 CALL SROT( ILAST+1-IFRSTM, A( IFRSTM, ILAST-1 ), 1, 832 $ A( IFRSTM, ILAST ), 1, CR, SR ) 833* 834 IF( ILAST.LT.ILASTM ) 835 $ CALL SROT( ILASTM-ILAST, B( ILAST-1, ILAST+1 ), LDB, 836 $ B( ILAST, ILAST+1 ), LDA, CL, SL ) 837 IF( IFRSTM.LT.ILAST-1 ) 838 $ CALL SROT( IFIRST-IFRSTM, B( IFRSTM, ILAST-1 ), 1, 839 $ B( IFRSTM, ILAST ), 1, CR, SR ) 840* 841 IF( ILQ ) 842 $ CALL SROT( N, Q( 1, ILAST-1 ), 1, Q( 1, ILAST ), 1, CL, 843 $ SL ) 844 IF( ILZ ) 845 $ CALL SROT( N, Z( 1, ILAST-1 ), 1, Z( 1, ILAST ), 1, CR, 846 $ SR ) 847* 848 B( ILAST-1, ILAST-1 ) = B11 849 B( ILAST-1, ILAST ) = ZERO 850 B( ILAST, ILAST-1 ) = ZERO 851 B( ILAST, ILAST ) = B22 852* 853* If B22 is negative, negate column ILAST 854* 855 IF( B22.LT.ZERO ) THEN 856 DO 210 J = IFRSTM, ILAST 857 A( J, ILAST ) = -A( J, ILAST ) 858 B( J, ILAST ) = -B( J, ILAST ) 859 210 CONTINUE 860* 861 IF( ILZ ) THEN 862 DO 220 J = 1, N 863 Z( J, ILAST ) = -Z( J, ILAST ) 864 220 CONTINUE 865 END IF 866 END IF 867* 868* Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.) 869* 870* Recompute shift 871* 872 CALL SLAG2( A( ILAST-1, ILAST-1 ), LDA, 873 $ B( ILAST-1, ILAST-1 ), LDB, SAFMIN*SAFETY, S1, 874 $ TEMP, WR, TEMP2, WI ) 875* 876* ------------------- Begin Timing Code ---------------------- 877 OPST = OPST + REAL( 103+12*( ILASTM+ILAST-IFIRST-IFRSTM )+6* 878 $ ( NQ+NZ ) ) 879* -------------------- End Timing Code ----------------------- 880* 881* If standardization has perturbed the shift onto real line, 882* do another (real single-shift) QR step. 883* 884 IF( WI.EQ.ZERO ) 885 $ GO TO 350 886 S1INV = ONE / S1 887* 888* Do EISPACK (QZVAL) computation of alpha and beta 889* 890 A11 = A( ILAST-1, ILAST-1 ) 891 A21 = A( ILAST, ILAST-1 ) 892 A12 = A( ILAST-1, ILAST ) 893 A22 = A( ILAST, ILAST ) 894* 895* Compute complex Givens rotation on right 896* (Assume some element of C = (sA - wB) > unfl ) 897* __ 898* (sA - wB) ( CZ -SZ ) 899* ( SZ CZ ) 900* 901 C11R = S1*A11 - WR*B11 902 C11I = -WI*B11 903 C12 = S1*A12 904 C21 = S1*A21 905 C22R = S1*A22 - WR*B22 906 C22I = -WI*B22 907* 908 IF( ABS( C11R )+ABS( C11I )+ABS( C12 ).GT.ABS( C21 )+ 909 $ ABS( C22R )+ABS( C22I ) ) THEN 910 T = SLAPY3( C12, C11R, C11I ) 911 CZ = C12 / T 912 SZR = -C11R / T 913 SZI = -C11I / T 914 ELSE 915 CZ = SLAPY2( C22R, C22I ) 916 IF( CZ.LE.SAFMIN ) THEN 917 CZ = ZERO 918 SZR = ONE 919 SZI = ZERO 920 ELSE 921 TEMPR = C22R / CZ 922 TEMPI = C22I / CZ 923 T = SLAPY2( CZ, C21 ) 924 CZ = CZ / T 925 SZR = -C21*TEMPR / T 926 SZI = C21*TEMPI / T 927 END IF 928 END IF 929* 930* Compute Givens rotation on left 931* 932* ( CQ SQ ) 933* ( __ ) A or B 934* ( -SQ CQ ) 935* 936 AN = ABS( A11 ) + ABS( A12 ) + ABS( A21 ) + ABS( A22 ) 937 BN = ABS( B11 ) + ABS( B22 ) 938 WABS = ABS( WR ) + ABS( WI ) 939 IF( S1*AN.GT.WABS*BN ) THEN 940 CQ = CZ*B11 941 SQR = SZR*B22 942 SQI = -SZI*B22 943 ELSE 944 A1R = CZ*A11 + SZR*A12 945 A1I = SZI*A12 946 A2R = CZ*A21 + SZR*A22 947 A2I = SZI*A22 948 CQ = SLAPY2( A1R, A1I ) 949 IF( CQ.LE.SAFMIN ) THEN 950 CQ = ZERO 951 SQR = ONE 952 SQI = ZERO 953 ELSE 954 TEMPR = A1R / CQ 955 TEMPI = A1I / CQ 956 SQR = TEMPR*A2R + TEMPI*A2I 957 SQI = TEMPI*A2R - TEMPR*A2I 958 END IF 959 END IF 960 T = SLAPY3( CQ, SQR, SQI ) 961 CQ = CQ / T 962 SQR = SQR / T 963 SQI = SQI / T 964* 965* Compute diagonal elements of QBZ 966* 967 TEMPR = SQR*SZR - SQI*SZI 968 TEMPI = SQR*SZI + SQI*SZR 969 B1R = CQ*CZ*B11 + TEMPR*B22 970 B1I = TEMPI*B22 971 B1A = SLAPY2( B1R, B1I ) 972 B2R = CQ*CZ*B22 + TEMPR*B11 973 B2I = -TEMPI*B11 974 B2A = SLAPY2( B2R, B2I ) 975* 976* Normalize so beta > 0, and Im( alpha1 ) > 0 977* 978 BETA( ILAST-1 ) = B1A 979 BETA( ILAST ) = B2A 980 ALPHAR( ILAST-1 ) = ( WR*B1A )*S1INV 981 ALPHAI( ILAST-1 ) = ( WI*B1A )*S1INV 982 ALPHAR( ILAST ) = ( WR*B2A )*S1INV 983 ALPHAI( ILAST ) = -( WI*B2A )*S1INV 984* 985* ------------------- Begin Timing Code ---------------------- 986 OPST = OPST + REAL( 93 ) 987* -------------------- End Timing Code ----------------------- 988* 989* Step 3: Go to next block -- exit if finished. 990* 991 ILAST = IFIRST - 1 992 IF( ILAST.LT.ILO ) 993 $ GO TO 380 994* 995* Reset counters 996* 997 IITER = 0 998 ESHIFT = ZERO 999 IF( .NOT.ILSCHR ) THEN 1000 ILASTM = ILAST 1001 IF( IFRSTM.GT.ILAST ) 1002 $ IFRSTM = ILO 1003 END IF 1004 GO TO 350 1005 ELSE 1006* 1007* Usual case: 3x3 or larger block, using Francis implicit 1008* double-shift 1009* 1010* 2 1011* Eigenvalue equation is w - c w + d = 0, 1012* 1013* -1 2 -1 1014* so compute 1st column of (A B ) - c A B + d 1015* using the formula in QZIT (from EISPACK) 1016* 1017* We assume that the block is at least 3x3 1018* 1019 AD11 = ( ASCALE*A( ILAST-1, ILAST-1 ) ) / 1020 $ ( BSCALE*B( ILAST-1, ILAST-1 ) ) 1021 AD21 = ( ASCALE*A( ILAST, ILAST-1 ) ) / 1022 $ ( BSCALE*B( ILAST-1, ILAST-1 ) ) 1023 AD12 = ( ASCALE*A( ILAST-1, ILAST ) ) / 1024 $ ( BSCALE*B( ILAST, ILAST ) ) 1025 AD22 = ( ASCALE*A( ILAST, ILAST ) ) / 1026 $ ( BSCALE*B( ILAST, ILAST ) ) 1027 U12 = B( ILAST-1, ILAST ) / B( ILAST, ILAST ) 1028 AD11L = ( ASCALE*A( IFIRST, IFIRST ) ) / 1029 $ ( BSCALE*B( IFIRST, IFIRST ) ) 1030 AD21L = ( ASCALE*A( IFIRST+1, IFIRST ) ) / 1031 $ ( BSCALE*B( IFIRST, IFIRST ) ) 1032 AD12L = ( ASCALE*A( IFIRST, IFIRST+1 ) ) / 1033 $ ( BSCALE*B( IFIRST+1, IFIRST+1 ) ) 1034 AD22L = ( ASCALE*A( IFIRST+1, IFIRST+1 ) ) / 1035 $ ( BSCALE*B( IFIRST+1, IFIRST+1 ) ) 1036 AD32L = ( ASCALE*A( IFIRST+2, IFIRST+1 ) ) / 1037 $ ( BSCALE*B( IFIRST+1, IFIRST+1 ) ) 1038 U12L = B( IFIRST, IFIRST+1 ) / B( IFIRST+1, IFIRST+1 ) 1039* 1040 V( 1 ) = ( AD11-AD11L )*( AD22-AD11L ) - AD12*AD21 + 1041 $ AD21*U12*AD11L + ( AD12L-AD11L*U12L )*AD21L 1042 V( 2 ) = ( ( AD22L-AD11L )-AD21L*U12L-( AD11-AD11L )- 1043 $ ( AD22-AD11L )+AD21*U12 )*AD21L 1044 V( 3 ) = AD32L*AD21L 1045* 1046 ISTART = IFIRST 1047* 1048 CALL SLARFG( 3, V( 1 ), V( 2 ), 1, TAU ) 1049 V( 1 ) = ONE 1050* 1051* Sweep 1052* 1053 DO 290 J = ISTART, ILAST - 2 1054* 1055* All but last elements: use 3x3 Householder transforms. 1056* 1057* Zero (j-1)st column of A 1058* 1059 IF( J.GT.ISTART ) THEN 1060 V( 1 ) = A( J, J-1 ) 1061 V( 2 ) = A( J+1, J-1 ) 1062 V( 3 ) = A( J+2, J-1 ) 1063* 1064 CALL SLARFG( 3, A( J, J-1 ), V( 2 ), 1, TAU ) 1065 V( 1 ) = ONE 1066 A( J+1, J-1 ) = ZERO 1067 A( J+2, J-1 ) = ZERO 1068 END IF 1069* 1070 DO 230 JC = J, ILASTM 1071 TEMP = TAU*( A( J, JC )+V( 2 )*A( J+1, JC )+V( 3 )* 1072 $ A( J+2, JC ) ) 1073 A( J, JC ) = A( J, JC ) - TEMP 1074 A( J+1, JC ) = A( J+1, JC ) - TEMP*V( 2 ) 1075 A( J+2, JC ) = A( J+2, JC ) - TEMP*V( 3 ) 1076 TEMP2 = TAU*( B( J, JC )+V( 2 )*B( J+1, JC )+V( 3 )* 1077 $ B( J+2, JC ) ) 1078 B( J, JC ) = B( J, JC ) - TEMP2 1079 B( J+1, JC ) = B( J+1, JC ) - TEMP2*V( 2 ) 1080 B( J+2, JC ) = B( J+2, JC ) - TEMP2*V( 3 ) 1081 230 CONTINUE 1082 IF( ILQ ) THEN 1083 DO 240 JR = 1, N 1084 TEMP = TAU*( Q( JR, J )+V( 2 )*Q( JR, J+1 )+V( 3 )* 1085 $ Q( JR, J+2 ) ) 1086 Q( JR, J ) = Q( JR, J ) - TEMP 1087 Q( JR, J+1 ) = Q( JR, J+1 ) - TEMP*V( 2 ) 1088 Q( JR, J+2 ) = Q( JR, J+2 ) - TEMP*V( 3 ) 1089 240 CONTINUE 1090 END IF 1091* 1092* Zero j-th column of B (see SLAGBC for details) 1093* 1094* Swap rows to pivot 1095* 1096 ILPIVT = .FALSE. 1097 TEMP = MAX( ABS( B( J+1, J+1 ) ), ABS( B( J+1, J+2 ) ) ) 1098 TEMP2 = MAX( ABS( B( J+2, J+1 ) ), ABS( B( J+2, J+2 ) ) ) 1099 IF( MAX( TEMP, TEMP2 ).LT.SAFMIN ) THEN 1100 SCALE = ZERO 1101 U1 = ONE 1102 U2 = ZERO 1103 GO TO 250 1104 ELSE IF( TEMP.GE.TEMP2 ) THEN 1105 W11 = B( J+1, J+1 ) 1106 W21 = B( J+2, J+1 ) 1107 W12 = B( J+1, J+2 ) 1108 W22 = B( J+2, J+2 ) 1109 U1 = B( J+1, J ) 1110 U2 = B( J+2, J ) 1111 ELSE 1112 W21 = B( J+1, J+1 ) 1113 W11 = B( J+2, J+1 ) 1114 W22 = B( J+1, J+2 ) 1115 W12 = B( J+2, J+2 ) 1116 U2 = B( J+1, J ) 1117 U1 = B( J+2, J ) 1118 END IF 1119* 1120* Swap columns if nec. 1121* 1122 IF( ABS( W12 ).GT.ABS( W11 ) ) THEN 1123 ILPIVT = .TRUE. 1124 TEMP = W12 1125 TEMP2 = W22 1126 W12 = W11 1127 W22 = W21 1128 W11 = TEMP 1129 W21 = TEMP2 1130 END IF 1131* 1132* LU-factor 1133* 1134 TEMP = W21 / W11 1135 U2 = U2 - TEMP*U1 1136 W22 = W22 - TEMP*W12 1137 W21 = ZERO 1138* 1139* Compute SCALE 1140* 1141 SCALE = ONE 1142 IF( ABS( W22 ).LT.SAFMIN ) THEN 1143 SCALE = ZERO 1144 U2 = ONE 1145 U1 = -W12 / W11 1146 GO TO 250 1147 END IF 1148 IF( ABS( W22 ).LT.ABS( U2 ) ) 1149 $ SCALE = ABS( W22 / U2 ) 1150 IF( ABS( W11 ).LT.ABS( U1 ) ) 1151 $ SCALE = MIN( SCALE, ABS( W11 / U1 ) ) 1152* 1153* Solve 1154* 1155 U2 = ( SCALE*U2 ) / W22 1156 U1 = ( SCALE*U1-W12*U2 ) / W11 1157* 1158 250 CONTINUE 1159 IF( ILPIVT ) THEN 1160 TEMP = U2 1161 U2 = U1 1162 U1 = TEMP 1163 END IF 1164* 1165* Compute Householder Vector 1166* 1167 T = SQRT( SCALE**2+U1**2+U2**2 ) 1168 TAU = ONE + SCALE / T 1169 VS = -ONE / ( SCALE+T ) 1170 V( 1 ) = ONE 1171 V( 2 ) = VS*U1 1172 V( 3 ) = VS*U2 1173* 1174* Apply transformations from the right. 1175* 1176 DO 260 JR = IFRSTM, MIN( J+3, ILAST ) 1177 TEMP = TAU*( A( JR, J )+V( 2 )*A( JR, J+1 )+V( 3 )* 1178 $ A( JR, J+2 ) ) 1179 A( JR, J ) = A( JR, J ) - TEMP 1180 A( JR, J+1 ) = A( JR, J+1 ) - TEMP*V( 2 ) 1181 A( JR, J+2 ) = A( JR, J+2 ) - TEMP*V( 3 ) 1182 260 CONTINUE 1183 DO 270 JR = IFRSTM, J + 2 1184 TEMP = TAU*( B( JR, J )+V( 2 )*B( JR, J+1 )+V( 3 )* 1185 $ B( JR, J+2 ) ) 1186 B( JR, J ) = B( JR, J ) - TEMP 1187 B( JR, J+1 ) = B( JR, J+1 ) - TEMP*V( 2 ) 1188 B( JR, J+2 ) = B( JR, J+2 ) - TEMP*V( 3 ) 1189 270 CONTINUE 1190 IF( ILZ ) THEN 1191 DO 280 JR = 1, N 1192 TEMP = TAU*( Z( JR, J )+V( 2 )*Z( JR, J+1 )+V( 3 )* 1193 $ Z( JR, J+2 ) ) 1194 Z( JR, J ) = Z( JR, J ) - TEMP 1195 Z( JR, J+1 ) = Z( JR, J+1 ) - TEMP*V( 2 ) 1196 Z( JR, J+2 ) = Z( JR, J+2 ) - TEMP*V( 3 ) 1197 280 CONTINUE 1198 END IF 1199 B( J+1, J ) = ZERO 1200 B( J+2, J ) = ZERO 1201 290 CONTINUE 1202* 1203* Last elements: Use Givens rotations 1204* 1205* Rotations from the left 1206* 1207 J = ILAST - 1 1208 TEMP = A( J, J-1 ) 1209 CALL SLARTG( TEMP, A( J+1, J-1 ), C, S, A( J, J-1 ) ) 1210 A( J+1, J-1 ) = ZERO 1211* 1212 DO 300 JC = J, ILASTM 1213 TEMP = C*A( J, JC ) + S*A( J+1, JC ) 1214 A( J+1, JC ) = -S*A( J, JC ) + C*A( J+1, JC ) 1215 A( J, JC ) = TEMP 1216 TEMP2 = C*B( J, JC ) + S*B( J+1, JC ) 1217 B( J+1, JC ) = -S*B( J, JC ) + C*B( J+1, JC ) 1218 B( J, JC ) = TEMP2 1219 300 CONTINUE 1220 IF( ILQ ) THEN 1221 DO 310 JR = 1, N 1222 TEMP = C*Q( JR, J ) + S*Q( JR, J+1 ) 1223 Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 ) 1224 Q( JR, J ) = TEMP 1225 310 CONTINUE 1226 END IF 1227* 1228* Rotations from the right. 1229* 1230 TEMP = B( J+1, J+1 ) 1231 CALL SLARTG( TEMP, B( J+1, J ), C, S, B( J+1, J+1 ) ) 1232 B( J+1, J ) = ZERO 1233* 1234 DO 320 JR = IFRSTM, ILAST 1235 TEMP = C*A( JR, J+1 ) + S*A( JR, J ) 1236 A( JR, J ) = -S*A( JR, J+1 ) + C*A( JR, J ) 1237 A( JR, J+1 ) = TEMP 1238 320 CONTINUE 1239 DO 330 JR = IFRSTM, ILAST - 1 1240 TEMP = C*B( JR, J+1 ) + S*B( JR, J ) 1241 B( JR, J ) = -S*B( JR, J+1 ) + C*B( JR, J ) 1242 B( JR, J+1 ) = TEMP 1243 330 CONTINUE 1244 IF( ILZ ) THEN 1245 DO 340 JR = 1, N 1246 TEMP = C*Z( JR, J+1 ) + S*Z( JR, J ) 1247 Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J ) 1248 Z( JR, J+1 ) = TEMP 1249 340 CONTINUE 1250 END IF 1251* 1252* ------------------- Begin Timing Code ---------------------- 1253 OPST = OPST + ( REAL( 14+30-10+52+12*( ILASTM-IFRSTM )+6* 1254 $ ( NQ+NZ ) )+REAL( ILAST-1-ISTART )* 1255 $ REAL( 14+24+90+20*( ILASTM-IFRSTM )+10*( NQ+NZ ) ) ) 1256* -------------------- End Timing Code ----------------------- 1257* 1258* End of Double-Shift code 1259* 1260 END IF 1261* 1262 GO TO 350 1263* 1264* End of iteration loop 1265* 1266 350 CONTINUE 1267* --------------------- Begin Timing Code ----------------------- 1268 OPS = OPS + OPST 1269 OPST = ZERO 1270* ---------------------- End Timing Code ------------------------ 1271* 1272* 1273 360 CONTINUE 1274* 1275* Drop-through = non-convergence 1276* 1277 370 CONTINUE 1278* ---------------------- Begin Timing Code ------------------------- 1279 OPS = OPS + OPST 1280 OPST = ZERO 1281* ----------------------- End Timing Code -------------------------- 1282* 1283 INFO = ILAST 1284 GO TO 420 1285* 1286* Successful completion of all QZ steps 1287* 1288 380 CONTINUE 1289* 1290* Set Eigenvalues 1:ILO-1 1291* 1292 DO 410 J = 1, ILO - 1 1293 IF( B( J, J ).LT.ZERO ) THEN 1294 IF( ILSCHR ) THEN 1295 DO 390 JR = 1, J 1296 A( JR, J ) = -A( JR, J ) 1297 B( JR, J ) = -B( JR, J ) 1298 390 CONTINUE 1299 ELSE 1300 A( J, J ) = -A( J, J ) 1301 B( J, J ) = -B( J, J ) 1302 END IF 1303 IF( ILZ ) THEN 1304 DO 400 JR = 1, N 1305 Z( JR, J ) = -Z( JR, J ) 1306 400 CONTINUE 1307 END IF 1308 END IF 1309 ALPHAR( J ) = A( J, J ) 1310 ALPHAI( J ) = ZERO 1311 BETA( J ) = B( J, J ) 1312 410 CONTINUE 1313* 1314* Normal Termination 1315* 1316 INFO = 0 1317* 1318* Exit (other than argument error) -- return optimal workspace size 1319* 1320 420 CONTINUE 1321* 1322* ---------------------- Begin Timing Code ------------------------- 1323 OPS = OPS + OPST 1324 OPST = ZERO 1325 ITCNT = JITER 1326* ----------------------- End Timing Code -------------------------- 1327* 1328 WORK( 1 ) = REAL( N ) 1329 RETURN 1330* 1331* End of SHGEQZ 1332* 1333 END 1334