1*> \brief \b ZHGEQZ 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZHGEQZ + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhgeqz.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhgeqz.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhgeqz.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, 22* ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, 23* RWORK, INFO ) 24* 25* .. Scalar Arguments .. 26* CHARACTER COMPQ, COMPZ, JOB 27* INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N 28* .. 29* .. Array Arguments .. 30* DOUBLE PRECISION RWORK( * ) 31* COMPLEX*16 ALPHA( * ), BETA( * ), H( LDH, * ), 32* $ Q( LDQ, * ), T( LDT, * ), WORK( * ), 33* $ Z( LDZ, * ) 34* .. 35* 36* 37*> \par Purpose: 38* ============= 39*> 40*> \verbatim 41*> 42*> ZHGEQZ computes the eigenvalues of a complex matrix pair (H,T), 43*> where H is an upper Hessenberg matrix and T is upper triangular, 44*> using the single-shift QZ method. 45*> Matrix pairs of this type are produced by the reduction to 46*> generalized upper Hessenberg form of a complex matrix pair (A,B): 47*> 48*> A = Q1*H*Z1**H, B = Q1*T*Z1**H, 49*> 50*> as computed by ZGGHRD. 51*> 52*> If JOB='S', then the Hessenberg-triangular pair (H,T) is 53*> also reduced to generalized Schur form, 54*> 55*> H = Q*S*Z**H, T = Q*P*Z**H, 56*> 57*> where Q and Z are unitary matrices and S and P are upper triangular. 58*> 59*> Optionally, the unitary matrix Q from the generalized Schur 60*> factorization may be postmultiplied into an input matrix Q1, and the 61*> unitary matrix Z may be postmultiplied into an input matrix Z1. 62*> If Q1 and Z1 are the unitary matrices from ZGGHRD that reduced 63*> the matrix pair (A,B) to generalized Hessenberg form, then the output 64*> matrices Q1*Q and Z1*Z are the unitary factors from the generalized 65*> Schur factorization of (A,B): 66*> 67*> A = (Q1*Q)*S*(Z1*Z)**H, B = (Q1*Q)*P*(Z1*Z)**H. 68*> 69*> To avoid overflow, eigenvalues of the matrix pair (H,T) 70*> (equivalently, of (A,B)) are computed as a pair of complex values 71*> (alpha,beta). If beta is nonzero, lambda = alpha / beta is an 72*> eigenvalue of the generalized nonsymmetric eigenvalue problem (GNEP) 73*> A*x = lambda*B*x 74*> and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the 75*> alternate form of the GNEP 76*> mu*A*y = B*y. 77*> The values of alpha and beta for the i-th eigenvalue can be read 78*> directly from the generalized Schur form: alpha = S(i,i), 79*> beta = P(i,i). 80*> 81*> Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix 82*> Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973), 83*> pp. 241--256. 84*> \endverbatim 85* 86* Arguments: 87* ========== 88* 89*> \param[in] JOB 90*> \verbatim 91*> JOB is CHARACTER*1 92*> = 'E': Compute eigenvalues only; 93*> = 'S': Computer eigenvalues and the Schur form. 94*> \endverbatim 95*> 96*> \param[in] COMPQ 97*> \verbatim 98*> COMPQ is CHARACTER*1 99*> = 'N': Left Schur vectors (Q) are not computed; 100*> = 'I': Q is initialized to the unit matrix and the matrix Q 101*> of left Schur vectors of (H,T) is returned; 102*> = 'V': Q must contain a unitary matrix Q1 on entry and 103*> the product Q1*Q is returned. 104*> \endverbatim 105*> 106*> \param[in] COMPZ 107*> \verbatim 108*> COMPZ is CHARACTER*1 109*> = 'N': Right Schur vectors (Z) are not computed; 110*> = 'I': Q is initialized to the unit matrix and the matrix Z 111*> of right Schur vectors of (H,T) is returned; 112*> = 'V': Z must contain a unitary matrix Z1 on entry and 113*> the product Z1*Z is returned. 114*> \endverbatim 115*> 116*> \param[in] N 117*> \verbatim 118*> N is INTEGER 119*> The order of the matrices H, T, Q, and Z. N >= 0. 120*> \endverbatim 121*> 122*> \param[in] ILO 123*> \verbatim 124*> ILO is INTEGER 125*> \endverbatim 126*> 127*> \param[in] IHI 128*> \verbatim 129*> IHI is INTEGER 130*> ILO and IHI mark the rows and columns of H which are in 131*> Hessenberg form. It is assumed that A is already upper 132*> triangular in rows and columns 1:ILO-1 and IHI+1:N. 133*> If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0. 134*> \endverbatim 135*> 136*> \param[in,out] H 137*> \verbatim 138*> H is COMPLEX*16 array, dimension (LDH, N) 139*> On entry, the N-by-N upper Hessenberg matrix H. 140*> On exit, if JOB = 'S', H contains the upper triangular 141*> matrix S from the generalized Schur factorization. 142*> If JOB = 'E', the diagonal of H matches that of S, but 143*> the rest of H is unspecified. 144*> \endverbatim 145*> 146*> \param[in] LDH 147*> \verbatim 148*> LDH is INTEGER 149*> The leading dimension of the array H. LDH >= max( 1, N ). 150*> \endverbatim 151*> 152*> \param[in,out] T 153*> \verbatim 154*> T is COMPLEX*16 array, dimension (LDT, N) 155*> On entry, the N-by-N upper triangular matrix T. 156*> On exit, if JOB = 'S', T contains the upper triangular 157*> matrix P from the generalized Schur factorization. 158*> If JOB = 'E', the diagonal of T matches that of P, but 159*> the rest of T is unspecified. 160*> \endverbatim 161*> 162*> \param[in] LDT 163*> \verbatim 164*> LDT is INTEGER 165*> The leading dimension of the array T. LDT >= max( 1, N ). 166*> \endverbatim 167*> 168*> \param[out] ALPHA 169*> \verbatim 170*> ALPHA is COMPLEX*16 array, dimension (N) 171*> The complex scalars alpha that define the eigenvalues of 172*> GNEP. ALPHA(i) = S(i,i) in the generalized Schur 173*> factorization. 174*> \endverbatim 175*> 176*> \param[out] BETA 177*> \verbatim 178*> BETA is COMPLEX*16 array, dimension (N) 179*> The real non-negative scalars beta that define the 180*> eigenvalues of GNEP. BETA(i) = P(i,i) in the generalized 181*> Schur factorization. 182*> 183*> Together, the quantities alpha = ALPHA(j) and beta = BETA(j) 184*> represent the j-th eigenvalue of the matrix pair (A,B), in 185*> one of the forms lambda = alpha/beta or mu = beta/alpha. 186*> Since either lambda or mu may overflow, they should not, 187*> in general, be computed. 188*> \endverbatim 189*> 190*> \param[in,out] Q 191*> \verbatim 192*> Q is COMPLEX*16 array, dimension (LDQ, N) 193*> On entry, if COMPQ = 'V', the unitary matrix Q1 used in the 194*> reduction of (A,B) to generalized Hessenberg form. 195*> On exit, if COMPQ = 'I', the unitary matrix of left Schur 196*> vectors of (H,T), and if COMPQ = 'V', the unitary matrix of 197*> left Schur vectors of (A,B). 198*> Not referenced if COMPQ = 'N'. 199*> \endverbatim 200*> 201*> \param[in] LDQ 202*> \verbatim 203*> LDQ is INTEGER 204*> The leading dimension of the array Q. LDQ >= 1. 205*> If COMPQ='V' or 'I', then LDQ >= N. 206*> \endverbatim 207*> 208*> \param[in,out] Z 209*> \verbatim 210*> Z is COMPLEX*16 array, dimension (LDZ, N) 211*> On entry, if COMPZ = 'V', the unitary matrix Z1 used in the 212*> reduction of (A,B) to generalized Hessenberg form. 213*> On exit, if COMPZ = 'I', the unitary matrix of right Schur 214*> vectors of (H,T), and if COMPZ = 'V', the unitary matrix of 215*> right Schur vectors of (A,B). 216*> Not referenced if COMPZ = 'N'. 217*> \endverbatim 218*> 219*> \param[in] LDZ 220*> \verbatim 221*> LDZ is INTEGER 222*> The leading dimension of the array Z. LDZ >= 1. 223*> If COMPZ='V' or 'I', then LDZ >= N. 224*> \endverbatim 225*> 226*> \param[out] WORK 227*> \verbatim 228*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) 229*> On exit, if INFO >= 0, WORK(1) returns the optimal LWORK. 230*> \endverbatim 231*> 232*> \param[in] LWORK 233*> \verbatim 234*> LWORK is INTEGER 235*> The dimension of the array WORK. LWORK >= max(1,N). 236*> 237*> If LWORK = -1, then a workspace query is assumed; the routine 238*> only calculates the optimal size of the WORK array, returns 239*> this value as the first entry of the WORK array, and no error 240*> message related to LWORK is issued by XERBLA. 241*> \endverbatim 242*> 243*> \param[out] RWORK 244*> \verbatim 245*> RWORK is DOUBLE PRECISION array, dimension (N) 246*> \endverbatim 247*> 248*> \param[out] INFO 249*> \verbatim 250*> INFO is INTEGER 251*> = 0: successful exit 252*> < 0: if INFO = -i, the i-th argument had an illegal value 253*> = 1,...,N: the QZ iteration did not converge. (H,T) is not 254*> in Schur form, but ALPHA(i) and BETA(i), 255*> i=INFO+1,...,N should be correct. 256*> = N+1,...,2*N: the shift calculation failed. (H,T) is not 257*> in Schur form, but ALPHA(i) and BETA(i), 258*> i=INFO-N+1,...,N should be correct. 259*> \endverbatim 260* 261* Authors: 262* ======== 263* 264*> \author Univ. of Tennessee 265*> \author Univ. of California Berkeley 266*> \author Univ. of Colorado Denver 267*> \author NAG Ltd. 268* 269*> \ingroup complex16GEcomputational 270* 271*> \par Further Details: 272* ===================== 273*> 274*> \verbatim 275*> 276*> We assume that complex ABS works as long as its value is less than 277*> overflow. 278*> \endverbatim 279*> 280* ===================================================================== 281 SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, 282 $ ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, 283 $ RWORK, INFO ) 284* 285* -- LAPACK computational routine -- 286* -- LAPACK is a software package provided by Univ. of Tennessee, -- 287* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 288* 289* .. Scalar Arguments .. 290 CHARACTER COMPQ, COMPZ, JOB 291 INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N 292* .. 293* .. Array Arguments .. 294 DOUBLE PRECISION RWORK( * ) 295 COMPLEX*16 ALPHA( * ), BETA( * ), H( LDH, * ), 296 $ Q( LDQ, * ), T( LDT, * ), WORK( * ), 297 $ Z( LDZ, * ) 298* .. 299* 300* ===================================================================== 301* 302* .. Parameters .. 303 COMPLEX*16 CZERO, CONE 304 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 305 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 306 DOUBLE PRECISION ZERO, ONE 307 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 308 DOUBLE PRECISION HALF 309 PARAMETER ( HALF = 0.5D+0 ) 310* .. 311* .. Local Scalars .. 312 LOGICAL ILAZR2, ILAZRO, ILQ, ILSCHR, ILZ, LQUERY 313 INTEGER ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST, 314 $ ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER, 315 $ JR, MAXIT 316 DOUBLE PRECISION ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL, 317 $ C, SAFMIN, TEMP, TEMP2, TEMPR, ULP 318 COMPLEX*16 ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2, 319 $ CTEMP3, ESHIFT, S, SHIFT, SIGNBC, 320 $ U12, X, ABI12, Y 321* .. 322* .. External Functions .. 323 COMPLEX*16 ZLADIV 324 LOGICAL LSAME 325 DOUBLE PRECISION DLAMCH, ZLANHS 326 EXTERNAL ZLADIV, LSAME, DLAMCH, ZLANHS 327* .. 328* .. External Subroutines .. 329 EXTERNAL XERBLA, ZLARTG, ZLASET, ZROT, ZSCAL 330* .. 331* .. Intrinsic Functions .. 332 INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN, 333 $ SQRT 334* .. 335* .. Statement Functions .. 336 DOUBLE PRECISION ABS1 337* .. 338* .. Statement Function definitions .. 339 ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) ) 340* .. 341* .. Executable Statements .. 342* 343* Decode JOB, COMPQ, COMPZ 344* 345 IF( LSAME( JOB, 'E' ) ) THEN 346 ILSCHR = .FALSE. 347 ISCHUR = 1 348 ELSE IF( LSAME( JOB, 'S' ) ) THEN 349 ILSCHR = .TRUE. 350 ISCHUR = 2 351 ELSE 352 ILSCHR = .TRUE. 353 ISCHUR = 0 354 END IF 355* 356 IF( LSAME( COMPQ, 'N' ) ) THEN 357 ILQ = .FALSE. 358 ICOMPQ = 1 359 ELSE IF( LSAME( COMPQ, 'V' ) ) THEN 360 ILQ = .TRUE. 361 ICOMPQ = 2 362 ELSE IF( LSAME( COMPQ, 'I' ) ) THEN 363 ILQ = .TRUE. 364 ICOMPQ = 3 365 ELSE 366 ILQ = .TRUE. 367 ICOMPQ = 0 368 END IF 369* 370 IF( LSAME( COMPZ, 'N' ) ) THEN 371 ILZ = .FALSE. 372 ICOMPZ = 1 373 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN 374 ILZ = .TRUE. 375 ICOMPZ = 2 376 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN 377 ILZ = .TRUE. 378 ICOMPZ = 3 379 ELSE 380 ILZ = .TRUE. 381 ICOMPZ = 0 382 END IF 383* 384* Check Argument Values 385* 386 INFO = 0 387 WORK( 1 ) = MAX( 1, N ) 388 LQUERY = ( LWORK.EQ.-1 ) 389 IF( ISCHUR.EQ.0 ) THEN 390 INFO = -1 391 ELSE IF( ICOMPQ.EQ.0 ) THEN 392 INFO = -2 393 ELSE IF( ICOMPZ.EQ.0 ) THEN 394 INFO = -3 395 ELSE IF( N.LT.0 ) THEN 396 INFO = -4 397 ELSE IF( ILO.LT.1 ) THEN 398 INFO = -5 399 ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN 400 INFO = -6 401 ELSE IF( LDH.LT.N ) THEN 402 INFO = -8 403 ELSE IF( LDT.LT.N ) THEN 404 INFO = -10 405 ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN 406 INFO = -14 407 ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN 408 INFO = -16 409 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN 410 INFO = -18 411 END IF 412 IF( INFO.NE.0 ) THEN 413 CALL XERBLA( 'ZHGEQZ', -INFO ) 414 RETURN 415 ELSE IF( LQUERY ) THEN 416 RETURN 417 END IF 418* 419* Quick return if possible 420* 421* WORK( 1 ) = CMPLX( 1 ) 422 IF( N.LE.0 ) THEN 423 WORK( 1 ) = DCMPLX( 1 ) 424 RETURN 425 END IF 426* 427* Initialize Q and Z 428* 429 IF( ICOMPQ.EQ.3 ) 430 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ ) 431 IF( ICOMPZ.EQ.3 ) 432 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ ) 433* 434* Machine Constants 435* 436 IN = IHI + 1 - ILO 437 SAFMIN = DLAMCH( 'S' ) 438 ULP = DLAMCH( 'E' )*DLAMCH( 'B' ) 439 ANORM = ZLANHS( 'F', IN, H( ILO, ILO ), LDH, RWORK ) 440 BNORM = ZLANHS( 'F', IN, T( ILO, ILO ), LDT, RWORK ) 441 ATOL = MAX( SAFMIN, ULP*ANORM ) 442 BTOL = MAX( SAFMIN, ULP*BNORM ) 443 ASCALE = ONE / MAX( SAFMIN, ANORM ) 444 BSCALE = ONE / MAX( SAFMIN, BNORM ) 445* 446* 447* Set Eigenvalues IHI+1:N 448* 449 DO 10 J = IHI + 1, N 450 ABSB = ABS( T( J, J ) ) 451 IF( ABSB.GT.SAFMIN ) THEN 452 SIGNBC = DCONJG( T( J, J ) / ABSB ) 453 T( J, J ) = ABSB 454 IF( ILSCHR ) THEN 455 CALL ZSCAL( J-1, SIGNBC, T( 1, J ), 1 ) 456 CALL ZSCAL( J, SIGNBC, H( 1, J ), 1 ) 457 ELSE 458 CALL ZSCAL( 1, SIGNBC, H( J, J ), 1 ) 459 END IF 460 IF( ILZ ) 461 $ CALL ZSCAL( N, SIGNBC, Z( 1, J ), 1 ) 462 ELSE 463 T( J, J ) = CZERO 464 END IF 465 ALPHA( J ) = H( J, J ) 466 BETA( J ) = T( J, J ) 467 10 CONTINUE 468* 469* If IHI < ILO, skip QZ steps 470* 471 IF( IHI.LT.ILO ) 472 $ GO TO 190 473* 474* MAIN QZ ITERATION LOOP 475* 476* Initialize dynamic indices 477* 478* Eigenvalues ILAST+1:N have been found. 479* Column operations modify rows IFRSTM:whatever 480* Row operations modify columns whatever:ILASTM 481* 482* If only eigenvalues are being computed, then 483* IFRSTM is the row of the last splitting row above row ILAST; 484* this is always at least ILO. 485* IITER counts iterations since the last eigenvalue was found, 486* to tell when to use an extraordinary shift. 487* MAXIT is the maximum number of QZ sweeps allowed. 488* 489 ILAST = IHI 490 IF( ILSCHR ) THEN 491 IFRSTM = 1 492 ILASTM = N 493 ELSE 494 IFRSTM = ILO 495 ILASTM = IHI 496 END IF 497 IITER = 0 498 ESHIFT = CZERO 499 MAXIT = 30*( IHI-ILO+1 ) 500* 501 DO 170 JITER = 1, MAXIT 502* 503* Check for too many iterations. 504* 505 IF( JITER.GT.MAXIT ) 506 $ GO TO 180 507* 508* Split the matrix if possible. 509* 510* Two tests: 511* 1: H(j,j-1)=0 or j=ILO 512* 2: T(j,j)=0 513* 514* Special case: j=ILAST 515* 516 IF( ILAST.EQ.ILO ) THEN 517 GO TO 60 518 ELSE 519 IF( ABS1( H( ILAST, ILAST-1 ) ).LE.MAX( SAFMIN, ULP*( 520 $ ABS1( H( ILAST, ILAST ) ) + ABS1( H( ILAST-1, ILAST-1 ) 521 $ ) ) ) ) THEN 522 H( ILAST, ILAST-1 ) = CZERO 523 GO TO 60 524 END IF 525 END IF 526* 527 IF( ABS( T( ILAST, ILAST ) ).LE.MAX( SAFMIN, ULP*( 528 $ ABS( T( ILAST - 1, ILAST ) ) + ABS( T( ILAST-1, ILAST-1 ) 529 $ ) ) ) ) THEN 530 T( ILAST, ILAST ) = CZERO 531 GO TO 50 532 END IF 533* 534* General case: j<ILAST 535* 536 DO 40 J = ILAST - 1, ILO, -1 537* 538* Test 1: for H(j,j-1)=0 or j=ILO 539* 540 IF( J.EQ.ILO ) THEN 541 ILAZRO = .TRUE. 542 ELSE 543 IF( ABS1( H( J, J-1 ) ).LE.MAX( SAFMIN, ULP*( 544 $ ABS1( H( J, J ) ) + ABS1( H( J-1, J-1 ) ) 545 $ ) ) ) THEN 546 H( J, J-1 ) = CZERO 547 ILAZRO = .TRUE. 548 ELSE 549 ILAZRO = .FALSE. 550 END IF 551 END IF 552* 553* Test 2: for T(j,j)=0 554* 555 TEMP = ABS ( T( J, J + 1 ) ) 556 IF ( J .GT. ILO ) 557 $ TEMP = TEMP + ABS ( T( J - 1, J ) ) 558 IF( ABS( T( J, J ) ).LT.MAX( SAFMIN,ULP*TEMP ) ) THEN 559 T( J, J ) = CZERO 560* 561* Test 1a: Check for 2 consecutive small subdiagonals in A 562* 563 ILAZR2 = .FALSE. 564 IF( .NOT.ILAZRO ) THEN 565 IF( ABS1( H( J, J-1 ) )*( ASCALE*ABS1( H( J+1, 566 $ J ) ) ).LE.ABS1( H( J, J ) )*( ASCALE*ATOL ) ) 567 $ ILAZR2 = .TRUE. 568 END IF 569* 570* If both tests pass (1 & 2), i.e., the leading diagonal 571* element of B in the block is zero, split a 1x1 block off 572* at the top. (I.e., at the J-th row/column) The leading 573* diagonal element of the remainder can also be zero, so 574* this may have to be done repeatedly. 575* 576 IF( ILAZRO .OR. ILAZR2 ) THEN 577 DO 20 JCH = J, ILAST - 1 578 CTEMP = H( JCH, JCH ) 579 CALL ZLARTG( CTEMP, H( JCH+1, JCH ), C, S, 580 $ H( JCH, JCH ) ) 581 H( JCH+1, JCH ) = CZERO 582 CALL ZROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH, 583 $ H( JCH+1, JCH+1 ), LDH, C, S ) 584 CALL ZROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT, 585 $ T( JCH+1, JCH+1 ), LDT, C, S ) 586 IF( ILQ ) 587 $ CALL ZROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1, 588 $ C, DCONJG( S ) ) 589 IF( ILAZR2 ) 590 $ H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C 591 ILAZR2 = .FALSE. 592 IF( ABS1( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN 593 IF( JCH+1.GE.ILAST ) THEN 594 GO TO 60 595 ELSE 596 IFIRST = JCH + 1 597 GO TO 70 598 END IF 599 END IF 600 T( JCH+1, JCH+1 ) = CZERO 601 20 CONTINUE 602 GO TO 50 603 ELSE 604* 605* Only test 2 passed -- chase the zero to T(ILAST,ILAST) 606* Then process as in the case T(ILAST,ILAST)=0 607* 608 DO 30 JCH = J, ILAST - 1 609 CTEMP = T( JCH, JCH+1 ) 610 CALL ZLARTG( CTEMP, T( JCH+1, JCH+1 ), C, S, 611 $ T( JCH, JCH+1 ) ) 612 T( JCH+1, JCH+1 ) = CZERO 613 IF( JCH.LT.ILASTM-1 ) 614 $ CALL ZROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT, 615 $ T( JCH+1, JCH+2 ), LDT, C, S ) 616 CALL ZROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH, 617 $ H( JCH+1, JCH-1 ), LDH, C, S ) 618 IF( ILQ ) 619 $ CALL ZROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1, 620 $ C, DCONJG( S ) ) 621 CTEMP = H( JCH+1, JCH ) 622 CALL ZLARTG( CTEMP, H( JCH+1, JCH-1 ), C, S, 623 $ H( JCH+1, JCH ) ) 624 H( JCH+1, JCH-1 ) = CZERO 625 CALL ZROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1, 626 $ H( IFRSTM, JCH-1 ), 1, C, S ) 627 CALL ZROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1, 628 $ T( IFRSTM, JCH-1 ), 1, C, S ) 629 IF( ILZ ) 630 $ CALL ZROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1, 631 $ C, S ) 632 30 CONTINUE 633 GO TO 50 634 END IF 635 ELSE IF( ILAZRO ) THEN 636* 637* Only test 1 passed -- work on J:ILAST 638* 639 IFIRST = J 640 GO TO 70 641 END IF 642* 643* Neither test passed -- try next J 644* 645 40 CONTINUE 646* 647* (Drop-through is "impossible") 648* 649 INFO = 2*N + 1 650 GO TO 210 651* 652* T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a 653* 1x1 block. 654* 655 50 CONTINUE 656 CTEMP = H( ILAST, ILAST ) 657 CALL ZLARTG( CTEMP, H( ILAST, ILAST-1 ), C, S, 658 $ H( ILAST, ILAST ) ) 659 H( ILAST, ILAST-1 ) = CZERO 660 CALL ZROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1, 661 $ H( IFRSTM, ILAST-1 ), 1, C, S ) 662 CALL ZROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1, 663 $ T( IFRSTM, ILAST-1 ), 1, C, S ) 664 IF( ILZ ) 665 $ CALL ZROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S ) 666* 667* H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHA and BETA 668* 669 60 CONTINUE 670 ABSB = ABS( T( ILAST, ILAST ) ) 671 IF( ABSB.GT.SAFMIN ) THEN 672 SIGNBC = DCONJG( T( ILAST, ILAST ) / ABSB ) 673 T( ILAST, ILAST ) = ABSB 674 IF( ILSCHR ) THEN 675 CALL ZSCAL( ILAST-IFRSTM, SIGNBC, T( IFRSTM, ILAST ), 1 ) 676 CALL ZSCAL( ILAST+1-IFRSTM, SIGNBC, H( IFRSTM, ILAST ), 677 $ 1 ) 678 ELSE 679 CALL ZSCAL( 1, SIGNBC, H( ILAST, ILAST ), 1 ) 680 END IF 681 IF( ILZ ) 682 $ CALL ZSCAL( N, SIGNBC, Z( 1, ILAST ), 1 ) 683 ELSE 684 T( ILAST, ILAST ) = CZERO 685 END IF 686 ALPHA( ILAST ) = H( ILAST, ILAST ) 687 BETA( ILAST ) = T( ILAST, ILAST ) 688* 689* Go to next block -- exit if finished. 690* 691 ILAST = ILAST - 1 692 IF( ILAST.LT.ILO ) 693 $ GO TO 190 694* 695* Reset counters 696* 697 IITER = 0 698 ESHIFT = CZERO 699 IF( .NOT.ILSCHR ) THEN 700 ILASTM = ILAST 701 IF( IFRSTM.GT.ILAST ) 702 $ IFRSTM = ILO 703 END IF 704 GO TO 160 705* 706* QZ step 707* 708* This iteration only involves rows/columns IFIRST:ILAST. We 709* assume IFIRST < ILAST, and that the diagonal of B is non-zero. 710* 711 70 CONTINUE 712 IITER = IITER + 1 713 IF( .NOT.ILSCHR ) THEN 714 IFRSTM = IFIRST 715 END IF 716* 717* Compute the Shift. 718* 719* At this point, IFIRST < ILAST, and the diagonal elements of 720* T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in 721* magnitude) 722* 723 IF( ( IITER / 10 )*10.NE.IITER ) THEN 724* 725* The Wilkinson shift (AEP p.512), i.e., the eigenvalue of 726* the bottom-right 2x2 block of A inv(B) which is nearest to 727* the bottom-right element. 728* 729* We factor B as U*D, where U has unit diagonals, and 730* compute (A*inv(D))*inv(U). 731* 732 U12 = ( BSCALE*T( ILAST-1, ILAST ) ) / 733 $ ( BSCALE*T( ILAST, ILAST ) ) 734 AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) / 735 $ ( BSCALE*T( ILAST-1, ILAST-1 ) ) 736 AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) / 737 $ ( BSCALE*T( ILAST-1, ILAST-1 ) ) 738 AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) / 739 $ ( BSCALE*T( ILAST, ILAST ) ) 740 AD22 = ( ASCALE*H( ILAST, ILAST ) ) / 741 $ ( BSCALE*T( ILAST, ILAST ) ) 742 ABI22 = AD22 - U12*AD21 743 ABI12 = AD12 - U12*AD11 744* 745 SHIFT = ABI22 746 CTEMP = SQRT( ABI12 )*SQRT( AD21 ) 747 TEMP = ABS1( CTEMP ) 748 IF( CTEMP.NE.ZERO ) THEN 749 X = HALF*( AD11-SHIFT ) 750 TEMP2 = ABS1( X ) 751 TEMP = MAX( TEMP, ABS1( X ) ) 752 Y = TEMP*SQRT( ( X / TEMP )**2+( CTEMP / TEMP )**2 ) 753 IF( TEMP2.GT.ZERO ) THEN 754 IF( DBLE( X / TEMP2 )*DBLE( Y )+ 755 $ DIMAG( X / TEMP2 )*DIMAG( Y ).LT.ZERO )Y = -Y 756 END IF 757 SHIFT = SHIFT - CTEMP*ZLADIV( CTEMP, ( X+Y ) ) 758 END IF 759 ELSE 760* 761* Exceptional shift. Chosen for no particularly good reason. 762* 763 IF( ( IITER / 20 )*20.EQ.IITER .AND. 764 $ BSCALE*ABS1(T( ILAST, ILAST )).GT.SAFMIN ) THEN 765 ESHIFT = ESHIFT + ( ASCALE*H( ILAST, 766 $ ILAST ) )/( BSCALE*T( ILAST, ILAST ) ) 767 ELSE 768 ESHIFT = ESHIFT + ( ASCALE*H( ILAST, 769 $ ILAST-1 ) )/( BSCALE*T( ILAST-1, ILAST-1 ) ) 770 END IF 771 SHIFT = ESHIFT 772 END IF 773* 774* Now check for two consecutive small subdiagonals. 775* 776 DO 80 J = ILAST - 1, IFIRST + 1, -1 777 ISTART = J 778 CTEMP = ASCALE*H( J, J ) - SHIFT*( BSCALE*T( J, J ) ) 779 TEMP = ABS1( CTEMP ) 780 TEMP2 = ASCALE*ABS1( H( J+1, J ) ) 781 TEMPR = MAX( TEMP, TEMP2 ) 782 IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN 783 TEMP = TEMP / TEMPR 784 TEMP2 = TEMP2 / TEMPR 785 END IF 786 IF( ABS1( H( J, J-1 ) )*TEMP2.LE.TEMP*ATOL ) 787 $ GO TO 90 788 80 CONTINUE 789* 790 ISTART = IFIRST 791 CTEMP = ASCALE*H( IFIRST, IFIRST ) - 792 $ SHIFT*( BSCALE*T( IFIRST, IFIRST ) ) 793 90 CONTINUE 794* 795* Do an implicit-shift QZ sweep. 796* 797* Initial Q 798* 799 CTEMP2 = ASCALE*H( ISTART+1, ISTART ) 800 CALL ZLARTG( CTEMP, CTEMP2, C, S, CTEMP3 ) 801* 802* Sweep 803* 804 DO 150 J = ISTART, ILAST - 1 805 IF( J.GT.ISTART ) THEN 806 CTEMP = H( J, J-1 ) 807 CALL ZLARTG( CTEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) ) 808 H( J+1, J-1 ) = CZERO 809 END IF 810* 811 DO 100 JC = J, ILASTM 812 CTEMP = C*H( J, JC ) + S*H( J+1, JC ) 813 H( J+1, JC ) = -DCONJG( S )*H( J, JC ) + C*H( J+1, JC ) 814 H( J, JC ) = CTEMP 815 CTEMP2 = C*T( J, JC ) + S*T( J+1, JC ) 816 T( J+1, JC ) = -DCONJG( S )*T( J, JC ) + C*T( J+1, JC ) 817 T( J, JC ) = CTEMP2 818 100 CONTINUE 819 IF( ILQ ) THEN 820 DO 110 JR = 1, N 821 CTEMP = C*Q( JR, J ) + DCONJG( S )*Q( JR, J+1 ) 822 Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 ) 823 Q( JR, J ) = CTEMP 824 110 CONTINUE 825 END IF 826* 827 CTEMP = T( J+1, J+1 ) 828 CALL ZLARTG( CTEMP, T( J+1, J ), C, S, T( J+1, J+1 ) ) 829 T( J+1, J ) = CZERO 830* 831 DO 120 JR = IFRSTM, MIN( J+2, ILAST ) 832 CTEMP = C*H( JR, J+1 ) + S*H( JR, J ) 833 H( JR, J ) = -DCONJG( S )*H( JR, J+1 ) + C*H( JR, J ) 834 H( JR, J+1 ) = CTEMP 835 120 CONTINUE 836 DO 130 JR = IFRSTM, J 837 CTEMP = C*T( JR, J+1 ) + S*T( JR, J ) 838 T( JR, J ) = -DCONJG( S )*T( JR, J+1 ) + C*T( JR, J ) 839 T( JR, J+1 ) = CTEMP 840 130 CONTINUE 841 IF( ILZ ) THEN 842 DO 140 JR = 1, N 843 CTEMP = C*Z( JR, J+1 ) + S*Z( JR, J ) 844 Z( JR, J ) = -DCONJG( S )*Z( JR, J+1 ) + C*Z( JR, J ) 845 Z( JR, J+1 ) = CTEMP 846 140 CONTINUE 847 END IF 848 150 CONTINUE 849* 850 160 CONTINUE 851* 852 170 CONTINUE 853* 854* Drop-through = non-convergence 855* 856 180 CONTINUE 857 INFO = ILAST 858 GO TO 210 859* 860* Successful completion of all QZ steps 861* 862 190 CONTINUE 863* 864* Set Eigenvalues 1:ILO-1 865* 866 DO 200 J = 1, ILO - 1 867 ABSB = ABS( T( J, J ) ) 868 IF( ABSB.GT.SAFMIN ) THEN 869 SIGNBC = DCONJG( T( J, J ) / ABSB ) 870 T( J, J ) = ABSB 871 IF( ILSCHR ) THEN 872 CALL ZSCAL( J-1, SIGNBC, T( 1, J ), 1 ) 873 CALL ZSCAL( J, SIGNBC, H( 1, J ), 1 ) 874 ELSE 875 CALL ZSCAL( 1, SIGNBC, H( J, J ), 1 ) 876 END IF 877 IF( ILZ ) 878 $ CALL ZSCAL( N, SIGNBC, Z( 1, J ), 1 ) 879 ELSE 880 T( J, J ) = CZERO 881 END IF 882 ALPHA( J ) = H( J, J ) 883 BETA( J ) = T( J, J ) 884 200 CONTINUE 885* 886* Normal Termination 887* 888 INFO = 0 889* 890* Exit (other than argument error) -- return optimal workspace size 891* 892 210 CONTINUE 893 WORK( 1 ) = DCMPLX( N ) 894 RETURN 895* 896* End of ZHGEQZ 897* 898 END 899