1*> \brief \b ZLAQZ0 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZLAQZ0 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaqz0.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaqz0.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaqz0.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A, LDA, B, 22* $ LDB, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, REC, 23* $ INFO ) 24* IMPLICIT NONE 25* 26* Arguments 27* CHARACTER, INTENT( IN ) :: WANTS, WANTQ, WANTZ 28* INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK, 29* $ REC 30* INTEGER, INTENT( OUT ) :: INFO 31* COMPLEX*16, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, 32* $ * ), Z( LDZ, * ), ALPHA( * ), BETA( * ), WORK( * ) 33* DOUBLE PRECISION, INTENT( OUT ) :: RWORK( * ) 34* .. 35* 36* 37*> \par Purpose: 38* ============= 39*> 40*> \verbatim 41*> 42*> ZLAQZ0 computes the eigenvalues of a real matrix pair (H,T), 43*> where H is an upper Hessenberg matrix and T is upper triangular, 44*> using the double-shift QZ method. 45*> Matrix pairs of this type are produced by the reduction to 46*> generalized upper Hessenberg form of a real 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, P and S are an upper triangular 58*> matrices. 59*> 60*> Optionally, the unitary matrix Q from the generalized Schur 61*> factorization may be postmultiplied into an input matrix Q1, and the 62*> unitary matrix Z may be postmultiplied into an input matrix Z1. 63*> If Q1 and Z1 are the unitary matrices from ZGGHRD that reduced 64*> the matrix pair (A,B) to generalized upper Hessenberg form, then the 65*> output matrices Q1*Q and Z1*Z are the unitary factors from the 66*> generalized Schur factorization of (A,B): 67*> 68*> A = (Q1*Q)*S*(Z1*Z)**H, B = (Q1*Q)*P*(Z1*Z)**H. 69*> 70*> To avoid overflow, eigenvalues of the matrix pair (H,T) (equivalently, 71*> of (A,B)) are computed as a pair of values (alpha,beta), where alpha is 72*> complex and beta real. 73*> If beta is nonzero, lambda = alpha / beta is an eigenvalue of the 74*> generalized nonsymmetric eigenvalue problem (GNEP) 75*> A*x = lambda*B*x 76*> and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the 77*> alternate form of the GNEP 78*> mu*A*y = B*y. 79*> Eigenvalues can be read directly from the generalized Schur 80*> form: 81*> alpha = S(i,i), beta = P(i,i). 82*> 83*> Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix 84*> Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973), 85*> pp. 241--256. 86*> 87*> Ref: B. Kagstrom, D. Kressner, "Multishift Variants of the QZ 88*> Algorithm with Aggressive Early Deflation", SIAM J. Numer. 89*> Anal., 29(2006), pp. 199--227. 90*> 91*> Ref: T. Steel, D. Camps, K. Meerbergen, R. Vandebril "A multishift, 92*> multipole rational QZ method with agressive early deflation" 93*> \endverbatim 94* 95* Arguments: 96* ========== 97* 98*> \param[in] WANTS 99*> \verbatim 100*> WANTS is CHARACTER*1 101*> = 'E': Compute eigenvalues only; 102*> = 'S': Compute eigenvalues and the Schur form. 103*> \endverbatim 104*> 105*> \param[in] WANTQ 106*> \verbatim 107*> WANTQ is CHARACTER*1 108*> = 'N': Left Schur vectors (Q) are not computed; 109*> = 'I': Q is initialized to the unit matrix and the matrix Q 110*> of left Schur vectors of (A,B) is returned; 111*> = 'V': Q must contain an unitary matrix Q1 on entry and 112*> the product Q1*Q is returned. 113*> \endverbatim 114*> 115*> \param[in] WANTZ 116*> \verbatim 117*> WANTZ is CHARACTER*1 118*> = 'N': Right Schur vectors (Z) are not computed; 119*> = 'I': Z is initialized to the unit matrix and the matrix Z 120*> of right Schur vectors of (A,B) is returned; 121*> = 'V': Z must contain an unitary matrix Z1 on entry and 122*> the product Z1*Z is returned. 123*> \endverbatim 124*> 125*> \param[in] N 126*> \verbatim 127*> N is INTEGER 128*> The order of the matrices A, B, Q, and Z. N >= 0. 129*> \endverbatim 130*> 131*> \param[in] ILO 132*> \verbatim 133*> ILO is INTEGER 134*> \endverbatim 135*> 136*> \param[in] IHI 137*> \verbatim 138*> IHI is INTEGER 139*> ILO and IHI mark the rows and columns of A which are in 140*> Hessenberg form. It is assumed that A is already upper 141*> triangular in rows and columns 1:ILO-1 and IHI+1:N. 142*> If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0. 143*> \endverbatim 144*> 145*> \param[in,out] A 146*> \verbatim 147*> A is COMPLEX*16 array, dimension (LDA, N) 148*> On entry, the N-by-N upper Hessenberg matrix A. 149*> On exit, if JOB = 'S', A contains the upper triangular 150*> matrix S from the generalized Schur factorization. 151*> If JOB = 'E', the diagonal blocks of A match those of S, but 152*> the rest of A is unspecified. 153*> \endverbatim 154*> 155*> \param[in] LDA 156*> \verbatim 157*> LDA is INTEGER 158*> The leading dimension of the array A. LDA >= max( 1, N ). 159*> \endverbatim 160*> 161*> \param[in,out] B 162*> \verbatim 163*> B is COMPLEX*16 array, dimension (LDB, N) 164*> On entry, the N-by-N upper triangular matrix B. 165*> On exit, if JOB = 'S', B contains the upper triangular 166*> matrix P from the generalized Schur factorization; 167*> If JOB = 'E', the diagonal blocks of B match those of P, but 168*> the rest of B is unspecified. 169*> \endverbatim 170*> 171*> \param[in] LDB 172*> \verbatim 173*> LDB is INTEGER 174*> The leading dimension of the array B. LDB >= max( 1, N ). 175*> \endverbatim 176*> 177*> \param[out] ALPHA 178*> \verbatim 179*> ALPHA is COMPLEX*16 array, dimension (N) 180*> Each scalar alpha defining an eigenvalue 181*> of GNEP. 182*> \endverbatim 183*> 184*> \param[out] BETA 185*> \verbatim 186*> BETA is COMPLEX*16 array, dimension (N) 187*> The scalars beta that define the eigenvalues of GNEP. 188*> Together, the quantities alpha = ALPHA(j) and 189*> beta = BETA(j) represent the j-th eigenvalue of the matrix 190*> pair (A,B), in one of the forms lambda = alpha/beta or 191*> mu = beta/alpha. Since either lambda or mu may overflow, 192*> they should not, in general, be computed. 193*> \endverbatim 194*> 195*> \param[in,out] Q 196*> \verbatim 197*> Q is COMPLEX*16 array, dimension (LDQ, N) 198*> On entry, if COMPQ = 'V', the unitary matrix Q1 used in 199*> the reduction of (A,B) to generalized Hessenberg form. 200*> On exit, if COMPQ = 'I', the unitary matrix of left Schur 201*> vectors of (A,B), and if COMPQ = 'V', the unitary matrix 202*> of left Schur vectors of (A,B). 203*> Not referenced if COMPQ = 'N'. 204*> \endverbatim 205*> 206*> \param[in] LDQ 207*> \verbatim 208*> LDQ is INTEGER 209*> The leading dimension of the array Q. LDQ >= 1. 210*> If COMPQ='V' or 'I', then LDQ >= N. 211*> \endverbatim 212*> 213*> \param[in,out] Z 214*> \verbatim 215*> Z is COMPLEX*16 array, dimension (LDZ, N) 216*> On entry, if COMPZ = 'V', the unitary matrix Z1 used in 217*> the reduction of (A,B) to generalized Hessenberg form. 218*> On exit, if COMPZ = 'I', the unitary matrix of 219*> right Schur vectors of (H,T), and if COMPZ = 'V', the 220*> unitary matrix of right Schur vectors of (A,B). 221*> Not referenced if COMPZ = 'N'. 222*> \endverbatim 223*> 224*> \param[in] LDZ 225*> \verbatim 226*> LDZ is INTEGER 227*> The leading dimension of the array Z. LDZ >= 1. 228*> If COMPZ='V' or 'I', then LDZ >= N. 229*> \endverbatim 230*> 231*> \param[out] WORK 232*> \verbatim 233*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) 234*> On exit, if INFO >= 0, WORK(1) returns the optimal LWORK. 235*> \endverbatim 236*> 237*> \param[out] RWORK 238*> \verbatim 239*> RWORK is DOUBLE PRECISION array, dimension (N) 240*> \endverbatim 241*> 242*> \param[in] LWORK 243*> \verbatim 244*> LWORK is INTEGER 245*> The dimension of the array WORK. LWORK >= max(1,N). 246*> 247*> If LWORK = -1, then a workspace query is assumed; the routine 248*> only calculates the optimal size of the WORK array, returns 249*> this value as the first entry of the WORK array, and no error 250*> message related to LWORK is issued by XERBLA. 251*> \endverbatim 252*> 253*> \param[in] REC 254*> \verbatim 255*> REC is INTEGER 256*> REC indicates the current recursion level. Should be set 257*> to 0 on first call. 258*> \endverbatim 259*> 260*> \param[out] INFO 261*> \verbatim 262*> INFO is INTEGER 263*> = 0: successful exit 264*> < 0: if INFO = -i, the i-th argument had an illegal value 265*> = 1,...,N: the QZ iteration did not converge. (A,B) is not 266*> in Schur form, but ALPHA(i) and 267*> BETA(i), i=INFO+1,...,N should be correct. 268*> \endverbatim 269* 270* Authors: 271* ======== 272* 273*> \author Thijs Steel, KU Leuven 274* 275*> \date May 2020 276* 277*> \ingroup complex16GEcomputational 278*> 279* ===================================================================== 280 RECURSIVE SUBROUTINE ZLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A, 281 $ LDA, B, LDB, ALPHA, BETA, Q, LDQ, Z, 282 $ LDZ, WORK, LWORK, RWORK, REC, 283 $ INFO ) 284 IMPLICIT NONE 285 286* Arguments 287 CHARACTER, INTENT( IN ) :: WANTS, WANTQ, WANTZ 288 INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK, 289 $ REC 290 INTEGER, INTENT( OUT ) :: INFO 291 COMPLEX*16, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, 292 $ * ), Z( LDZ, * ), ALPHA( * ), BETA( * ), WORK( * ) 293 DOUBLE PRECISION, INTENT( OUT ) :: RWORK( * ) 294 295* Parameters 296 COMPLEX*16 CZERO, CONE 297 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), CONE = ( 1.0D+0, 298 $ 0.0D+0 ) ) 299 DOUBLE PRECISION :: ZERO, ONE, HALF 300 PARAMETER( ZERO = 0.0D0, ONE = 1.0D0, HALF = 0.5D0 ) 301 302* Local scalars 303 DOUBLE PRECISION :: SMLNUM, ULP, SAFMIN, SAFMAX, C1, TEMPR 304 COMPLEX*16 :: ESHIFT, S1, TEMP 305 INTEGER :: ISTART, ISTOP, IITER, MAXIT, ISTART2, K, LD, NSHIFTS, 306 $ NBLOCK, NW, NMIN, NIBBLE, N_UNDEFLATED, N_DEFLATED, 307 $ NS, SWEEP_INFO, SHIFTPOS, LWORKREQ, K2, ISTARTM, 308 $ ISTOPM, IWANTS, IWANTQ, IWANTZ, NORM_INFO, AED_INFO, 309 $ NWR, NBR, NSR, ITEMP1, ITEMP2, RCOST 310 LOGICAL :: ILSCHUR, ILQ, ILZ 311 CHARACTER :: JBCMPZ*3 312 313* External Functions 314 EXTERNAL :: XERBLA, ZHGEQZ, ZLAQZ2, ZLAQZ3, ZLASET, DLABAD, 315 $ ZLARTG, ZROT 316 DOUBLE PRECISION, EXTERNAL :: DLAMCH 317 LOGICAL, EXTERNAL :: LSAME 318 INTEGER, EXTERNAL :: ILAENV 319 320* 321* Decode wantS,wantQ,wantZ 322* 323 IF( LSAME( WANTS, 'E' ) ) THEN 324 ILSCHUR = .FALSE. 325 IWANTS = 1 326 ELSE IF( LSAME( WANTS, 'S' ) ) THEN 327 ILSCHUR = .TRUE. 328 IWANTS = 2 329 ELSE 330 IWANTS = 0 331 END IF 332 333 IF( LSAME( WANTQ, 'N' ) ) THEN 334 ILQ = .FALSE. 335 IWANTQ = 1 336 ELSE IF( LSAME( WANTQ, 'V' ) ) THEN 337 ILQ = .TRUE. 338 IWANTQ = 2 339 ELSE IF( LSAME( WANTQ, 'I' ) ) THEN 340 ILQ = .TRUE. 341 IWANTQ = 3 342 ELSE 343 IWANTQ = 0 344 END IF 345 346 IF( LSAME( WANTZ, 'N' ) ) THEN 347 ILZ = .FALSE. 348 IWANTZ = 1 349 ELSE IF( LSAME( WANTZ, 'V' ) ) THEN 350 ILZ = .TRUE. 351 IWANTZ = 2 352 ELSE IF( LSAME( WANTZ, 'I' ) ) THEN 353 ILZ = .TRUE. 354 IWANTZ = 3 355 ELSE 356 IWANTZ = 0 357 END IF 358* 359* Check Argument Values 360* 361 INFO = 0 362 IF( IWANTS.EQ.0 ) THEN 363 INFO = -1 364 ELSE IF( IWANTQ.EQ.0 ) THEN 365 INFO = -2 366 ELSE IF( IWANTZ.EQ.0 ) THEN 367 INFO = -3 368 ELSE IF( N.LT.0 ) THEN 369 INFO = -4 370 ELSE IF( ILO.LT.1 ) THEN 371 INFO = -5 372 ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN 373 INFO = -6 374 ELSE IF( LDA.LT.N ) THEN 375 INFO = -8 376 ELSE IF( LDB.LT.N ) THEN 377 INFO = -10 378 ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN 379 INFO = -15 380 ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN 381 INFO = -17 382 END IF 383 IF( INFO.NE.0 ) THEN 384 CALL XERBLA( 'ZLAQZ0', -INFO ) 385 RETURN 386 END IF 387 388* 389* Quick return if possible 390* 391 IF( N.LE.0 ) THEN 392 WORK( 1 ) = DBLE( 1 ) 393 RETURN 394 END IF 395 396* 397* Get the parameters 398* 399 JBCMPZ( 1:1 ) = WANTS 400 JBCMPZ( 2:2 ) = WANTQ 401 JBCMPZ( 3:3 ) = WANTZ 402 403 NMIN = ILAENV( 12, 'ZLAQZ0', JBCMPZ, N, ILO, IHI, LWORK ) 404 405 NWR = ILAENV( 13, 'ZLAQZ0', JBCMPZ, N, ILO, IHI, LWORK ) 406 NWR = MAX( 2, NWR ) 407 NWR = MIN( IHI-ILO+1, ( N-1 ) / 3, NWR ) 408 409 NIBBLE = ILAENV( 14, 'ZLAQZ0', JBCMPZ, N, ILO, IHI, LWORK ) 410 411 NSR = ILAENV( 15, 'ZLAQZ0', JBCMPZ, N, ILO, IHI, LWORK ) 412 NSR = MIN( NSR, ( N+6 ) / 9, IHI-ILO ) 413 NSR = MAX( 2, NSR-MOD( NSR, 2 ) ) 414 415 RCOST = ILAENV( 17, 'ZLAQZ0', JBCMPZ, N, ILO, IHI, LWORK ) 416 ITEMP1 = INT( NSR/SQRT( 1+2*NSR/( DBLE( RCOST )/100*N ) ) ) 417 ITEMP1 = ( ( ITEMP1-1 )/4 )*4+4 418 NBR = NSR+ITEMP1 419 420 IF( N .LT. NMIN .OR. REC .GE. 2 ) THEN 421 CALL ZHGEQZ( WANTS, WANTQ, WANTZ, N, ILO, IHI, A, LDA, B, LDB, 422 $ ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, 423 $ INFO ) 424 RETURN 425 END IF 426 427* 428* Find out required workspace 429* 430 431* Workspace query to ZLAQZ2 432 NW = MAX( NWR, NMIN ) 433 CALL ZLAQZ2( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW, A, LDA, B, LDB, 434 $ Q, LDQ, Z, LDZ, N_UNDEFLATED, N_DEFLATED, ALPHA, 435 $ BETA, WORK, NW, WORK, NW, WORK, -1, RWORK, REC, 436 $ AED_INFO ) 437 ITEMP1 = INT( WORK( 1 ) ) 438* Workspace query to ZLAQZ3 439 CALL ZLAQZ3( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSR, NBR, ALPHA, 440 $ BETA, A, LDA, B, LDB, Q, LDQ, Z, LDZ, WORK, NBR, 441 $ WORK, NBR, WORK, -1, SWEEP_INFO ) 442 ITEMP2 = INT( WORK( 1 ) ) 443 444 LWORKREQ = MAX( ITEMP1+2*NW**2, ITEMP2+2*NBR**2 ) 445 IF ( LWORK .EQ.-1 ) THEN 446 WORK( 1 ) = DBLE( LWORKREQ ) 447 RETURN 448 ELSE IF ( LWORK .LT. LWORKREQ ) THEN 449 INFO = -19 450 END IF 451 IF( INFO.NE.0 ) THEN 452 CALL XERBLA( 'ZLAQZ0', INFO ) 453 RETURN 454 END IF 455* 456* Initialize Q and Z 457* 458 IF( IWANTQ.EQ.3 ) CALL ZLASET( 'FULL', N, N, CZERO, CONE, Q, 459 $ LDQ ) 460 IF( IWANTZ.EQ.3 ) CALL ZLASET( 'FULL', N, N, CZERO, CONE, Z, 461 $ LDZ ) 462 463* Get machine constants 464 SAFMIN = DLAMCH( 'SAFE MINIMUM' ) 465 SAFMAX = ONE/SAFMIN 466 CALL DLABAD( SAFMIN, SAFMAX ) 467 ULP = DLAMCH( 'PRECISION' ) 468 SMLNUM = SAFMIN*( DBLE( N )/ULP ) 469 470 ISTART = ILO 471 ISTOP = IHI 472 MAXIT = 30*( IHI-ILO+1 ) 473 LD = 0 474 475 DO IITER = 1, MAXIT 476 IF( IITER .GE. MAXIT ) THEN 477 INFO = ISTOP+1 478 GOTO 80 479 END IF 480 IF ( ISTART+1 .GE. ISTOP ) THEN 481 ISTOP = ISTART 482 EXIT 483 END IF 484 485* Check deflations at the end 486 IF ( ABS( A( ISTOP, ISTOP-1 ) ) .LE. MAX( SMLNUM, 487 $ ULP*( ABS( A( ISTOP, ISTOP ) )+ABS( A( ISTOP-1, 488 $ ISTOP-1 ) ) ) ) ) THEN 489 A( ISTOP, ISTOP-1 ) = CZERO 490 ISTOP = ISTOP-1 491 LD = 0 492 ESHIFT = CZERO 493 END IF 494* Check deflations at the start 495 IF ( ABS( A( ISTART+1, ISTART ) ) .LE. MAX( SMLNUM, 496 $ ULP*( ABS( A( ISTART, ISTART ) )+ABS( A( ISTART+1, 497 $ ISTART+1 ) ) ) ) ) THEN 498 A( ISTART+1, ISTART ) = CZERO 499 ISTART = ISTART+1 500 LD = 0 501 ESHIFT = CZERO 502 END IF 503 504 IF ( ISTART+1 .GE. ISTOP ) THEN 505 EXIT 506 END IF 507 508* Check interior deflations 509 ISTART2 = ISTART 510 DO K = ISTOP, ISTART+1, -1 511 IF ( ABS( A( K, K-1 ) ) .LE. MAX( SMLNUM, ULP*( ABS( A( K, 512 $ K ) )+ABS( A( K-1, K-1 ) ) ) ) ) THEN 513 A( K, K-1 ) = CZERO 514 ISTART2 = K 515 EXIT 516 END IF 517 END DO 518 519* Get range to apply rotations to 520 IF ( ILSCHUR ) THEN 521 ISTARTM = 1 522 ISTOPM = N 523 ELSE 524 ISTARTM = ISTART2 525 ISTOPM = ISTOP 526 END IF 527 528* Check infinite eigenvalues, this is done without blocking so might 529* slow down the method when many infinite eigenvalues are present 530 K = ISTOP 531 DO WHILE ( K.GE.ISTART2 ) 532 TEMPR = ZERO 533 IF( K .LT. ISTOP ) THEN 534 TEMPR = TEMPR+ABS( B( K, K+1 ) ) 535 END IF 536 IF( K .GT. ISTART2 ) THEN 537 TEMPR = TEMPR+ABS( B( K-1, K ) ) 538 END IF 539 540 IF( ABS( B( K, K ) ) .LT. MAX( SMLNUM, ULP*TEMPR ) ) THEN 541* A diagonal element of B is negligable, move it 542* to the top and deflate it 543 544 DO K2 = K, ISTART2+1, -1 545 CALL ZLARTG( B( K2-1, K2 ), B( K2-1, K2-1 ), C1, S1, 546 $ TEMP ) 547 B( K2-1, K2 ) = TEMP 548 B( K2-1, K2-1 ) = CZERO 549 550 CALL ZROT( K2-2-ISTARTM+1, B( ISTARTM, K2 ), 1, 551 $ B( ISTARTM, K2-1 ), 1, C1, S1 ) 552 CALL ZROT( MIN( K2+1, ISTOP )-ISTARTM+1, A( ISTARTM, 553 $ K2 ), 1, A( ISTARTM, K2-1 ), 1, C1, S1 ) 554 IF ( ILZ ) THEN 555 CALL ZROT( N, Z( 1, K2 ), 1, Z( 1, K2-1 ), 1, C1, 556 $ S1 ) 557 END IF 558 559 IF( K2.LT.ISTOP ) THEN 560 CALL ZLARTG( A( K2, K2-1 ), A( K2+1, K2-1 ), C1, 561 $ S1, TEMP ) 562 A( K2, K2-1 ) = TEMP 563 A( K2+1, K2-1 ) = CZERO 564 565 CALL ZROT( ISTOPM-K2+1, A( K2, K2 ), LDA, A( K2+1, 566 $ K2 ), LDA, C1, S1 ) 567 CALL ZROT( ISTOPM-K2+1, B( K2, K2 ), LDB, B( K2+1, 568 $ K2 ), LDB, C1, S1 ) 569 IF( ILQ ) THEN 570 CALL ZROT( N, Q( 1, K2 ), 1, Q( 1, K2+1 ), 1, 571 $ C1, DCONJG( S1 ) ) 572 END IF 573 END IF 574 575 END DO 576 577 IF( ISTART2.LT.ISTOP )THEN 578 CALL ZLARTG( A( ISTART2, ISTART2 ), A( ISTART2+1, 579 $ ISTART2 ), C1, S1, TEMP ) 580 A( ISTART2, ISTART2 ) = TEMP 581 A( ISTART2+1, ISTART2 ) = CZERO 582 583 CALL ZROT( ISTOPM-( ISTART2+1 )+1, A( ISTART2, 584 $ ISTART2+1 ), LDA, A( ISTART2+1, 585 $ ISTART2+1 ), LDA, C1, S1 ) 586 CALL ZROT( ISTOPM-( ISTART2+1 )+1, B( ISTART2, 587 $ ISTART2+1 ), LDB, B( ISTART2+1, 588 $ ISTART2+1 ), LDB, C1, S1 ) 589 IF( ILQ ) THEN 590 CALL ZROT( N, Q( 1, ISTART2 ), 1, Q( 1, 591 $ ISTART2+1 ), 1, C1, DCONJG( S1 ) ) 592 END IF 593 END IF 594 595 ISTART2 = ISTART2+1 596 597 END IF 598 K = K-1 599 END DO 600 601* istart2 now points to the top of the bottom right 602* unreduced Hessenberg block 603 IF ( ISTART2 .GE. ISTOP ) THEN 604 ISTOP = ISTART2-1 605 LD = 0 606 ESHIFT = CZERO 607 CYCLE 608 END IF 609 610 NW = NWR 611 NSHIFTS = NSR 612 NBLOCK = NBR 613 614 IF ( ISTOP-ISTART2+1 .LT. NMIN ) THEN 615* Setting nw to the size of the subblock will make AED deflate 616* all the eigenvalues. This is slightly more efficient than just 617* using qz_small because the off diagonal part gets updated via BLAS. 618 IF ( ISTOP-ISTART+1 .LT. NMIN ) THEN 619 NW = ISTOP-ISTART+1 620 ISTART2 = ISTART 621 ELSE 622 NW = ISTOP-ISTART2+1 623 END IF 624 END IF 625 626* 627* Time for AED 628* 629 CALL ZLAQZ2( ILSCHUR, ILQ, ILZ, N, ISTART2, ISTOP, NW, A, LDA, 630 $ B, LDB, Q, LDQ, Z, LDZ, N_UNDEFLATED, N_DEFLATED, 631 $ ALPHA, BETA, WORK, NW, WORK( NW**2+1 ), NW, 632 $ WORK( 2*NW**2+1 ), LWORK-2*NW**2, RWORK, REC, 633 $ AED_INFO ) 634 635 IF ( N_DEFLATED > 0 ) THEN 636 ISTOP = ISTOP-N_DEFLATED 637 LD = 0 638 ESHIFT = CZERO 639 END IF 640 641 IF ( 100*N_DEFLATED > NIBBLE*( N_DEFLATED+N_UNDEFLATED ) .OR. 642 $ ISTOP-ISTART2+1 .LT. NMIN ) THEN 643* AED has uncovered many eigenvalues. Skip a QZ sweep and run 644* AED again. 645 CYCLE 646 END IF 647 648 LD = LD+1 649 650 NS = MIN( NSHIFTS, ISTOP-ISTART2 ) 651 NS = MIN( NS, N_UNDEFLATED ) 652 SHIFTPOS = ISTOP-N_DEFLATED-N_UNDEFLATED+1 653 654 IF ( MOD( LD, 6 ) .EQ. 0 ) THEN 655* 656* Exceptional shift. Chosen for no particularly good reason. 657* 658 IF( ( DBLE( MAXIT )*SAFMIN )*ABS( A( ISTOP, 659 $ ISTOP-1 ) ).LT.ABS( A( ISTOP-1, ISTOP-1 ) ) ) THEN 660 ESHIFT = A( ISTOP, ISTOP-1 )/B( ISTOP-1, ISTOP-1 ) 661 ELSE 662 ESHIFT = ESHIFT+CONE/( SAFMIN*DBLE( MAXIT ) ) 663 END IF 664 ALPHA( SHIFTPOS ) = CONE 665 BETA( SHIFTPOS ) = ESHIFT 666 NS = 1 667 END IF 668 669* 670* Time for a QZ sweep 671* 672 CALL ZLAQZ3( ILSCHUR, ILQ, ILZ, N, ISTART2, ISTOP, NS, NBLOCK, 673 $ ALPHA( SHIFTPOS ), BETA( SHIFTPOS ), A, LDA, B, 674 $ LDB, Q, LDQ, Z, LDZ, WORK, NBLOCK, WORK( NBLOCK** 675 $ 2+1 ), NBLOCK, WORK( 2*NBLOCK**2+1 ), 676 $ LWORK-2*NBLOCK**2, SWEEP_INFO ) 677 678 END DO 679 680* 681* Call ZHGEQZ to normalize the eigenvalue blocks and set the eigenvalues 682* If all the eigenvalues have been found, ZHGEQZ will not do any iterations 683* and only normalize the blocks. In case of a rare convergence failure, 684* the single shift might perform better. 685* 686 80 CALL ZHGEQZ( WANTS, WANTQ, WANTZ, N, ILO, IHI, A, LDA, B, LDB, 687 $ ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, 688 $ NORM_INFO ) 689 690 INFO = NORM_INFO 691 692 END SUBROUTINE 693