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 procedure) 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*> \date November 2011 324* 325*> \ingroup complex16GEeigen 326* 327* ===================================================================== 328 SUBROUTINE ZGGESX( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA, 329 $ B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, 330 $ LDVSR, RCONDE, RCONDV, WORK, LWORK, RWORK, 331 $ IWORK, LIWORK, BWORK, INFO ) 332* 333* -- LAPACK driver routine (version 3.4.0) -- 334* -- LAPACK is a software package provided by Univ. of Tennessee, -- 335* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 336* November 2011 337* 338* .. Scalar Arguments .. 339 CHARACTER JOBVSL, JOBVSR, SENSE, SORT 340 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N, 341 $ SDIM 342* .. 343* .. Array Arguments .. 344 LOGICAL BWORK( * ) 345 INTEGER IWORK( * ) 346 DOUBLE PRECISION RCONDE( 2 ), RCONDV( 2 ), RWORK( * ) 347 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ), 348 $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ), 349 $ WORK( * ) 350* .. 351* .. Function Arguments .. 352 LOGICAL SELCTG 353 EXTERNAL SELCTG 354* .. 355* 356* ===================================================================== 357* 358* .. Parameters .. 359 DOUBLE PRECISION ZERO, ONE 360 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 361 COMPLEX*16 CZERO, CONE 362 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 363 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 364* .. 365* .. Local Scalars .. 366 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL, 367 $ LQUERY, WANTSB, WANTSE, WANTSN, WANTST, WANTSV 368 INTEGER I, ICOLS, IERR, IHI, IJOB, IJOBVL, IJOBVR, 369 $ ILEFT, ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK, 370 $ LIWMIN, LWRK, MAXWRK, MINWRK 371 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PL, 372 $ PR, SMLNUM 373* .. 374* .. Local Arrays .. 375 DOUBLE PRECISION DIF( 2 ) 376* .. 377* .. External Subroutines .. 378 EXTERNAL DLABAD, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD, 379 $ ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGSEN, ZUNGQR, 380 $ ZUNMQR 381* .. 382* .. External Functions .. 383 LOGICAL LSAME 384 INTEGER ILAENV 385 DOUBLE PRECISION DLAMCH, ZLANGE 386 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE 387* .. 388* .. Intrinsic Functions .. 389 INTRINSIC MAX, SQRT 390* .. 391* .. Executable Statements .. 392* 393* Decode the input arguments 394* 395 IF( LSAME( JOBVSL, 'N' ) ) THEN 396 IJOBVL = 1 397 ILVSL = .FALSE. 398 ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN 399 IJOBVL = 2 400 ILVSL = .TRUE. 401 ELSE 402 IJOBVL = -1 403 ILVSL = .FALSE. 404 END IF 405* 406 IF( LSAME( JOBVSR, 'N' ) ) THEN 407 IJOBVR = 1 408 ILVSR = .FALSE. 409 ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN 410 IJOBVR = 2 411 ILVSR = .TRUE. 412 ELSE 413 IJOBVR = -1 414 ILVSR = .FALSE. 415 END IF 416* 417 WANTST = LSAME( SORT, 'S' ) 418 WANTSN = LSAME( SENSE, 'N' ) 419 WANTSE = LSAME( SENSE, 'E' ) 420 WANTSV = LSAME( SENSE, 'V' ) 421 WANTSB = LSAME( SENSE, 'B' ) 422 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 423 IF( WANTSN ) THEN 424 IJOB = 0 425 ELSE IF( WANTSE ) THEN 426 IJOB = 1 427 ELSE IF( WANTSV ) THEN 428 IJOB = 2 429 ELSE IF( WANTSB ) THEN 430 IJOB = 4 431 END IF 432* 433* Test the input arguments 434* 435 INFO = 0 436 IF( IJOBVL.LE.0 ) THEN 437 INFO = -1 438 ELSE IF( IJOBVR.LE.0 ) THEN 439 INFO = -2 440 ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN 441 INFO = -3 442 ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR. 443 $ ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN 444 INFO = -5 445 ELSE IF( N.LT.0 ) THEN 446 INFO = -6 447 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 448 INFO = -8 449 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 450 INFO = -10 451 ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN 452 INFO = -15 453 ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN 454 INFO = -17 455 END IF 456* 457* Compute workspace 458* (Note: Comments in the code beginning "Workspace:" describe the 459* minimal amount of workspace needed at that point in the code, 460* as well as the preferred amount for good performance. 461* NB refers to the optimal block size for the immediately 462* following subroutine, as returned by ILAENV.) 463* 464 IF( INFO.EQ.0 ) THEN 465 IF( N.GT.0) THEN 466 MINWRK = 2*N 467 MAXWRK = N*(1 + ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) ) 468 MAXWRK = MAX( MAXWRK, N*( 1 + 469 $ ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, -1 ) ) ) 470 IF( ILVSL ) THEN 471 MAXWRK = MAX( MAXWRK, N*( 1 + 472 $ ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, -1 ) ) ) 473 END IF 474 LWRK = MAXWRK 475 IF( IJOB.GE.1 ) 476 $ LWRK = MAX( LWRK, N*N/2 ) 477 ELSE 478 MINWRK = 1 479 MAXWRK = 1 480 LWRK = 1 481 END IF 482 WORK( 1 ) = LWRK 483 IF( WANTSN .OR. N.EQ.0 ) THEN 484 LIWMIN = 1 485 ELSE 486 LIWMIN = N + 2 487 END IF 488 IWORK( 1 ) = LIWMIN 489* 490 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 491 INFO = -21 492 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY) THEN 493 INFO = -24 494 END IF 495 END IF 496* 497 IF( INFO.NE.0 ) THEN 498 CALL XERBLA( 'ZGGESX', -INFO ) 499 RETURN 500 ELSE IF (LQUERY) THEN 501 RETURN 502 END IF 503* 504* Quick return if possible 505* 506 IF( N.EQ.0 ) THEN 507 SDIM = 0 508 RETURN 509 END IF 510* 511* Get machine constants 512* 513 EPS = DLAMCH( 'P' ) 514 SMLNUM = DLAMCH( 'S' ) 515 BIGNUM = ONE / SMLNUM 516 CALL DLABAD( SMLNUM, BIGNUM ) 517 SMLNUM = SQRT( SMLNUM ) / EPS 518 BIGNUM = ONE / SMLNUM 519* 520* Scale A if max element outside range [SMLNUM,BIGNUM] 521* 522 ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK ) 523 ILASCL = .FALSE. 524 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 525 ANRMTO = SMLNUM 526 ILASCL = .TRUE. 527 ELSE IF( ANRM.GT.BIGNUM ) THEN 528 ANRMTO = BIGNUM 529 ILASCL = .TRUE. 530 END IF 531 IF( ILASCL ) 532 $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR ) 533* 534* Scale B if max element outside range [SMLNUM,BIGNUM] 535* 536 BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK ) 537 ILBSCL = .FALSE. 538 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 539 BNRMTO = SMLNUM 540 ILBSCL = .TRUE. 541 ELSE IF( BNRM.GT.BIGNUM ) THEN 542 BNRMTO = BIGNUM 543 ILBSCL = .TRUE. 544 END IF 545 IF( ILBSCL ) 546 $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR ) 547* 548* Permute the matrix to make it more nearly triangular 549* (Real Workspace: need 6*N) 550* 551 ILEFT = 1 552 IRIGHT = N + 1 553 IRWRK = IRIGHT + N 554 CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ), 555 $ RWORK( IRIGHT ), RWORK( IRWRK ), IERR ) 556* 557* Reduce B to triangular form (QR decomposition of B) 558* (Complex Workspace: need N, prefer N*NB) 559* 560 IROWS = IHI + 1 - ILO 561 ICOLS = N + 1 - ILO 562 ITAU = 1 563 IWRK = ITAU + IROWS 564 CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ), 565 $ WORK( IWRK ), LWORK+1-IWRK, IERR ) 566* 567* Apply the unitary transformation to matrix A 568* (Complex Workspace: need N, prefer N*NB) 569* 570 CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB, 571 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ), 572 $ LWORK+1-IWRK, IERR ) 573* 574* Initialize VSL 575* (Complex Workspace: need N, prefer N*NB) 576* 577 IF( ILVSL ) THEN 578 CALL ZLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL ) 579 IF( IROWS.GT.1 ) THEN 580 CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB, 581 $ VSL( ILO+1, ILO ), LDVSL ) 582 END IF 583 CALL ZUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL, 584 $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR ) 585 END IF 586* 587* Initialize VSR 588* 589 IF( ILVSR ) 590 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR ) 591* 592* Reduce to generalized Hessenberg form 593* (Workspace: none needed) 594* 595 CALL ZGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL, 596 $ LDVSL, VSR, LDVSR, IERR ) 597* 598 SDIM = 0 599* 600* Perform QZ algorithm, computing Schur vectors if desired 601* (Complex Workspace: need N) 602* (Real Workspace: need N) 603* 604 IWRK = ITAU 605 CALL ZHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, 606 $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWRK ), 607 $ LWORK+1-IWRK, RWORK( IRWRK ), IERR ) 608 IF( IERR.NE.0 ) THEN 609 IF( IERR.GT.0 .AND. IERR.LE.N ) THEN 610 INFO = IERR 611 ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN 612 INFO = IERR - N 613 ELSE 614 INFO = N + 1 615 END IF 616 GO TO 40 617 END IF 618* 619* Sort eigenvalues ALPHA/BETA and compute the reciprocal of 620* condition number(s) 621* 622 IF( WANTST ) THEN 623* 624* Undo scaling on eigenvalues before SELCTGing 625* 626 IF( ILASCL ) 627 $ CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR ) 628 IF( ILBSCL ) 629 $ CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR ) 630* 631* Select eigenvalues 632* 633 DO 10 I = 1, N 634 BWORK( I ) = SELCTG( ALPHA( I ), BETA( I ) ) 635 10 CONTINUE 636* 637* Reorder eigenvalues, transform Generalized Schur vectors, and 638* compute reciprocal condition numbers 639* (Complex Workspace: If IJOB >= 1, need MAX(1, 2*SDIM*(N-SDIM)) 640* otherwise, need 1 ) 641* 642 CALL ZTGSEN( IJOB, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, 643 $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PL, PR, 644 $ DIF, WORK( IWRK ), LWORK-IWRK+1, IWORK, LIWORK, 645 $ IERR ) 646* 647 IF( IJOB.GE.1 ) 648 $ MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) ) 649 IF( IERR.EQ.-21 ) THEN 650* 651* not enough complex workspace 652* 653 INFO = -21 654 ELSE 655 IF( IJOB.EQ.1 .OR. IJOB.EQ.4 ) THEN 656 RCONDE( 1 ) = PL 657 RCONDE( 2 ) = PR 658 END IF 659 IF( IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN 660 RCONDV( 1 ) = DIF( 1 ) 661 RCONDV( 2 ) = DIF( 2 ) 662 END IF 663 IF( IERR.EQ.1 ) 664 $ INFO = N + 3 665 END IF 666* 667 END IF 668* 669* Apply permutation to VSL and VSR 670* (Workspace: none needed) 671* 672 IF( ILVSL ) 673 $ CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ), 674 $ RWORK( IRIGHT ), N, VSL, LDVSL, IERR ) 675* 676 IF( ILVSR ) 677 $ CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ), 678 $ RWORK( IRIGHT ), N, VSR, LDVSR, IERR ) 679* 680* Undo scaling 681* 682 IF( ILASCL ) THEN 683 CALL ZLASCL( 'U', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR ) 684 CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR ) 685 END IF 686* 687 IF( ILBSCL ) THEN 688 CALL ZLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR ) 689 CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR ) 690 END IF 691* 692 IF( WANTST ) THEN 693* 694* Check if reordering is correct 695* 696 LASTSL = .TRUE. 697 SDIM = 0 698 DO 30 I = 1, N 699 CURSL = SELCTG( ALPHA( I ), BETA( I ) ) 700 IF( CURSL ) 701 $ SDIM = SDIM + 1 702 IF( CURSL .AND. .NOT.LASTSL ) 703 $ INFO = N + 2 704 LASTSL = CURSL 705 30 CONTINUE 706* 707 END IF 708* 709 40 CONTINUE 710* 711 WORK( 1 ) = MAXWRK 712 IWORK( 1 ) = LIWMIN 713* 714 RETURN 715* 716* End of ZGGESX 717* 718 END 719