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