1*> \brief \b CLAQZ0 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download CLAQZ0 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/CLAQZ0.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/CLAQZ0.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/CLAQZ0.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE CLAQZ0( 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, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 32* $ Z( LDZ, * ), ALPHA( * ), BETA( * ), WORK( * ) 33* REAL, INTENT( OUT ) :: RWORK( * ) 34* .. 35* 36* 37*> \par Purpose: 38* ============= 39*> 40*> \verbatim 41*> 42*> CLAQZ0 computes the eigenvalues of a 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 matrix pair (A,B): 47*> 48*> A = Q1*H*Z1**H, B = Q1*T*Z1**H, 49*> 50*> as computed by CGGHRD. 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 CGGHRD 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 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 of A matches that 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 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 of B matches that 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 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 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 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 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 array, dimension (MAX(1,LWORK)) 234*> On exit, if INFO >= 0, WORK(1) returns the optimal LWORK. 235*> \endverbatim 236*> 237*> \param[in] LWORK 238*> \verbatim 239*> LWORK is INTEGER 240*> The dimension of the array WORK. LWORK >= max(1,N). 241*> 242*> If LWORK = -1, then a workspace query is assumed; the routine 243*> only calculates the optimal size of the WORK array, returns 244*> this value as the first entry of the WORK array, and no error 245*> message related to LWORK is issued by XERBLA. 246*> \endverbatim 247*> 248*> \param[out] RWORK 249*> \verbatim 250*> RWORK is REAL array, dimension (N) 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 complexGEcomputational 278*> 279* ===================================================================== 280 RECURSIVE SUBROUTINE CLAQZ0( 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, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 292 $ Z( LDZ, * ), ALPHA( * ), BETA( * ), WORK( * ) 293 REAL, INTENT( OUT ) :: RWORK( * ) 294 295* Parameters 296 COMPLEX CZERO, CONE 297 PARAMETER ( CZERO = ( 0.0, 0.0 ), CONE = ( 1.0, 0.0 ) ) 298 REAL :: ZERO, ONE, HALF 299 PARAMETER( ZERO = 0.0, ONE = 1.0, HALF = 0.5 ) 300 301* Local scalars 302 REAL :: SMLNUM, ULP, SAFMIN, SAFMAX, C1, TEMPR 303 COMPLEX :: ESHIFT, S1, TEMP 304 INTEGER :: ISTART, ISTOP, IITER, MAXIT, ISTART2, K, LD, NSHIFTS, 305 $ NBLOCK, NW, NMIN, NIBBLE, N_UNDEFLATED, N_DEFLATED, 306 $ NS, SWEEP_INFO, SHIFTPOS, LWORKREQ, K2, ISTARTM, 307 $ ISTOPM, IWANTS, IWANTQ, IWANTZ, NORM_INFO, AED_INFO, 308 $ NWR, NBR, NSR, ITEMP1, ITEMP2, RCOST 309 LOGICAL :: ILSCHUR, ILQ, ILZ 310 CHARACTER :: JBCMPZ*3 311 312* External Functions 313 EXTERNAL :: XERBLA, CHGEQZ, CLAQZ2, CLAQZ3, CLASET, SLABAD, 314 $ CLARTG, CROT 315 REAL, EXTERNAL :: SLAMCH 316 LOGICAL, EXTERNAL :: LSAME 317 INTEGER, EXTERNAL :: ILAENV 318 319* 320* Decode wantS,wantQ,wantZ 321* 322 IF( LSAME( WANTS, 'E' ) ) THEN 323 ILSCHUR = .FALSE. 324 IWANTS = 1 325 ELSE IF( LSAME( WANTS, 'S' ) ) THEN 326 ILSCHUR = .TRUE. 327 IWANTS = 2 328 ELSE 329 IWANTS = 0 330 END IF 331 332 IF( LSAME( WANTQ, 'N' ) ) THEN 333 ILQ = .FALSE. 334 IWANTQ = 1 335 ELSE IF( LSAME( WANTQ, 'V' ) ) THEN 336 ILQ = .TRUE. 337 IWANTQ = 2 338 ELSE IF( LSAME( WANTQ, 'I' ) ) THEN 339 ILQ = .TRUE. 340 IWANTQ = 3 341 ELSE 342 IWANTQ = 0 343 END IF 344 345 IF( LSAME( WANTZ, 'N' ) ) THEN 346 ILZ = .FALSE. 347 IWANTZ = 1 348 ELSE IF( LSAME( WANTZ, 'V' ) ) THEN 349 ILZ = .TRUE. 350 IWANTZ = 2 351 ELSE IF( LSAME( WANTZ, 'I' ) ) THEN 352 ILZ = .TRUE. 353 IWANTZ = 3 354 ELSE 355 IWANTZ = 0 356 END IF 357* 358* Check Argument Values 359* 360 INFO = 0 361 IF( IWANTS.EQ.0 ) THEN 362 INFO = -1 363 ELSE IF( IWANTQ.EQ.0 ) THEN 364 INFO = -2 365 ELSE IF( IWANTZ.EQ.0 ) THEN 366 INFO = -3 367 ELSE IF( N.LT.0 ) THEN 368 INFO = -4 369 ELSE IF( ILO.LT.1 ) THEN 370 INFO = -5 371 ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN 372 INFO = -6 373 ELSE IF( LDA.LT.N ) THEN 374 INFO = -8 375 ELSE IF( LDB.LT.N ) THEN 376 INFO = -10 377 ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN 378 INFO = -15 379 ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN 380 INFO = -17 381 END IF 382 IF( INFO.NE.0 ) THEN 383 CALL XERBLA( 'CLAQZ0', -INFO ) 384 RETURN 385 END IF 386 387* 388* Quick return if possible 389* 390 IF( N.LE.0 ) THEN 391 WORK( 1 ) = REAL( 1 ) 392 RETURN 393 END IF 394 395* 396* Get the parameters 397* 398 JBCMPZ( 1:1 ) = WANTS 399 JBCMPZ( 2:2 ) = WANTQ 400 JBCMPZ( 3:3 ) = WANTZ 401 402 NMIN = ILAENV( 12, 'CLAQZ0', JBCMPZ, N, ILO, IHI, LWORK ) 403 404 NWR = ILAENV( 13, 'CLAQZ0', JBCMPZ, N, ILO, IHI, LWORK ) 405 NWR = MAX( 2, NWR ) 406 NWR = MIN( IHI-ILO+1, ( N-1 ) / 3, NWR ) 407 408 NIBBLE = ILAENV( 14, 'CLAQZ0', JBCMPZ, N, ILO, IHI, LWORK ) 409 410 NSR = ILAENV( 15, 'CLAQZ0', JBCMPZ, N, ILO, IHI, LWORK ) 411 NSR = MIN( NSR, ( N+6 ) / 9, IHI-ILO ) 412 NSR = MAX( 2, NSR-MOD( NSR, 2 ) ) 413 414 RCOST = ILAENV( 17, 'CLAQZ0', JBCMPZ, N, ILO, IHI, LWORK ) 415 ITEMP1 = INT( NSR/SQRT( 1+2*NSR/( REAL( RCOST )/100*N ) ) ) 416 ITEMP1 = ( ( ITEMP1-1 )/4 )*4+4 417 NBR = NSR+ITEMP1 418 419 IF( N .LT. NMIN .OR. REC .GE. 2 ) THEN 420 CALL CHGEQZ( WANTS, WANTQ, WANTZ, N, ILO, IHI, A, LDA, B, LDB, 421 $ ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, 422 $ INFO ) 423 RETURN 424 END IF 425 426* 427* Find out required workspace 428* 429 430* Workspace query to CLAQZ2 431 NW = MAX( NWR, NMIN ) 432 CALL CLAQZ2( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW, A, LDA, B, LDB, 433 $ Q, LDQ, Z, LDZ, N_UNDEFLATED, N_DEFLATED, ALPHA, 434 $ BETA, WORK, NW, WORK, NW, WORK, -1, RWORK, REC, 435 $ AED_INFO ) 436 ITEMP1 = INT( WORK( 1 ) ) 437* Workspace query to CLAQZ3 438 CALL CLAQZ3( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSR, NBR, ALPHA, 439 $ BETA, A, LDA, B, LDB, Q, LDQ, Z, LDZ, WORK, NBR, 440 $ WORK, NBR, WORK, -1, SWEEP_INFO ) 441 ITEMP2 = INT( WORK( 1 ) ) 442 443 LWORKREQ = MAX( ITEMP1+2*NW**2, ITEMP2+2*NBR**2 ) 444 IF ( LWORK .EQ.-1 ) THEN 445 WORK( 1 ) = REAL( LWORKREQ ) 446 RETURN 447 ELSE IF ( LWORK .LT. LWORKREQ ) THEN 448 INFO = -19 449 END IF 450 IF( INFO.NE.0 ) THEN 451 CALL XERBLA( 'CLAQZ0', INFO ) 452 RETURN 453 END IF 454* 455* Initialize Q and Z 456* 457 IF( IWANTQ.EQ.3 ) CALL CLASET( 'FULL', N, N, CZERO, CONE, Q, 458 $ LDQ ) 459 IF( IWANTZ.EQ.3 ) CALL CLASET( 'FULL', N, N, CZERO, CONE, Z, 460 $ LDZ ) 461 462* Get machine constants 463 SAFMIN = SLAMCH( 'SAFE MINIMUM' ) 464 SAFMAX = ONE/SAFMIN 465 CALL SLABAD( SAFMIN, SAFMAX ) 466 ULP = SLAMCH( 'PRECISION' ) 467 SMLNUM = SAFMIN*( REAL( N )/ULP ) 468 469 ISTART = ILO 470 ISTOP = IHI 471 MAXIT = 30*( IHI-ILO+1 ) 472 LD = 0 473 474 DO IITER = 1, MAXIT 475 IF( IITER .GE. MAXIT ) THEN 476 INFO = ISTOP+1 477 GOTO 80 478 END IF 479 IF ( ISTART+1 .GE. ISTOP ) THEN 480 ISTOP = ISTART 481 EXIT 482 END IF 483 484* Check deflations at the end 485 IF ( ABS( A( ISTOP, ISTOP-1 ) ) .LE. MAX( SMLNUM, 486 $ ULP*( ABS( A( ISTOP, ISTOP ) )+ABS( A( ISTOP-1, 487 $ ISTOP-1 ) ) ) ) ) THEN 488 A( ISTOP, ISTOP-1 ) = CZERO 489 ISTOP = ISTOP-1 490 LD = 0 491 ESHIFT = CZERO 492 END IF 493* Check deflations at the start 494 IF ( ABS( A( ISTART+1, ISTART ) ) .LE. MAX( SMLNUM, 495 $ ULP*( ABS( A( ISTART, ISTART ) )+ABS( A( ISTART+1, 496 $ ISTART+1 ) ) ) ) ) THEN 497 A( ISTART+1, ISTART ) = CZERO 498 ISTART = ISTART+1 499 LD = 0 500 ESHIFT = CZERO 501 END IF 502 503 IF ( ISTART+1 .GE. ISTOP ) THEN 504 EXIT 505 END IF 506 507* Check interior deflations 508 ISTART2 = ISTART 509 DO K = ISTOP, ISTART+1, -1 510 IF ( ABS( A( K, K-1 ) ) .LE. MAX( SMLNUM, ULP*( ABS( A( K, 511 $ K ) )+ABS( A( K-1, K-1 ) ) ) ) ) THEN 512 A( K, K-1 ) = CZERO 513 ISTART2 = K 514 EXIT 515 END IF 516 END DO 517 518* Get range to apply rotations to 519 IF ( ILSCHUR ) THEN 520 ISTARTM = 1 521 ISTOPM = N 522 ELSE 523 ISTARTM = ISTART2 524 ISTOPM = ISTOP 525 END IF 526 527* Check infinite eigenvalues, this is done without blocking so might 528* slow down the method when many infinite eigenvalues are present 529 K = ISTOP 530 DO WHILE ( K.GE.ISTART2 ) 531 TEMPR = ZERO 532 IF( K .LT. ISTOP ) THEN 533 TEMPR = TEMPR+ABS( B( K, K+1 ) ) 534 END IF 535 IF( K .GT. ISTART2 ) THEN 536 TEMPR = TEMPR+ABS( B( K-1, K ) ) 537 END IF 538 539 IF( ABS( B( K, K ) ) .LT. MAX( SMLNUM, ULP*TEMPR ) ) THEN 540* A diagonal element of B is negligable, move it 541* to the top and deflate it 542 543 DO K2 = K, ISTART2+1, -1 544 CALL CLARTG( B( K2-1, K2 ), B( K2-1, K2-1 ), C1, S1, 545 $ TEMP ) 546 B( K2-1, K2 ) = TEMP 547 B( K2-1, K2-1 ) = CZERO 548 549 CALL CROT( K2-2-ISTARTM+1, B( ISTARTM, K2 ), 1, 550 $ B( ISTARTM, K2-1 ), 1, C1, S1 ) 551 CALL CROT( MIN( K2+1, ISTOP )-ISTARTM+1, A( ISTARTM, 552 $ K2 ), 1, A( ISTARTM, K2-1 ), 1, C1, S1 ) 553 IF ( ILZ ) THEN 554 CALL CROT( N, Z( 1, K2 ), 1, Z( 1, K2-1 ), 1, C1, 555 $ S1 ) 556 END IF 557 558 IF( K2.LT.ISTOP ) THEN 559 CALL CLARTG( A( K2, K2-1 ), A( K2+1, K2-1 ), C1, 560 $ S1, TEMP ) 561 A( K2, K2-1 ) = TEMP 562 A( K2+1, K2-1 ) = CZERO 563 564 CALL CROT( ISTOPM-K2+1, A( K2, K2 ), LDA, A( K2+1, 565 $ K2 ), LDA, C1, S1 ) 566 CALL CROT( ISTOPM-K2+1, B( K2, K2 ), LDB, B( K2+1, 567 $ K2 ), LDB, C1, S1 ) 568 IF( ILQ ) THEN 569 CALL CROT( N, Q( 1, K2 ), 1, Q( 1, K2+1 ), 1, 570 $ C1, CONJG( S1 ) ) 571 END IF 572 END IF 573 574 END DO 575 576 IF( ISTART2.LT.ISTOP )THEN 577 CALL CLARTG( A( ISTART2, ISTART2 ), A( ISTART2+1, 578 $ ISTART2 ), C1, S1, TEMP ) 579 A( ISTART2, ISTART2 ) = TEMP 580 A( ISTART2+1, ISTART2 ) = CZERO 581 582 CALL CROT( ISTOPM-( ISTART2+1 )+1, A( ISTART2, 583 $ ISTART2+1 ), LDA, A( ISTART2+1, 584 $ ISTART2+1 ), LDA, C1, S1 ) 585 CALL CROT( ISTOPM-( ISTART2+1 )+1, B( ISTART2, 586 $ ISTART2+1 ), LDB, B( ISTART2+1, 587 $ ISTART2+1 ), LDB, C1, S1 ) 588 IF( ILQ ) THEN 589 CALL CROT( N, Q( 1, ISTART2 ), 1, Q( 1, 590 $ ISTART2+1 ), 1, C1, CONJG( S1 ) ) 591 END IF 592 END IF 593 594 ISTART2 = ISTART2+1 595 596 END IF 597 K = K-1 598 END DO 599 600* istart2 now points to the top of the bottom right 601* unreduced Hessenberg block 602 IF ( ISTART2 .GE. ISTOP ) THEN 603 ISTOP = ISTART2-1 604 LD = 0 605 ESHIFT = CZERO 606 CYCLE 607 END IF 608 609 NW = NWR 610 NSHIFTS = NSR 611 NBLOCK = NBR 612 613 IF ( ISTOP-ISTART2+1 .LT. NMIN ) THEN 614* Setting nw to the size of the subblock will make AED deflate 615* all the eigenvalues. This is slightly more efficient than just 616* using CHGEQZ because the off diagonal part gets updated via BLAS. 617 IF ( ISTOP-ISTART+1 .LT. NMIN ) THEN 618 NW = ISTOP-ISTART+1 619 ISTART2 = ISTART 620 ELSE 621 NW = ISTOP-ISTART2+1 622 END IF 623 END IF 624 625* 626* Time for AED 627* 628 CALL CLAQZ2( ILSCHUR, ILQ, ILZ, N, ISTART2, ISTOP, NW, A, LDA, 629 $ B, LDB, Q, LDQ, Z, LDZ, N_UNDEFLATED, N_DEFLATED, 630 $ ALPHA, BETA, WORK, NW, WORK( NW**2+1 ), NW, 631 $ WORK( 2*NW**2+1 ), LWORK-2*NW**2, RWORK, REC, 632 $ AED_INFO ) 633 634 IF ( N_DEFLATED > 0 ) THEN 635 ISTOP = ISTOP-N_DEFLATED 636 LD = 0 637 ESHIFT = CZERO 638 END IF 639 640 IF ( 100*N_DEFLATED > NIBBLE*( N_DEFLATED+N_UNDEFLATED ) .OR. 641 $ ISTOP-ISTART2+1 .LT. NMIN ) THEN 642* AED has uncovered many eigenvalues. Skip a QZ sweep and run 643* AED again. 644 CYCLE 645 END IF 646 647 LD = LD+1 648 649 NS = MIN( NSHIFTS, ISTOP-ISTART2 ) 650 NS = MIN( NS, N_UNDEFLATED ) 651 SHIFTPOS = ISTOP-N_DEFLATED-N_UNDEFLATED+1 652 653 IF ( MOD( LD, 6 ) .EQ. 0 ) THEN 654* 655* Exceptional shift. Chosen for no particularly good reason. 656* 657 IF( ( REAL( MAXIT )*SAFMIN )*ABS( A( ISTOP, 658 $ ISTOP-1 ) ).LT.ABS( A( ISTOP-1, ISTOP-1 ) ) ) THEN 659 ESHIFT = A( ISTOP, ISTOP-1 )/B( ISTOP-1, ISTOP-1 ) 660 ELSE 661 ESHIFT = ESHIFT+CONE/( SAFMIN*REAL( MAXIT ) ) 662 END IF 663 ALPHA( SHIFTPOS ) = CONE 664 BETA( SHIFTPOS ) = ESHIFT 665 NS = 1 666 END IF 667 668* 669* Time for a QZ sweep 670* 671 CALL CLAQZ3( ILSCHUR, ILQ, ILZ, N, ISTART2, ISTOP, NS, NBLOCK, 672 $ ALPHA( SHIFTPOS ), BETA( SHIFTPOS ), A, LDA, B, 673 $ LDB, Q, LDQ, Z, LDZ, WORK, NBLOCK, WORK( NBLOCK** 674 $ 2+1 ), NBLOCK, WORK( 2*NBLOCK**2+1 ), 675 $ LWORK-2*NBLOCK**2, SWEEP_INFO ) 676 677 END DO 678 679* 680* Call CHGEQZ to normalize the eigenvalue blocks and set the eigenvalues 681* If all the eigenvalues have been found, CHGEQZ will not do any iterations 682* and only normalize the blocks. In case of a rare convergence failure, 683* the single shift might perform better. 684* 685 80 CALL CHGEQZ( WANTS, WANTQ, WANTZ, N, ILO, IHI, A, LDA, B, LDB, 686 $ ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, 687 $ NORM_INFO ) 688 689 INFO = NORM_INFO 690 691 END SUBROUTINE 692