1*> \brief <b> DGGESX 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 DGGESX + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggesx.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggesx.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggesx.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DGGESX( 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* DOUBLE PRECISION 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*> DGGESX 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 procedure) LOGICAL FUNCTION of three DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N) 187*> \endverbatim 188*> 189*> \param[out] ALPHAI 190*> \verbatim 191*> ALPHAI is DOUBLE PRECISION array, dimension (N) 192*> \endverbatim 193*> 194*> \param[out] BETA 195*> \verbatim 196*> BETA is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DHGEQZ 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 DTGSEN. 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*> \date November 2011 341* 342*> \ingroup doubleGEeigen 343* 344*> \par Further Details: 345* ===================== 346*> 347*> \verbatim 348*> 349*> An approximate (asymptotic) bound on the average absolute error of 350*> the selected eigenvalues is 351*> 352*> EPS * norm((A, B)) / RCONDE( 1 ). 353*> 354*> An approximate (asymptotic) bound on the maximum angular error in 355*> the computed deflating subspaces is 356*> 357*> EPS * norm((A, B)) / RCONDV( 2 ). 358*> 359*> See LAPACK User's Guide, section 4.11 for more information. 360*> \endverbatim 361*> 362* ===================================================================== 363 SUBROUTINE DGGESX( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA, 364 $ B, LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, 365 $ VSR, LDVSR, RCONDE, RCONDV, WORK, LWORK, IWORK, 366 $ LIWORK, BWORK, INFO ) 367* 368* -- LAPACK driver routine (version 3.4.0) -- 369* -- LAPACK is a software package provided by Univ. of Tennessee, -- 370* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 371* November 2011 372* 373* .. Scalar Arguments .. 374 CHARACTER JOBVSL, JOBVSR, SENSE, SORT 375 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N, 376 $ SDIM 377* .. 378* .. Array Arguments .. 379 LOGICAL BWORK( * ) 380 INTEGER IWORK( * ) 381 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ), 382 $ B( LDB, * ), BETA( * ), RCONDE( 2 ), 383 $ RCONDV( 2 ), VSL( LDVSL, * ), VSR( LDVSR, * ), 384 $ WORK( * ) 385* .. 386* .. Function Arguments .. 387 LOGICAL SELCTG 388 EXTERNAL SELCTG 389* .. 390* 391* ===================================================================== 392* 393* .. Parameters .. 394 DOUBLE PRECISION ZERO, ONE 395 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 396* .. 397* .. Local Scalars .. 398 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL, 399 $ LQUERY, LST2SL, WANTSB, WANTSE, WANTSN, WANTST, 400 $ WANTSV 401 INTEGER I, ICOLS, IERR, IHI, IJOB, IJOBVL, IJOBVR, 402 $ ILEFT, ILO, IP, IRIGHT, IROWS, ITAU, IWRK, 403 $ LIWMIN, LWRK, MAXWRK, MINWRK 404 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PL, 405 $ PR, SAFMAX, SAFMIN, SMLNUM 406* .. 407* .. Local Arrays .. 408 DOUBLE PRECISION DIF( 2 ) 409* .. 410* .. External Subroutines .. 411 EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLABAD, 412 $ DLACPY, DLASCL, DLASET, DORGQR, DORMQR, DTGSEN, 413 $ XERBLA 414* .. 415* .. External Functions .. 416 LOGICAL LSAME 417 INTEGER ILAENV 418 DOUBLE PRECISION DLAMCH, DLANGE 419 EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE 420* .. 421* .. Intrinsic Functions .. 422 INTRINSIC ABS, MAX, SQRT 423* .. 424* .. Executable Statements .. 425* 426* Decode the input arguments 427* 428 IF( LSAME( JOBVSL, 'N' ) ) THEN 429 IJOBVL = 1 430 ILVSL = .FALSE. 431 ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN 432 IJOBVL = 2 433 ILVSL = .TRUE. 434 ELSE 435 IJOBVL = -1 436 ILVSL = .FALSE. 437 END IF 438* 439 IF( LSAME( JOBVSR, 'N' ) ) THEN 440 IJOBVR = 1 441 ILVSR = .FALSE. 442 ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN 443 IJOBVR = 2 444 ILVSR = .TRUE. 445 ELSE 446 IJOBVR = -1 447 ILVSR = .FALSE. 448 END IF 449* 450 WANTST = LSAME( SORT, 'S' ) 451 WANTSN = LSAME( SENSE, 'N' ) 452 WANTSE = LSAME( SENSE, 'E' ) 453 WANTSV = LSAME( SENSE, 'V' ) 454 WANTSB = LSAME( SENSE, 'B' ) 455 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 456 IF( WANTSN ) THEN 457 IJOB = 0 458 ELSE IF( WANTSE ) THEN 459 IJOB = 1 460 ELSE IF( WANTSV ) THEN 461 IJOB = 2 462 ELSE IF( WANTSB ) THEN 463 IJOB = 4 464 END IF 465* 466* Test the input arguments 467* 468 INFO = 0 469 IF( IJOBVL.LE.0 ) THEN 470 INFO = -1 471 ELSE IF( IJOBVR.LE.0 ) THEN 472 INFO = -2 473 ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN 474 INFO = -3 475 ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR. 476 $ ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN 477 INFO = -5 478 ELSE IF( N.LT.0 ) THEN 479 INFO = -6 480 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 481 INFO = -8 482 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 483 INFO = -10 484 ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN 485 INFO = -16 486 ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN 487 INFO = -18 488 END IF 489* 490* Compute workspace 491* (Note: Comments in the code beginning "Workspace:" describe the 492* minimal amount of workspace needed at that point in the code, 493* as well as the preferred amount for good performance. 494* NB refers to the optimal block size for the immediately 495* following subroutine, as returned by ILAENV.) 496* 497 IF( INFO.EQ.0 ) THEN 498 IF( N.GT.0) THEN 499 MINWRK = MAX( 8*N, 6*N + 16 ) 500 MAXWRK = MINWRK - N + 501 $ N*ILAENV( 1, 'DGEQRF', ' ', N, 1, N, 0 ) 502 MAXWRK = MAX( MAXWRK, MINWRK - N + 503 $ N*ILAENV( 1, 'DORMQR', ' ', N, 1, N, -1 ) ) 504 IF( ILVSL ) THEN 505 MAXWRK = MAX( MAXWRK, MINWRK - N + 506 $ N*ILAENV( 1, 'DORGQR', ' ', N, 1, N, -1 ) ) 507 END IF 508 LWRK = MAXWRK 509 IF( IJOB.GE.1 ) 510 $ LWRK = MAX( LWRK, N*N/2 ) 511 ELSE 512 MINWRK = 1 513 MAXWRK = 1 514 LWRK = 1 515 END IF 516 WORK( 1 ) = LWRK 517 IF( WANTSN .OR. N.EQ.0 ) THEN 518 LIWMIN = 1 519 ELSE 520 LIWMIN = N + 6 521 END IF 522 IWORK( 1 ) = LIWMIN 523* 524 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 525 INFO = -22 526 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 527 INFO = -24 528 END IF 529 END IF 530* 531 IF( INFO.NE.0 ) THEN 532 CALL XERBLA( 'DGGESX', -INFO ) 533 RETURN 534 ELSE IF (LQUERY) THEN 535 RETURN 536 END IF 537* 538* Quick return if possible 539* 540 IF( N.EQ.0 ) THEN 541 SDIM = 0 542 RETURN 543 END IF 544* 545* Get machine constants 546* 547 EPS = DLAMCH( 'P' ) 548 SAFMIN = DLAMCH( 'S' ) 549 SAFMAX = ONE / SAFMIN 550 CALL DLABAD( SAFMIN, SAFMAX ) 551 SMLNUM = SQRT( SAFMIN ) / EPS 552 BIGNUM = ONE / SMLNUM 553* 554* Scale A if max element outside range [SMLNUM,BIGNUM] 555* 556 ANRM = DLANGE( 'M', N, N, A, LDA, WORK ) 557 ILASCL = .FALSE. 558 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 559 ANRMTO = SMLNUM 560 ILASCL = .TRUE. 561 ELSE IF( ANRM.GT.BIGNUM ) THEN 562 ANRMTO = BIGNUM 563 ILASCL = .TRUE. 564 END IF 565 IF( ILASCL ) 566 $ CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR ) 567* 568* Scale B if max element outside range [SMLNUM,BIGNUM] 569* 570 BNRM = DLANGE( 'M', N, N, B, LDB, WORK ) 571 ILBSCL = .FALSE. 572 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 573 BNRMTO = SMLNUM 574 ILBSCL = .TRUE. 575 ELSE IF( BNRM.GT.BIGNUM ) THEN 576 BNRMTO = BIGNUM 577 ILBSCL = .TRUE. 578 END IF 579 IF( ILBSCL ) 580 $ CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR ) 581* 582* Permute the matrix to make it more nearly triangular 583* (Workspace: need 6*N + 2*N for permutation parameters) 584* 585 ILEFT = 1 586 IRIGHT = N + 1 587 IWRK = IRIGHT + N 588 CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ), 589 $ WORK( IRIGHT ), WORK( IWRK ), IERR ) 590* 591* Reduce B to triangular form (QR decomposition of B) 592* (Workspace: need N, prefer N*NB) 593* 594 IROWS = IHI + 1 - ILO 595 ICOLS = N + 1 - ILO 596 ITAU = IWRK 597 IWRK = ITAU + IROWS 598 CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ), 599 $ WORK( IWRK ), LWORK+1-IWRK, IERR ) 600* 601* Apply the orthogonal transformation to matrix A 602* (Workspace: need N, prefer N*NB) 603* 604 CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB, 605 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ), 606 $ LWORK+1-IWRK, IERR ) 607* 608* Initialize VSL 609* (Workspace: need N, prefer N*NB) 610* 611 IF( ILVSL ) THEN 612 CALL DLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL ) 613 IF( IROWS.GT.1 ) THEN 614 CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB, 615 $ VSL( ILO+1, ILO ), LDVSL ) 616 END IF 617 CALL DORGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL, 618 $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR ) 619 END IF 620* 621* Initialize VSR 622* 623 IF( ILVSR ) 624 $ CALL DLASET( 'Full', N, N, ZERO, ONE, VSR, LDVSR ) 625* 626* Reduce to generalized Hessenberg form 627* (Workspace: none needed) 628* 629 CALL DGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL, 630 $ LDVSL, VSR, LDVSR, IERR ) 631* 632 SDIM = 0 633* 634* Perform QZ algorithm, computing Schur vectors if desired 635* (Workspace: need N) 636* 637 IWRK = ITAU 638 CALL DHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, 639 $ ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, 640 $ WORK( IWRK ), LWORK+1-IWRK, IERR ) 641 IF( IERR.NE.0 ) THEN 642 IF( IERR.GT.0 .AND. IERR.LE.N ) THEN 643 INFO = IERR 644 ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN 645 INFO = IERR - N 646 ELSE 647 INFO = N + 1 648 END IF 649 GO TO 60 650 END IF 651* 652* Sort eigenvalues ALPHA/BETA and compute the reciprocal of 653* condition number(s) 654* (Workspace: If IJOB >= 1, need MAX( 8*(N+1), 2*SDIM*(N-SDIM) ) 655* otherwise, need 8*(N+1) ) 656* 657 IF( WANTST ) THEN 658* 659* Undo scaling on eigenvalues before SELCTGing 660* 661 IF( ILASCL ) THEN 662 CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, 663 $ IERR ) 664 CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, 665 $ IERR ) 666 END IF 667 IF( ILBSCL ) 668 $ CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR ) 669* 670* Select eigenvalues 671* 672 DO 10 I = 1, N 673 BWORK( I ) = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) ) 674 10 CONTINUE 675* 676* Reorder eigenvalues, transform Generalized Schur vectors, and 677* compute reciprocal condition numbers 678* 679 CALL DTGSEN( IJOB, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, 680 $ ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, 681 $ SDIM, PL, PR, DIF, WORK( IWRK ), LWORK-IWRK+1, 682 $ IWORK, LIWORK, IERR ) 683* 684 IF( IJOB.GE.1 ) 685 $ MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) ) 686 IF( IERR.EQ.-22 ) THEN 687* 688* not enough real workspace 689* 690 INFO = -22 691 ELSE 692 IF( IJOB.EQ.1 .OR. IJOB.EQ.4 ) THEN 693 RCONDE( 1 ) = PL 694 RCONDE( 2 ) = PR 695 END IF 696 IF( IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN 697 RCONDV( 1 ) = DIF( 1 ) 698 RCONDV( 2 ) = DIF( 2 ) 699 END IF 700 IF( IERR.EQ.1 ) 701 $ INFO = N + 3 702 END IF 703* 704 END IF 705* 706* Apply permutation to VSL and VSR 707* (Workspace: none needed) 708* 709 IF( ILVSL ) 710 $ CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ), 711 $ WORK( IRIGHT ), N, VSL, LDVSL, IERR ) 712* 713 IF( ILVSR ) 714 $ CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ), 715 $ WORK( IRIGHT ), N, VSR, LDVSR, IERR ) 716* 717* Check if unscaling would cause over/underflow, if so, rescale 718* (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of 719* B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I) 720* 721 IF( ILASCL ) THEN 722 DO 20 I = 1, N 723 IF( ALPHAI( I ).NE.ZERO ) THEN 724 IF( ( ALPHAR( I ) / SAFMAX ).GT.( ANRMTO / ANRM ) .OR. 725 $ ( SAFMIN / ALPHAR( I ) ).GT.( ANRM / ANRMTO ) ) THEN 726 WORK( 1 ) = ABS( A( I, I ) / ALPHAR( I ) ) 727 BETA( I ) = BETA( I )*WORK( 1 ) 728 ALPHAR( I ) = ALPHAR( I )*WORK( 1 ) 729 ALPHAI( I ) = ALPHAI( I )*WORK( 1 ) 730 ELSE IF( ( ALPHAI( I ) / SAFMAX ).GT. 731 $ ( ANRMTO / ANRM ) .OR. 732 $ ( SAFMIN / ALPHAI( I ) ).GT.( ANRM / ANRMTO ) ) 733 $ THEN 734 WORK( 1 ) = ABS( A( I, I+1 ) / ALPHAI( I ) ) 735 BETA( I ) = BETA( I )*WORK( 1 ) 736 ALPHAR( I ) = ALPHAR( I )*WORK( 1 ) 737 ALPHAI( I ) = ALPHAI( I )*WORK( 1 ) 738 END IF 739 END IF 740 20 CONTINUE 741 END IF 742* 743 IF( ILBSCL ) THEN 744 DO 30 I = 1, N 745 IF( ALPHAI( I ).NE.ZERO ) THEN 746 IF( ( BETA( I ) / SAFMAX ).GT.( BNRMTO / BNRM ) .OR. 747 $ ( SAFMIN / BETA( I ) ).GT.( BNRM / BNRMTO ) ) THEN 748 WORK( 1 ) = ABS( B( I, I ) / BETA( I ) ) 749 BETA( I ) = BETA( I )*WORK( 1 ) 750 ALPHAR( I ) = ALPHAR( I )*WORK( 1 ) 751 ALPHAI( I ) = ALPHAI( I )*WORK( 1 ) 752 END IF 753 END IF 754 30 CONTINUE 755 END IF 756* 757* Undo scaling 758* 759 IF( ILASCL ) THEN 760 CALL DLASCL( 'H', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR ) 761 CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR ) 762 CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR ) 763 END IF 764* 765 IF( ILBSCL ) THEN 766 CALL DLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR ) 767 CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR ) 768 END IF 769* 770 IF( WANTST ) THEN 771* 772* Check if reordering is correct 773* 774 LASTSL = .TRUE. 775 LST2SL = .TRUE. 776 SDIM = 0 777 IP = 0 778 DO 50 I = 1, N 779 CURSL = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) ) 780 IF( ALPHAI( I ).EQ.ZERO ) THEN 781 IF( CURSL ) 782 $ SDIM = SDIM + 1 783 IP = 0 784 IF( CURSL .AND. .NOT.LASTSL ) 785 $ INFO = N + 2 786 ELSE 787 IF( IP.EQ.1 ) THEN 788* 789* Last eigenvalue of conjugate pair 790* 791 CURSL = CURSL .OR. LASTSL 792 LASTSL = CURSL 793 IF( CURSL ) 794 $ SDIM = SDIM + 2 795 IP = -1 796 IF( CURSL .AND. .NOT.LST2SL ) 797 $ INFO = N + 2 798 ELSE 799* 800* First eigenvalue of conjugate pair 801* 802 IP = 1 803 END IF 804 END IF 805 LST2SL = LASTSL 806 LASTSL = CURSL 807 50 CONTINUE 808* 809 END IF 810* 811 60 CONTINUE 812* 813 WORK( 1 ) = MAXWRK 814 IWORK( 1 ) = LIWMIN 815* 816 RETURN 817* 818* End of DGGESX 819* 820 END 821