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