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