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