1*> \brief <b> ZGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices</b> 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZGGESX + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zggesx.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zggesx.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zggesx.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZGGESX( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA, 22* B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, 23* LDVSR, RCONDE, RCONDV, WORK, LWORK, RWORK, 24* IWORK, LIWORK, BWORK, INFO ) 25* 26* .. Scalar Arguments .. 27* CHARACTER JOBVSL, JOBVSR, SENSE, SORT 28* INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N, 29* $ SDIM 30* .. 31* .. Array Arguments .. 32* LOGICAL BWORK( * ) 33* INTEGER IWORK( * ) 34* DOUBLE PRECISION RCONDE( 2 ), RCONDV( 2 ), RWORK( * ) 35* COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ), 36* $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ), 37* $ WORK( * ) 38* .. 39* .. Function Arguments .. 40* LOGICAL SELCTG 41* EXTERNAL SELCTG 42* .. 43* 44* 45*> \par Purpose: 46* ============= 47*> 48*> \verbatim 49*> 50*> ZGGESX computes for a pair of N-by-N complex nonsymmetric matrices 51*> (A,B), the generalized eigenvalues, the complex Schur form (S,T), 52*> and, optionally, the left and/or right matrices of Schur vectors (VSL 53*> and VSR). This gives the generalized Schur factorization 54*> 55*> (A,B) = ( (VSL) S (VSR)**H, (VSL) T (VSR)**H ) 56*> 57*> where (VSR)**H is the conjugate-transpose of VSR. 58*> 59*> Optionally, it also orders the eigenvalues so that a selected cluster 60*> of eigenvalues appears in the leading diagonal blocks of the upper 61*> triangular matrix S and the upper triangular matrix T; computes 62*> a reciprocal condition number for the average of the selected 63*> eigenvalues (RCONDE); and computes a reciprocal condition number for 64*> the right and left deflating subspaces corresponding to the selected 65*> eigenvalues (RCONDV). The leading columns of VSL and VSR then form 66*> an orthonormal basis for the corresponding left and right eigenspaces 67*> (deflating subspaces). 68*> 69*> A generalized eigenvalue for a pair of matrices (A,B) is a scalar w 70*> or a ratio alpha/beta = w, such that A - w*B is singular. It is 71*> usually represented as the pair (alpha,beta), as there is a 72*> reasonable interpretation for beta=0 or for both being zero. 73*> 74*> A pair of matrices (S,T) is in generalized complex Schur form if T is 75*> upper triangular with non-negative diagonal and S is upper 76*> triangular. 77*> \endverbatim 78* 79* Arguments: 80* ========== 81* 82*> \param[in] JOBVSL 83*> \verbatim 84*> JOBVSL is CHARACTER*1 85*> = 'N': do not compute the left Schur vectors; 86*> = 'V': compute the left Schur vectors. 87*> \endverbatim 88*> 89*> \param[in] JOBVSR 90*> \verbatim 91*> JOBVSR is CHARACTER*1 92*> = 'N': do not compute the right Schur vectors; 93*> = 'V': compute the right Schur vectors. 94*> \endverbatim 95*> 96*> \param[in] SORT 97*> \verbatim 98*> SORT is CHARACTER*1 99*> Specifies whether or not to order the eigenvalues on the 100*> diagonal of the generalized Schur form. 101*> = 'N': Eigenvalues are not ordered; 102*> = 'S': Eigenvalues are ordered (see SELCTG). 103*> \endverbatim 104*> 105*> \param[in] SELCTG 106*> \verbatim 107*> SELCTG is a LOGICAL FUNCTION of two COMPLEX*16 arguments 108*> SELCTG must be declared EXTERNAL in the calling subroutine. 109*> If SORT = 'N', SELCTG is not referenced. 110*> If SORT = 'S', SELCTG is used to select eigenvalues to sort 111*> to the top left of the Schur form. 112*> Note that a selected complex eigenvalue may no longer satisfy 113*> SELCTG(ALPHA(j),BETA(j)) = .TRUE. after ordering, since 114*> ordering may change the value of complex eigenvalues 115*> (especially if the eigenvalue is ill-conditioned), in this 116*> case INFO is set to N+3 see INFO below). 117*> \endverbatim 118*> 119*> \param[in] SENSE 120*> \verbatim 121*> SENSE is CHARACTER*1 122*> Determines which reciprocal condition numbers are computed. 123*> = 'N': None are computed; 124*> = 'E': Computed for average of selected eigenvalues only; 125*> = 'V': Computed for selected deflating subspaces only; 126*> = 'B': Computed for both. 127*> If SENSE = 'E', 'V', or 'B', SORT must equal 'S'. 128*> \endverbatim 129*> 130*> \param[in] N 131*> \verbatim 132*> N is INTEGER 133*> The order of the matrices A, B, VSL, and VSR. N >= 0. 134*> \endverbatim 135*> 136*> \param[in,out] A 137*> \verbatim 138*> A is COMPLEX*16 array, dimension (LDA, N) 139*> On entry, the first of the pair of matrices. 140*> On exit, A has been overwritten by its generalized Schur 141*> form S. 142*> \endverbatim 143*> 144*> \param[in] LDA 145*> \verbatim 146*> LDA is INTEGER 147*> The leading dimension of A. LDA >= max(1,N). 148*> \endverbatim 149*> 150*> \param[in,out] B 151*> \verbatim 152*> B is COMPLEX*16 array, dimension (LDB, N) 153*> On entry, the second of the pair of matrices. 154*> On exit, B has been overwritten by its generalized Schur 155*> form T. 156*> \endverbatim 157*> 158*> \param[in] LDB 159*> \verbatim 160*> LDB is INTEGER 161*> The leading dimension of B. LDB >= max(1,N). 162*> \endverbatim 163*> 164*> \param[out] SDIM 165*> \verbatim 166*> SDIM is INTEGER 167*> If SORT = 'N', SDIM = 0. 168*> If SORT = 'S', SDIM = number of eigenvalues (after sorting) 169*> for which SELCTG is true. 170*> \endverbatim 171*> 172*> \param[out] ALPHA 173*> \verbatim 174*> ALPHA is COMPLEX*16 array, dimension (N) 175*> \endverbatim 176*> 177*> \param[out] BETA 178*> \verbatim 179*> BETA is COMPLEX*16 array, dimension (N) 180*> On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the 181*> generalized eigenvalues. ALPHA(j) and BETA(j),j=1,...,N are 182*> the diagonals of the complex Schur form (S,T). BETA(j) will 183*> be non-negative real. 184*> 185*> Note: the quotients ALPHA(j)/BETA(j) may easily over- or 186*> underflow, and BETA(j) may even be zero. Thus, the user 187*> should avoid naively computing the ratio alpha/beta. 188*> However, ALPHA will be always less than and usually 189*> comparable with norm(A) in magnitude, and BETA always less 190*> than and usually comparable with norm(B). 191*> \endverbatim 192*> 193*> \param[out] VSL 194*> \verbatim 195*> VSL is COMPLEX*16 array, dimension (LDVSL,N) 196*> If JOBVSL = 'V', VSL will contain the left Schur vectors. 197*> Not referenced if JOBVSL = 'N'. 198*> \endverbatim 199*> 200*> \param[in] LDVSL 201*> \verbatim 202*> LDVSL is INTEGER 203*> The leading dimension of the matrix VSL. LDVSL >=1, and 204*> if JOBVSL = 'V', LDVSL >= N. 205*> \endverbatim 206*> 207*> \param[out] VSR 208*> \verbatim 209*> VSR is COMPLEX*16 array, dimension (LDVSR,N) 210*> If JOBVSR = 'V', VSR will contain the right Schur vectors. 211*> Not referenced if JOBVSR = 'N'. 212*> \endverbatim 213*> 214*> \param[in] LDVSR 215*> \verbatim 216*> LDVSR is INTEGER 217*> The leading dimension of the matrix VSR. LDVSR >= 1, and 218*> if JOBVSR = 'V', LDVSR >= N. 219*> \endverbatim 220*> 221*> \param[out] RCONDE 222*> \verbatim 223*> RCONDE is DOUBLE PRECISION array, dimension ( 2 ) 224*> If SENSE = 'E' or 'B', RCONDE(1) and RCONDE(2) contain the 225*> reciprocal condition numbers for the average of the selected 226*> eigenvalues. 227*> Not referenced if SENSE = 'N' or 'V'. 228*> \endverbatim 229*> 230*> \param[out] RCONDV 231*> \verbatim 232*> RCONDV is DOUBLE PRECISION array, dimension ( 2 ) 233*> If SENSE = 'V' or 'B', RCONDV(1) and RCONDV(2) contain the 234*> reciprocal condition number for the selected deflating 235*> subspaces. 236*> Not referenced if SENSE = 'N' or 'E'. 237*> \endverbatim 238*> 239*> \param[out] WORK 240*> \verbatim 241*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) 242*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 243*> \endverbatim 244*> 245*> \param[in] LWORK 246*> \verbatim 247*> LWORK is INTEGER 248*> The dimension of the array WORK. 249*> If N = 0, LWORK >= 1, else if SENSE = 'E', 'V', or 'B', 250*> LWORK >= MAX(1,2*N,2*SDIM*(N-SDIM)), else 251*> LWORK >= MAX(1,2*N). Note that 2*SDIM*(N-SDIM) <= N*N/2. 252*> Note also that an error is only returned if 253*> LWORK < MAX(1,2*N), but if SENSE = 'E' or 'V' or 'B' this may 254*> not be large enough. 255*> 256*> If LWORK = -1, then a workspace query is assumed; the routine 257*> only calculates the bound on the optimal size of the WORK 258*> array and the minimum size of the IWORK array, returns these 259*> values as the first entries of the WORK and IWORK arrays, and 260*> no error message related to LWORK or LIWORK is issued by 261*> XERBLA. 262*> \endverbatim 263*> 264*> \param[out] RWORK 265*> \verbatim 266*> RWORK is DOUBLE PRECISION array, dimension ( 8*N ) 267*> Real workspace. 268*> \endverbatim 269*> 270*> \param[out] IWORK 271*> \verbatim 272*> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) 273*> On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK. 274*> \endverbatim 275*> 276*> \param[in] LIWORK 277*> \verbatim 278*> LIWORK is INTEGER 279*> The dimension of the array IWORK. 280*> If SENSE = 'N' or N = 0, LIWORK >= 1, otherwise 281*> LIWORK >= N+2. 282*> 283*> If LIWORK = -1, then a workspace query is assumed; the 284*> routine only calculates the bound on the optimal size of the 285*> WORK array and the minimum size of the IWORK array, returns 286*> these values as the first entries of the WORK and IWORK 287*> arrays, and no error message related to LWORK or LIWORK is 288*> issued by XERBLA. 289*> \endverbatim 290*> 291*> \param[out] BWORK 292*> \verbatim 293*> BWORK is LOGICAL array, dimension (N) 294*> Not referenced if SORT = 'N'. 295*> \endverbatim 296*> 297*> \param[out] INFO 298*> \verbatim 299*> INFO is INTEGER 300*> = 0: successful exit 301*> < 0: if INFO = -i, the i-th argument had an illegal value. 302*> = 1,...,N: 303*> The QZ iteration failed. (A,B) are not in Schur 304*> form, but ALPHA(j) and BETA(j) should be correct for 305*> j=INFO+1,...,N. 306*> > N: =N+1: other than QZ iteration failed in ZHGEQZ 307*> =N+2: after reordering, roundoff changed values of 308*> some complex eigenvalues so that leading 309*> eigenvalues in the Generalized Schur form no 310*> longer satisfy SELCTG=.TRUE. This could also 311*> be caused due to scaling. 312*> =N+3: reordering failed in ZTGSEN. 313*> \endverbatim 314* 315* Authors: 316* ======== 317* 318*> \author Univ. of Tennessee 319*> \author Univ. of California Berkeley 320*> \author Univ. of Colorado Denver 321*> \author NAG Ltd. 322* 323*> \ingroup complex16GEeigen 324* 325* ===================================================================== 326 SUBROUTINE ZGGESX( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA, 327 $ B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, 328 $ LDVSR, RCONDE, RCONDV, WORK, LWORK, RWORK, 329 $ IWORK, LIWORK, BWORK, INFO ) 330* 331* -- LAPACK driver routine -- 332* -- LAPACK is a software package provided by Univ. of Tennessee, -- 333* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 334* 335* .. Scalar Arguments .. 336 CHARACTER JOBVSL, JOBVSR, SENSE, SORT 337 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N, 338 $ SDIM 339* .. 340* .. Array Arguments .. 341 LOGICAL BWORK( * ) 342 INTEGER IWORK( * ) 343 DOUBLE PRECISION RCONDE( 2 ), RCONDV( 2 ), RWORK( * ) 344 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ), 345 $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ), 346 $ WORK( * ) 347* .. 348* .. Function Arguments .. 349 LOGICAL SELCTG 350 EXTERNAL SELCTG 351* .. 352* 353* ===================================================================== 354* 355* .. Parameters .. 356 DOUBLE PRECISION ZERO, ONE 357 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 358 COMPLEX*16 CZERO, CONE 359 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 360 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 361* .. 362* .. Local Scalars .. 363 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL, 364 $ LQUERY, WANTSB, WANTSE, WANTSN, WANTST, WANTSV 365 INTEGER I, ICOLS, IERR, IHI, IJOB, IJOBVL, IJOBVR, 366 $ ILEFT, ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK, 367 $ LIWMIN, LWRK, MAXWRK, MINWRK 368 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PL, 369 $ PR, SMLNUM 370* .. 371* .. Local Arrays .. 372 DOUBLE PRECISION DIF( 2 ) 373* .. 374* .. External Subroutines .. 375 EXTERNAL DLABAD, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD, 376 $ ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGSEN, ZUNGQR, 377 $ ZUNMQR 378* .. 379* .. External Functions .. 380 LOGICAL LSAME 381 INTEGER ILAENV 382 DOUBLE PRECISION DLAMCH, ZLANGE 383 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE 384* .. 385* .. Intrinsic Functions .. 386 INTRINSIC MAX, SQRT 387* .. 388* .. Executable Statements .. 389* 390* Decode the input arguments 391* 392 IF( LSAME( JOBVSL, 'N' ) ) THEN 393 IJOBVL = 1 394 ILVSL = .FALSE. 395 ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN 396 IJOBVL = 2 397 ILVSL = .TRUE. 398 ELSE 399 IJOBVL = -1 400 ILVSL = .FALSE. 401 END IF 402* 403 IF( LSAME( JOBVSR, 'N' ) ) THEN 404 IJOBVR = 1 405 ILVSR = .FALSE. 406 ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN 407 IJOBVR = 2 408 ILVSR = .TRUE. 409 ELSE 410 IJOBVR = -1 411 ILVSR = .FALSE. 412 END IF 413* 414 WANTST = LSAME( SORT, 'S' ) 415 WANTSN = LSAME( SENSE, 'N' ) 416 WANTSE = LSAME( SENSE, 'E' ) 417 WANTSV = LSAME( SENSE, 'V' ) 418 WANTSB = LSAME( SENSE, 'B' ) 419 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 420 IF( WANTSN ) THEN 421 IJOB = 0 422 ELSE IF( WANTSE ) THEN 423 IJOB = 1 424 ELSE IF( WANTSV ) THEN 425 IJOB = 2 426 ELSE IF( WANTSB ) THEN 427 IJOB = 4 428 END IF 429* 430* Test the input arguments 431* 432 INFO = 0 433 IF( IJOBVL.LE.0 ) THEN 434 INFO = -1 435 ELSE IF( IJOBVR.LE.0 ) THEN 436 INFO = -2 437 ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN 438 INFO = -3 439 ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR. 440 $ ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN 441 INFO = -5 442 ELSE IF( N.LT.0 ) THEN 443 INFO = -6 444 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 445 INFO = -8 446 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 447 INFO = -10 448 ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN 449 INFO = -15 450 ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN 451 INFO = -17 452 END IF 453* 454* Compute workspace 455* (Note: Comments in the code beginning "Workspace:" describe the 456* minimal amount of workspace needed at that point in the code, 457* as well as the preferred amount for good performance. 458* NB refers to the optimal block size for the immediately 459* following subroutine, as returned by ILAENV.) 460* 461 IF( INFO.EQ.0 ) THEN 462 IF( N.GT.0) THEN 463 MINWRK = 2*N 464 MAXWRK = N*(1 + ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) ) 465 MAXWRK = MAX( MAXWRK, N*( 1 + 466 $ ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, -1 ) ) ) 467 IF( ILVSL ) THEN 468 MAXWRK = MAX( MAXWRK, N*( 1 + 469 $ ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, -1 ) ) ) 470 END IF 471 LWRK = MAXWRK 472 IF( IJOB.GE.1 ) 473 $ LWRK = MAX( LWRK, N*N/2 ) 474 ELSE 475 MINWRK = 1 476 MAXWRK = 1 477 LWRK = 1 478 END IF 479 WORK( 1 ) = LWRK 480 IF( WANTSN .OR. N.EQ.0 ) THEN 481 LIWMIN = 1 482 ELSE 483 LIWMIN = N + 2 484 END IF 485 IWORK( 1 ) = LIWMIN 486* 487 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 488 INFO = -21 489 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY) THEN 490 INFO = -24 491 END IF 492 END IF 493* 494 IF( INFO.NE.0 ) THEN 495 CALL XERBLA( 'ZGGESX', -INFO ) 496 RETURN 497 ELSE IF (LQUERY) THEN 498 RETURN 499 END IF 500* 501* Quick return if possible 502* 503 IF( N.EQ.0 ) THEN 504 SDIM = 0 505 RETURN 506 END IF 507* 508* Get machine constants 509* 510 EPS = DLAMCH( 'P' ) 511 SMLNUM = DLAMCH( 'S' ) 512 BIGNUM = ONE / SMLNUM 513 CALL DLABAD( SMLNUM, BIGNUM ) 514 SMLNUM = SQRT( SMLNUM ) / EPS 515 BIGNUM = ONE / SMLNUM 516* 517* Scale A if max element outside range [SMLNUM,BIGNUM] 518* 519 ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK ) 520 ILASCL = .FALSE. 521 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 522 ANRMTO = SMLNUM 523 ILASCL = .TRUE. 524 ELSE IF( ANRM.GT.BIGNUM ) THEN 525 ANRMTO = BIGNUM 526 ILASCL = .TRUE. 527 END IF 528 IF( ILASCL ) 529 $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR ) 530* 531* Scale B if max element outside range [SMLNUM,BIGNUM] 532* 533 BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK ) 534 ILBSCL = .FALSE. 535 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 536 BNRMTO = SMLNUM 537 ILBSCL = .TRUE. 538 ELSE IF( BNRM.GT.BIGNUM ) THEN 539 BNRMTO = BIGNUM 540 ILBSCL = .TRUE. 541 END IF 542 IF( ILBSCL ) 543 $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR ) 544* 545* Permute the matrix to make it more nearly triangular 546* (Real Workspace: need 6*N) 547* 548 ILEFT = 1 549 IRIGHT = N + 1 550 IRWRK = IRIGHT + N 551 CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ), 552 $ RWORK( IRIGHT ), RWORK( IRWRK ), IERR ) 553* 554* Reduce B to triangular form (QR decomposition of B) 555* (Complex Workspace: need N, prefer N*NB) 556* 557 IROWS = IHI + 1 - ILO 558 ICOLS = N + 1 - ILO 559 ITAU = 1 560 IWRK = ITAU + IROWS 561 CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ), 562 $ WORK( IWRK ), LWORK+1-IWRK, IERR ) 563* 564* Apply the unitary transformation to matrix A 565* (Complex Workspace: need N, prefer N*NB) 566* 567 CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB, 568 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ), 569 $ LWORK+1-IWRK, IERR ) 570* 571* Initialize VSL 572* (Complex Workspace: need N, prefer N*NB) 573* 574 IF( ILVSL ) THEN 575 CALL ZLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL ) 576 IF( IROWS.GT.1 ) THEN 577 CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB, 578 $ VSL( ILO+1, ILO ), LDVSL ) 579 END IF 580 CALL ZUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL, 581 $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR ) 582 END IF 583* 584* Initialize VSR 585* 586 IF( ILVSR ) 587 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR ) 588* 589* Reduce to generalized Hessenberg form 590* (Workspace: none needed) 591* 592 CALL ZGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL, 593 $ LDVSL, VSR, LDVSR, IERR ) 594* 595 SDIM = 0 596* 597* Perform QZ algorithm, computing Schur vectors if desired 598* (Complex Workspace: need N) 599* (Real Workspace: need N) 600* 601 IWRK = ITAU 602 CALL ZHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, 603 $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWRK ), 604 $ LWORK+1-IWRK, RWORK( IRWRK ), IERR ) 605 IF( IERR.NE.0 ) THEN 606 IF( IERR.GT.0 .AND. IERR.LE.N ) THEN 607 INFO = IERR 608 ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN 609 INFO = IERR - N 610 ELSE 611 INFO = N + 1 612 END IF 613 GO TO 40 614 END IF 615* 616* Sort eigenvalues ALPHA/BETA and compute the reciprocal of 617* condition number(s) 618* 619 IF( WANTST ) THEN 620* 621* Undo scaling on eigenvalues before SELCTGing 622* 623 IF( ILASCL ) 624 $ CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR ) 625 IF( ILBSCL ) 626 $ CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR ) 627* 628* Select eigenvalues 629* 630 DO 10 I = 1, N 631 BWORK( I ) = SELCTG( ALPHA( I ), BETA( I ) ) 632 10 CONTINUE 633* 634* Reorder eigenvalues, transform Generalized Schur vectors, and 635* compute reciprocal condition numbers 636* (Complex Workspace: If IJOB >= 1, need MAX(1, 2*SDIM*(N-SDIM)) 637* otherwise, need 1 ) 638* 639 CALL ZTGSEN( IJOB, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, 640 $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PL, PR, 641 $ DIF, WORK( IWRK ), LWORK-IWRK+1, IWORK, LIWORK, 642 $ IERR ) 643* 644 IF( IJOB.GE.1 ) 645 $ MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) ) 646 IF( IERR.EQ.-21 ) THEN 647* 648* not enough complex workspace 649* 650 INFO = -21 651 ELSE 652 IF( IJOB.EQ.1 .OR. IJOB.EQ.4 ) THEN 653 RCONDE( 1 ) = PL 654 RCONDE( 2 ) = PR 655 END IF 656 IF( IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN 657 RCONDV( 1 ) = DIF( 1 ) 658 RCONDV( 2 ) = DIF( 2 ) 659 END IF 660 IF( IERR.EQ.1 ) 661 $ INFO = N + 3 662 END IF 663* 664 END IF 665* 666* Apply permutation to VSL and VSR 667* (Workspace: none needed) 668* 669 IF( ILVSL ) 670 $ CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ), 671 $ RWORK( IRIGHT ), N, VSL, LDVSL, IERR ) 672* 673 IF( ILVSR ) 674 $ CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ), 675 $ RWORK( IRIGHT ), N, VSR, LDVSR, IERR ) 676* 677* Undo scaling 678* 679 IF( ILASCL ) THEN 680 CALL ZLASCL( 'U', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR ) 681 CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR ) 682 END IF 683* 684 IF( ILBSCL ) THEN 685 CALL ZLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR ) 686 CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR ) 687 END IF 688* 689 IF( WANTST ) THEN 690* 691* Check if reordering is correct 692* 693 LASTSL = .TRUE. 694 SDIM = 0 695 DO 30 I = 1, N 696 CURSL = SELCTG( ALPHA( I ), BETA( I ) ) 697 IF( CURSL ) 698 $ SDIM = SDIM + 1 699 IF( CURSL .AND. .NOT.LASTSL ) 700 $ INFO = N + 2 701 LASTSL = CURSL 702 30 CONTINUE 703* 704 END IF 705* 706 40 CONTINUE 707* 708 WORK( 1 ) = MAXWRK 709 IWORK( 1 ) = LIWMIN 710* 711 RETURN 712* 713* End of ZGGESX 714* 715 END 716