1 SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, 2 $ ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, 3 $ RWORK, 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 REAL RWORK( * ) 16 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ), 17 $ BETA( * ), Q( LDQ, * ), WORK( * ), Z( LDZ, * ) 18* .. 19* 20* Purpose 21* ======= 22* 23* CHGEQZ implements a single-shift version of the QZ 24* method for finding the generalized eigenvalues w(i)=ALPHA(i)/BETA(i) 25* of the equation 26* 27* det( A - w(i) B ) = 0 28* 29* If JOB='S', then the pair (A,B) is simultaneously 30* reduced to Schur form (i.e., A and B are both upper triangular) by 31* applying one unitary tranformation (usually called Q) on the left and 32* another (usually called Z) on the right. The diagonal elements of 33* A are then ALPHA(1),...,ALPHA(N), and of B are BETA(1),...,BETA(N). 34* 35* If JOB='S' and COMPQ and COMPZ are 'V' or 'I', then the unitary 36* transformations used to reduce (A,B) are accumulated into the arrays 37* Q and Z s.t.: 38* 39* Q(in) A(in) Z(in)* = Q(out) A(out) Z(out)* 40* Q(in) B(in) Z(in)* = Q(out) B(out) Z(out)* 41* 42* Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix 43* Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973), 44* pp. 241--256. 45* 46* Arguments 47* ========= 48* 49* JOB (input) CHARACTER*1 50* = 'E': compute only ALPHA and BETA. A and B will not 51* necessarily be put into generalized Schur form. 52* = 'S': put A and B into generalized Schur form, as well 53* as computing ALPHA and BETA. 54* 55* COMPQ (input) CHARACTER*1 56* = 'N': do not modify Q. 57* = 'V': multiply the array Q on the right by the conjugate 58* transpose of the unitary tranformation that is 59* applied to the left side of A and B to reduce them 60* to Schur form. 61* = 'I': like COMPQ='V', except that Q will be initialized to 62* the identity first. 63* 64* COMPZ (input) CHARACTER*1 65* = 'N': do not modify Z. 66* = 'V': multiply the array Z on the right by the unitary 67* tranformation that is applied to the right side of 68* A and B to reduce them to Schur form. 69* = 'I': like COMPZ='V', except that Z will be initialized to 70* the identity first. 71* 72* N (input) INTEGER 73* The order of the matrices A, B, Q, and Z. N >= 0. 74* 75* ILO (input) INTEGER 76* IHI (input) INTEGER 77* It is assumed that A is already upper triangular in rows and 78* columns 1:ILO-1 and IHI+1:N. 79* 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. 80* 81* A (input/output) COMPLEX array, dimension (LDA, N) 82* On entry, the N-by-N upper Hessenberg matrix A. Elements 83* below the subdiagonal must be zero. 84* If JOB='S', then on exit A and B will have been 85* simultaneously reduced to upper triangular form. 86* If JOB='E', then on exit A will have been destroyed. 87* 88* LDA (input) INTEGER 89* The leading dimension of the array A. LDA >= max( 1, N ). 90* 91* B (input/output) COMPLEX array, dimension (LDB, N) 92* On entry, the N-by-N upper triangular matrix B. Elements 93* below the diagonal must be zero. 94* If JOB='S', then on exit A and B will have been 95* simultaneously reduced to upper triangular form. 96* If JOB='E', then on exit B will have been destroyed. 97* 98* LDB (input) INTEGER 99* The leading dimension of the array B. LDB >= max( 1, N ). 100* 101* ALPHA (output) COMPLEX array, dimension (N) 102* The diagonal elements of A when the pair (A,B) has been 103* reduced to Schur form. ALPHA(i)/BETA(i) i=1,...,N 104* are the generalized eigenvalues. 105* 106* BETA (output) COMPLEX array, dimension (N) 107* The diagonal elements of B when the pair (A,B) has been 108* reduced to Schur form. ALPHA(i)/BETA(i) i=1,...,N 109* are the generalized eigenvalues. A and B are normalized 110* so that BETA(1),...,BETA(N) are non-negative real numbers. 111* 112* Q (input/output) COMPLEX array, dimension (LDQ, N) 113* If COMPQ='N', then Q will not be referenced. 114* If COMPQ='V' or 'I', then the conjugate transpose of the 115* unitary transformations which are applied to A and B on 116* the left will be applied to the array Q on the right. 117* 118* LDQ (input) INTEGER 119* The leading dimension of the array Q. LDQ >= 1. 120* If COMPQ='V' or 'I', then LDQ >= N. 121* 122* Z (input/output) COMPLEX array, dimension (LDZ, N) 123* If COMPZ='N', then Z will not be referenced. 124* If COMPZ='V' or 'I', then the unitary transformations which 125* are applied to A and B on the right will be applied to the 126* array Z on the right. 127* 128* LDZ (input) INTEGER 129* The leading dimension of the array Z. LDZ >= 1. 130* If COMPZ='V' or 'I', then LDZ >= N. 131* 132* WORK (workspace/output) COMPLEX array, dimension (LWORK) 133* On exit, if INFO >= 0, WORK(1) returns the optimal LWORK. 134* 135* LWORK (input) INTEGER 136* The dimension of the array WORK. LWORK >= max(1,N). 137* 138* If LWORK = -1, then a workspace query is assumed; the routine 139* only calculates the optimal size of the WORK array, returns 140* this value as the first entry of the WORK array, and no error 141* message related to LWORK is issued by XERBLA. 142* 143* RWORK (workspace) REAL array, dimension (N) 144* 145* INFO (output) INTEGER 146* = 0: successful exit 147* < 0: if INFO = -i, the i-th argument had an illegal value 148* = 1,...,N: the QZ iteration did not converge. (A,B) is not 149* in Schur form, but ALPHA(i) and BETA(i), 150* i=INFO+1,...,N should be correct. 151* = N+1,...,2*N: the shift calculation failed. (A,B) is not 152* in Schur form, but ALPHA(i) and BETA(i), 153* i=INFO-N+1,...,N should be correct. 154* > 2*N: various "impossible" errors. 155* 156* Further Details 157* =============== 158* 159* We assume that complex ABS works as long as its value is less than 160* overflow. 161* 162* ===================================================================== 163* 164* .. Parameters .. 165 COMPLEX CZERO, CONE 166 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 167 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 168 REAL ZERO, ONE 169 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 170 REAL HALF 171 PARAMETER ( HALF = 0.5E+0 ) 172* .. 173* .. Local Scalars .. 174 LOGICAL ILAZR2, ILAZRO, ILQ, ILSCHR, ILZ, LQUERY 175 INTEGER ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST, 176 $ ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER, 177 $ JR, MAXIT 178 REAL ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL, 179 $ C, SAFMIN, TEMP, TEMP2, TEMPR, ULP 180 COMPLEX ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2, 181 $ CTEMP3, ESHIFT, RTDISC, S, SHIFT, SIGNBC, T, 182 $ U12, X 183* .. 184* .. External Functions .. 185 LOGICAL LSAME 186 REAL CLANHS, SLAMCH 187 EXTERNAL LSAME, CLANHS, SLAMCH 188* .. 189* .. External Subroutines .. 190 EXTERNAL CLARTG, CLASET, CROT, CSCAL, XERBLA 191* .. 192* .. Intrinsic Functions .. 193 INTRINSIC ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL, SQRT 194* .. 195* .. Statement Functions .. 196 REAL ABS1 197* .. 198* .. Statement Function definitions .. 199 ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) ) 200* .. 201* .. Executable Statements .. 202* 203* Decode JOB, COMPQ, COMPZ 204* 205 IF( LSAME( JOB, 'E' ) ) THEN 206 ILSCHR = .FALSE. 207 ISCHUR = 1 208 ELSE IF( LSAME( JOB, 'S' ) ) THEN 209 ILSCHR = .TRUE. 210 ISCHUR = 2 211 ELSE 212 ISCHUR = 0 213 END IF 214* 215 IF( LSAME( COMPQ, 'N' ) ) THEN 216 ILQ = .FALSE. 217 ICOMPQ = 1 218 ELSE IF( LSAME( COMPQ, 'V' ) ) THEN 219 ILQ = .TRUE. 220 ICOMPQ = 2 221 ELSE IF( LSAME( COMPQ, 'I' ) ) THEN 222 ILQ = .TRUE. 223 ICOMPQ = 3 224 ELSE 225 ICOMPQ = 0 226 END IF 227* 228 IF( LSAME( COMPZ, 'N' ) ) THEN 229 ILZ = .FALSE. 230 ICOMPZ = 1 231 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN 232 ILZ = .TRUE. 233 ICOMPZ = 2 234 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN 235 ILZ = .TRUE. 236 ICOMPZ = 3 237 ELSE 238 ICOMPZ = 0 239 END IF 240* 241* Check Argument Values 242* 243 INFO = 0 244 WORK( 1 ) = MAX( 1, N ) 245 LQUERY = ( LWORK.EQ.-1 ) 246 IF( ISCHUR.EQ.0 ) THEN 247 INFO = -1 248 ELSE IF( ICOMPQ.EQ.0 ) THEN 249 INFO = -2 250 ELSE IF( ICOMPZ.EQ.0 ) THEN 251 INFO = -3 252 ELSE IF( N.LT.0 ) THEN 253 INFO = -4 254 ELSE IF( ILO.LT.1 ) THEN 255 INFO = -5 256 ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN 257 INFO = -6 258 ELSE IF( LDA.LT.N ) THEN 259 INFO = -8 260 ELSE IF( LDB.LT.N ) THEN 261 INFO = -10 262 ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN 263 INFO = -14 264 ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN 265 INFO = -16 266 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN 267 INFO = -18 268 END IF 269 IF( INFO.NE.0 ) THEN 270 CALL XERBLA( 'CHGEQZ', -INFO ) 271 RETURN 272 ELSE IF( LQUERY ) THEN 273 RETURN 274 END IF 275* 276* Quick return if possible 277* 278c WORK( 1 ) = CMPLX( 1 ) 279 IF( N.LE.0 ) THEN 280 WORK( 1 ) = CMPLX( 1 ) 281 RETURN 282 END IF 283* 284* Initialize Q and Z 285* 286 IF( ICOMPQ.EQ.3 ) 287 $ CALL CLASET( 'Full', N, N, CZERO, CONE, Q, LDQ ) 288 IF( ICOMPZ.EQ.3 ) 289 $ CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDZ ) 290* 291* Machine Constants 292* 293 IN = IHI + 1 - ILO 294 SAFMIN = SLAMCH( 'S' ) 295 ULP = SLAMCH( 'E' )*SLAMCH( 'B' ) 296 ANORM = CLANHS( 'F', IN, A( ILO, ILO ), LDA, RWORK ) 297 BNORM = CLANHS( 'F', IN, B( ILO, ILO ), LDB, RWORK ) 298 ATOL = MAX( SAFMIN, ULP*ANORM ) 299 BTOL = MAX( SAFMIN, ULP*BNORM ) 300 ASCALE = ONE / MAX( SAFMIN, ANORM ) 301 BSCALE = ONE / MAX( SAFMIN, BNORM ) 302* 303* 304* Set Eigenvalues IHI+1:N 305* 306 DO 10 J = IHI + 1, N 307 ABSB = ABS( B( J, J ) ) 308 IF( ABSB.GT.SAFMIN ) THEN 309 SIGNBC = CONJG( B( J, J ) / ABSB ) 310 B( J, J ) = ABSB 311 IF( ILSCHR ) THEN 312 CALL CSCAL( J-1, SIGNBC, B( 1, J ), 1 ) 313 CALL CSCAL( J, SIGNBC, A( 1, J ), 1 ) 314 ELSE 315 A( J, J ) = A( J, J )*SIGNBC 316 END IF 317 IF( ILZ ) 318 $ CALL CSCAL( N, SIGNBC, Z( 1, J ), 1 ) 319 ELSE 320 B( J, J ) = CZERO 321 END IF 322 ALPHA( J ) = A( J, J ) 323 BETA( J ) = B( J, J ) 324 10 CONTINUE 325* 326* If IHI < ILO, skip QZ steps 327* 328 IF( IHI.LT.ILO ) 329 $ GO TO 190 330* 331* MAIN QZ ITERATION LOOP 332* 333* Initialize dynamic indices 334* 335* Eigenvalues ILAST+1:N have been found. 336* Column operations modify rows IFRSTM:whatever 337* Row operations modify columns whatever:ILASTM 338* 339* If only eigenvalues are being computed, then 340* IFRSTM is the row of the last splitting row above row ILAST; 341* this is always at least ILO. 342* IITER counts iterations since the last eigenvalue was found, 343* to tell when to use an extraordinary shift. 344* MAXIT is the maximum number of QZ sweeps allowed. 345* 346 ILAST = IHI 347 IF( ILSCHR ) THEN 348 IFRSTM = 1 349 ILASTM = N 350 ELSE 351 IFRSTM = ILO 352 ILASTM = IHI 353 END IF 354 IITER = 0 355 ESHIFT = CZERO 356 MAXIT = 30*( IHI-ILO+1 ) 357* 358 DO 170 JITER = 1, MAXIT 359* 360* Check for too many iterations. 361* 362 IF( JITER.GT.MAXIT ) 363 $ GO TO 180 364* 365* Split the matrix if possible. 366* 367* Two tests: 368* 1: A(j,j-1)=0 or j=ILO 369* 2: B(j,j)=0 370* 371* Special case: j=ILAST 372* 373 IF( ILAST.EQ.ILO ) THEN 374 GO TO 60 375 ELSE 376 IF( ABS1( A( ILAST, ILAST-1 ) ).LE.ATOL ) THEN 377 A( ILAST, ILAST-1 ) = CZERO 378 GO TO 60 379 END IF 380 END IF 381* 382 IF( ABS( B( ILAST, ILAST ) ).LE.BTOL ) THEN 383 B( ILAST, ILAST ) = CZERO 384 GO TO 50 385 END IF 386* 387* General case: j<ILAST 388* 389 DO 40 J = ILAST - 1, ILO, -1 390* 391* Test 1: for A(j,j-1)=0 or j=ILO 392* 393 IF( J.EQ.ILO ) THEN 394 ILAZRO = .TRUE. 395 ELSE 396 IF( ABS1( A( J, J-1 ) ).LE.ATOL ) THEN 397 A( J, J-1 ) = CZERO 398 ILAZRO = .TRUE. 399 ELSE 400 ILAZRO = .FALSE. 401 END IF 402 END IF 403* 404* Test 2: for B(j,j)=0 405* 406 IF( ABS( B( J, J ) ).LT.BTOL ) THEN 407 B( J, J ) = CZERO 408* 409* Test 1a: Check for 2 consecutive small subdiagonals in A 410* 411 ILAZR2 = .FALSE. 412 IF( .NOT.ILAZRO ) THEN 413 IF( ABS1( A( J, J-1 ) )*( ASCALE*ABS1( A( J+1, 414 $ J ) ) ).LE.ABS1( A( J, J ) )*( ASCALE*ATOL ) ) 415 $ ILAZR2 = .TRUE. 416 END IF 417* 418* If both tests pass (1 & 2), i.e., the leading diagonal 419* element of B in the block is zero, split a 1x1 block off 420* at the top. (I.e., at the J-th row/column) The leading 421* diagonal element of the remainder can also be zero, so 422* this may have to be done repeatedly. 423* 424 IF( ILAZRO .OR. ILAZR2 ) THEN 425 DO 20 JCH = J, ILAST - 1 426 CTEMP = A( JCH, JCH ) 427 CALL CLARTG( CTEMP, A( JCH+1, JCH ), C, S, 428 $ A( JCH, JCH ) ) 429 A( JCH+1, JCH ) = CZERO 430 CALL CROT( ILASTM-JCH, A( JCH, JCH+1 ), LDA, 431 $ A( JCH+1, JCH+1 ), LDA, C, S ) 432 CALL CROT( ILASTM-JCH, B( JCH, JCH+1 ), LDB, 433 $ B( JCH+1, JCH+1 ), LDB, C, S ) 434 IF( ILQ ) 435 $ CALL CROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1, 436 $ C, CONJG( S ) ) 437 IF( ILAZR2 ) 438 $ A( JCH, JCH-1 ) = A( JCH, JCH-1 )*C 439 ILAZR2 = .FALSE. 440 IF( ABS1( B( JCH+1, JCH+1 ) ).GE.BTOL ) THEN 441 IF( JCH+1.GE.ILAST ) THEN 442 GO TO 60 443 ELSE 444 IFIRST = JCH + 1 445 GO TO 70 446 END IF 447 END IF 448 B( JCH+1, JCH+1 ) = CZERO 449 20 CONTINUE 450 GO TO 50 451 ELSE 452* 453* Only test 2 passed -- chase the zero to B(ILAST,ILAST) 454* Then process as in the case B(ILAST,ILAST)=0 455* 456 DO 30 JCH = J, ILAST - 1 457 CTEMP = B( JCH, JCH+1 ) 458 CALL CLARTG( CTEMP, B( JCH+1, JCH+1 ), C, S, 459 $ B( JCH, JCH+1 ) ) 460 B( JCH+1, JCH+1 ) = CZERO 461 IF( JCH.LT.ILASTM-1 ) 462 $ CALL CROT( ILASTM-JCH-1, B( JCH, JCH+2 ), LDB, 463 $ B( JCH+1, JCH+2 ), LDB, C, S ) 464 CALL CROT( ILASTM-JCH+2, A( JCH, JCH-1 ), LDA, 465 $ A( JCH+1, JCH-1 ), LDA, C, S ) 466 IF( ILQ ) 467 $ CALL CROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1, 468 $ C, CONJG( S ) ) 469 CTEMP = A( JCH+1, JCH ) 470 CALL CLARTG( CTEMP, A( JCH+1, JCH-1 ), C, S, 471 $ A( JCH+1, JCH ) ) 472 A( JCH+1, JCH-1 ) = CZERO 473 CALL CROT( JCH+1-IFRSTM, A( IFRSTM, JCH ), 1, 474 $ A( IFRSTM, JCH-1 ), 1, C, S ) 475 CALL CROT( JCH-IFRSTM, B( IFRSTM, JCH ), 1, 476 $ B( IFRSTM, JCH-1 ), 1, C, S ) 477 IF( ILZ ) 478 $ CALL CROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1, 479 $ C, S ) 480 30 CONTINUE 481 GO TO 50 482 END IF 483 ELSE IF( ILAZRO ) THEN 484* 485* Only test 1 passed -- work on J:ILAST 486* 487 IFIRST = J 488 GO TO 70 489 END IF 490* 491* Neither test passed -- try next J 492* 493 40 CONTINUE 494* 495* (Drop-through is "impossible") 496* 497 INFO = 2*N + 1 498 GO TO 210 499* 500* B(ILAST,ILAST)=0 -- clear A(ILAST,ILAST-1) to split off a 501* 1x1 block. 502* 503 50 CONTINUE 504 CTEMP = A( ILAST, ILAST ) 505 CALL CLARTG( CTEMP, A( ILAST, ILAST-1 ), C, S, 506 $ A( ILAST, ILAST ) ) 507 A( ILAST, ILAST-1 ) = CZERO 508 CALL CROT( ILAST-IFRSTM, A( IFRSTM, ILAST ), 1, 509 $ A( IFRSTM, ILAST-1 ), 1, C, S ) 510 CALL CROT( ILAST-IFRSTM, B( IFRSTM, ILAST ), 1, 511 $ B( IFRSTM, ILAST-1 ), 1, C, S ) 512 IF( ILZ ) 513 $ CALL CROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S ) 514* 515* A(ILAST,ILAST-1)=0 -- Standardize B, set ALPHA and BETA 516* 517 60 CONTINUE 518 ABSB = ABS( B( ILAST, ILAST ) ) 519 IF( ABSB.GT.SAFMIN ) THEN 520 SIGNBC = CONJG( B( ILAST, ILAST ) / ABSB ) 521 B( ILAST, ILAST ) = ABSB 522 IF( ILSCHR ) THEN 523 CALL CSCAL( ILAST-IFRSTM, SIGNBC, B( IFRSTM, ILAST ), 1 ) 524 CALL CSCAL( ILAST+1-IFRSTM, SIGNBC, A( IFRSTM, ILAST ), 525 $ 1 ) 526 ELSE 527 A( ILAST, ILAST ) = A( ILAST, ILAST )*SIGNBC 528 END IF 529 IF( ILZ ) 530 $ CALL CSCAL( N, SIGNBC, Z( 1, ILAST ), 1 ) 531 ELSE 532 B( ILAST, ILAST ) = CZERO 533 END IF 534 ALPHA( ILAST ) = A( ILAST, ILAST ) 535 BETA( ILAST ) = B( ILAST, ILAST ) 536* 537* Go to next block -- exit if finished. 538* 539 ILAST = ILAST - 1 540 IF( ILAST.LT.ILO ) 541 $ GO TO 190 542* 543* Reset counters 544* 545 IITER = 0 546 ESHIFT = CZERO 547 IF( .NOT.ILSCHR ) THEN 548 ILASTM = ILAST 549 IF( IFRSTM.GT.ILAST ) 550 $ IFRSTM = ILO 551 END IF 552 GO TO 160 553* 554* QZ step 555* 556* This iteration only involves rows/columns IFIRST:ILAST. We 557* assume IFIRST < ILAST, and that the diagonal of B is non-zero. 558* 559 70 CONTINUE 560 IITER = IITER + 1 561 IF( .NOT.ILSCHR ) THEN 562 IFRSTM = IFIRST 563 END IF 564* 565* Compute the Shift. 566* 567* At this point, IFIRST < ILAST, and the diagonal elements of 568* B(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in 569* magnitude) 570* 571 IF( ( IITER / 10 )*10.NE.IITER ) THEN 572* 573* The Wilkinson shift (AEP p.512), i.e., the eigenvalue of 574* the bottom-right 2x2 block of A inv(B) which is nearest to 575* the bottom-right element. 576* 577* We factor B as U*D, where U has unit diagonals, and 578* compute (A*inv(D))*inv(U). 579* 580 U12 = ( BSCALE*B( ILAST-1, ILAST ) ) / 581 $ ( BSCALE*B( ILAST, ILAST ) ) 582 AD11 = ( ASCALE*A( ILAST-1, ILAST-1 ) ) / 583 $ ( BSCALE*B( ILAST-1, ILAST-1 ) ) 584 AD21 = ( ASCALE*A( ILAST, ILAST-1 ) ) / 585 $ ( BSCALE*B( ILAST-1, ILAST-1 ) ) 586 AD12 = ( ASCALE*A( ILAST-1, ILAST ) ) / 587 $ ( BSCALE*B( ILAST, ILAST ) ) 588 AD22 = ( ASCALE*A( ILAST, ILAST ) ) / 589 $ ( BSCALE*B( ILAST, ILAST ) ) 590 ABI22 = AD22 - U12*AD21 591* 592 T = HALF*( AD11+ABI22 ) 593 RTDISC = SQRT( T**2+AD12*AD21-AD11*AD22 ) 594 TEMP = REAL( T-ABI22 )*REAL( RTDISC ) + 595 $ AIMAG( T-ABI22 )*AIMAG( RTDISC ) 596 IF( TEMP.LE.ZERO ) THEN 597 SHIFT = T + RTDISC 598 ELSE 599 SHIFT = T - RTDISC 600 END IF 601 ELSE 602* 603* Exceptional shift. Chosen for no particularly good reason. 604* 605 ESHIFT = ESHIFT + CONJG( ( ASCALE*A( ILAST-1, ILAST ) ) / 606 $ ( BSCALE*B( ILAST-1, ILAST-1 ) ) ) 607 SHIFT = ESHIFT 608 END IF 609* 610* Now check for two consecutive small subdiagonals. 611* 612 DO 80 J = ILAST - 1, IFIRST + 1, -1 613 ISTART = J 614 CTEMP = ASCALE*A( J, J ) - SHIFT*( BSCALE*B( J, J ) ) 615 TEMP = ABS1( CTEMP ) 616 TEMP2 = ASCALE*ABS1( A( J+1, J ) ) 617 TEMPR = MAX( TEMP, TEMP2 ) 618 IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN 619 TEMP = TEMP / TEMPR 620 TEMP2 = TEMP2 / TEMPR 621 END IF 622 IF( ABS1( A( J, J-1 ) )*TEMP2.LE.TEMP*ATOL ) 623 $ GO TO 90 624 80 CONTINUE 625* 626 ISTART = IFIRST 627 CTEMP = ASCALE*A( IFIRST, IFIRST ) - 628 $ SHIFT*( BSCALE*B( IFIRST, IFIRST ) ) 629 90 CONTINUE 630* 631* Do an implicit-shift QZ sweep. 632* 633* Initial Q 634* 635 CTEMP2 = ASCALE*A( ISTART+1, ISTART ) 636 CALL CLARTG( CTEMP, CTEMP2, C, S, CTEMP3 ) 637* 638* Sweep 639* 640 DO 150 J = ISTART, ILAST - 1 641 IF( J.GT.ISTART ) THEN 642 CTEMP = A( J, J-1 ) 643 CALL CLARTG( CTEMP, A( J+1, J-1 ), C, S, A( J, J-1 ) ) 644 A( J+1, J-1 ) = CZERO 645 END IF 646* 647 DO 100 JC = J, ILASTM 648 CTEMP = C*A( J, JC ) + S*A( J+1, JC ) 649 A( J+1, JC ) = -CONJG( S )*A( J, JC ) + C*A( J+1, JC ) 650 A( J, JC ) = CTEMP 651 CTEMP2 = C*B( J, JC ) + S*B( J+1, JC ) 652 B( J+1, JC ) = -CONJG( S )*B( J, JC ) + C*B( J+1, JC ) 653 B( J, JC ) = CTEMP2 654 100 CONTINUE 655 IF( ILQ ) THEN 656 DO 110 JR = 1, N 657 CTEMP = C*Q( JR, J ) + CONJG( S )*Q( JR, J+1 ) 658 Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 ) 659 Q( JR, J ) = CTEMP 660 110 CONTINUE 661 END IF 662* 663 CTEMP = B( J+1, J+1 ) 664 CALL CLARTG( CTEMP, B( J+1, J ), C, S, B( J+1, J+1 ) ) 665 B( J+1, J ) = CZERO 666* 667 DO 120 JR = IFRSTM, MIN( J+2, ILAST ) 668 CTEMP = C*A( JR, J+1 ) + S*A( JR, J ) 669 A( JR, J ) = -CONJG( S )*A( JR, J+1 ) + C*A( JR, J ) 670 A( JR, J+1 ) = CTEMP 671 120 CONTINUE 672 DO 130 JR = IFRSTM, J 673 CTEMP = C*B( JR, J+1 ) + S*B( JR, J ) 674 B( JR, J ) = -CONJG( S )*B( JR, J+1 ) + C*B( JR, J ) 675 B( JR, J+1 ) = CTEMP 676 130 CONTINUE 677 IF( ILZ ) THEN 678 DO 140 JR = 1, N 679 CTEMP = C*Z( JR, J+1 ) + S*Z( JR, J ) 680 Z( JR, J ) = -CONJG( S )*Z( JR, J+1 ) + C*Z( JR, J ) 681 Z( JR, J+1 ) = CTEMP 682 140 CONTINUE 683 END IF 684 150 CONTINUE 685* 686 160 CONTINUE 687* 688 170 CONTINUE 689* 690* Drop-through = non-convergence 691* 692 180 CONTINUE 693 INFO = ILAST 694 GO TO 210 695* 696* Successful completion of all QZ steps 697* 698 190 CONTINUE 699* 700* Set Eigenvalues 1:ILO-1 701* 702 DO 200 J = 1, ILO - 1 703 ABSB = ABS( B( J, J ) ) 704 IF( ABSB.GT.SAFMIN ) THEN 705 SIGNBC = CONJG( B( J, J ) / ABSB ) 706 B( J, J ) = ABSB 707 IF( ILSCHR ) THEN 708 CALL CSCAL( J-1, SIGNBC, B( 1, J ), 1 ) 709 CALL CSCAL( J, SIGNBC, A( 1, J ), 1 ) 710 ELSE 711 A( J, J ) = A( J, J )*SIGNBC 712 END IF 713 IF( ILZ ) 714 $ CALL CSCAL( N, SIGNBC, Z( 1, J ), 1 ) 715 ELSE 716 B( J, J ) = CZERO 717 END IF 718 ALPHA( J ) = A( J, J ) 719 BETA( J ) = B( J, J ) 720 200 CONTINUE 721* 722* Normal Termination 723* 724 INFO = 0 725* 726* Exit (other than argument error) -- return optimal workspace size 727* 728 210 CONTINUE 729 WORK( 1 ) = CMPLX( N ) 730 RETURN 731* 732* End of CHGEQZ 733* 734 END 735