1*> \brief \b STGSNA 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download STGSNA + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stgsna.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stgsna.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stgsna.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE STGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, 22* LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK, 23* IWORK, INFO ) 24* 25* .. Scalar Arguments .. 26* CHARACTER HOWMNY, JOB 27* INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N 28* .. 29* .. Array Arguments .. 30* LOGICAL SELECT( * ) 31* INTEGER IWORK( * ) 32* REAL A( LDA, * ), B( LDB, * ), DIF( * ), S( * ), 33* $ VL( LDVL, * ), VR( LDVR, * ), WORK( * ) 34* .. 35* 36* 37*> \par Purpose: 38* ============= 39*> 40*> \verbatim 41*> 42*> STGSNA estimates reciprocal condition numbers for specified 43*> eigenvalues and/or eigenvectors of a matrix pair (A, B) in 44*> generalized real Schur canonical form (or of any matrix pair 45*> (Q*A*Z**T, Q*B*Z**T) with orthogonal matrices Q and Z, where 46*> Z**T denotes the transpose of Z. 47*> 48*> (A, B) must be in generalized real Schur form (as returned by SGGES), 49*> i.e. A is block upper triangular with 1-by-1 and 2-by-2 diagonal 50*> blocks. B is upper triangular. 51*> 52*> \endverbatim 53* 54* Arguments: 55* ========== 56* 57*> \param[in] JOB 58*> \verbatim 59*> JOB is CHARACTER*1 60*> Specifies whether condition numbers are required for 61*> eigenvalues (S) or eigenvectors (DIF): 62*> = 'E': for eigenvalues only (S); 63*> = 'V': for eigenvectors only (DIF); 64*> = 'B': for both eigenvalues and eigenvectors (S and DIF). 65*> \endverbatim 66*> 67*> \param[in] HOWMNY 68*> \verbatim 69*> HOWMNY is CHARACTER*1 70*> = 'A': compute condition numbers for all eigenpairs; 71*> = 'S': compute condition numbers for selected eigenpairs 72*> specified by the array SELECT. 73*> \endverbatim 74*> 75*> \param[in] SELECT 76*> \verbatim 77*> SELECT is LOGICAL array, dimension (N) 78*> If HOWMNY = 'S', SELECT specifies the eigenpairs for which 79*> condition numbers are required. To select condition numbers 80*> for the eigenpair corresponding to a real eigenvalue w(j), 81*> SELECT(j) must be set to .TRUE.. To select condition numbers 82*> corresponding to a complex conjugate pair of eigenvalues w(j) 83*> and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be 84*> set to .TRUE.. 85*> If HOWMNY = 'A', SELECT is not referenced. 86*> \endverbatim 87*> 88*> \param[in] N 89*> \verbatim 90*> N is INTEGER 91*> The order of the square matrix pair (A, B). N >= 0. 92*> \endverbatim 93*> 94*> \param[in] A 95*> \verbatim 96*> A is REAL array, dimension (LDA,N) 97*> The upper quasi-triangular matrix A in the pair (A,B). 98*> \endverbatim 99*> 100*> \param[in] LDA 101*> \verbatim 102*> LDA is INTEGER 103*> The leading dimension of the array A. LDA >= max(1,N). 104*> \endverbatim 105*> 106*> \param[in] B 107*> \verbatim 108*> B is REAL array, dimension (LDB,N) 109*> The upper triangular matrix B in the pair (A,B). 110*> \endverbatim 111*> 112*> \param[in] LDB 113*> \verbatim 114*> LDB is INTEGER 115*> The leading dimension of the array B. LDB >= max(1,N). 116*> \endverbatim 117*> 118*> \param[in] VL 119*> \verbatim 120*> VL is REAL array, dimension (LDVL,M) 121*> If JOB = 'E' or 'B', VL must contain left eigenvectors of 122*> (A, B), corresponding to the eigenpairs specified by HOWMNY 123*> and SELECT. The eigenvectors must be stored in consecutive 124*> columns of VL, as returned by STGEVC. 125*> If JOB = 'V', VL is not referenced. 126*> \endverbatim 127*> 128*> \param[in] LDVL 129*> \verbatim 130*> LDVL is INTEGER 131*> The leading dimension of the array VL. LDVL >= 1. 132*> If JOB = 'E' or 'B', LDVL >= N. 133*> \endverbatim 134*> 135*> \param[in] VR 136*> \verbatim 137*> VR is REAL array, dimension (LDVR,M) 138*> If JOB = 'E' or 'B', VR must contain right eigenvectors of 139*> (A, B), corresponding to the eigenpairs specified by HOWMNY 140*> and SELECT. The eigenvectors must be stored in consecutive 141*> columns ov VR, as returned by STGEVC. 142*> If JOB = 'V', VR is not referenced. 143*> \endverbatim 144*> 145*> \param[in] LDVR 146*> \verbatim 147*> LDVR is INTEGER 148*> The leading dimension of the array VR. LDVR >= 1. 149*> If JOB = 'E' or 'B', LDVR >= N. 150*> \endverbatim 151*> 152*> \param[out] S 153*> \verbatim 154*> S is REAL array, dimension (MM) 155*> If JOB = 'E' or 'B', the reciprocal condition numbers of the 156*> selected eigenvalues, stored in consecutive elements of the 157*> array. For a complex conjugate pair of eigenvalues two 158*> consecutive elements of S are set to the same value. Thus 159*> S(j), DIF(j), and the j-th columns of VL and VR all 160*> correspond to the same eigenpair (but not in general the 161*> j-th eigenpair, unless all eigenpairs are selected). 162*> If JOB = 'V', S is not referenced. 163*> \endverbatim 164*> 165*> \param[out] DIF 166*> \verbatim 167*> DIF is REAL array, dimension (MM) 168*> If JOB = 'V' or 'B', the estimated reciprocal condition 169*> numbers of the selected eigenvectors, stored in consecutive 170*> elements of the array. For a complex eigenvector two 171*> consecutive elements of DIF are set to the same value. If 172*> the eigenvalues cannot be reordered to compute DIF(j), DIF(j) 173*> is set to 0; this can only occur when the true value would be 174*> very small anyway. 175*> If JOB = 'E', DIF is not referenced. 176*> \endverbatim 177*> 178*> \param[in] MM 179*> \verbatim 180*> MM is INTEGER 181*> The number of elements in the arrays S and DIF. MM >= M. 182*> \endverbatim 183*> 184*> \param[out] M 185*> \verbatim 186*> M is INTEGER 187*> The number of elements of the arrays S and DIF used to store 188*> the specified condition numbers; for each selected real 189*> eigenvalue one element is used, and for each selected complex 190*> conjugate pair of eigenvalues, two elements are used. 191*> If HOWMNY = 'A', M is set to N. 192*> \endverbatim 193*> 194*> \param[out] WORK 195*> \verbatim 196*> WORK is REAL array, dimension (MAX(1,LWORK)) 197*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 198*> \endverbatim 199*> 200*> \param[in] LWORK 201*> \verbatim 202*> LWORK is INTEGER 203*> The dimension of the array WORK. LWORK >= max(1,N). 204*> If JOB = 'V' or 'B' LWORK >= 2*N*(N+2)+16. 205*> 206*> If LWORK = -1, then a workspace query is assumed; the routine 207*> only calculates the optimal size of the WORK array, returns 208*> this value as the first entry of the WORK array, and no error 209*> message related to LWORK is issued by XERBLA. 210*> \endverbatim 211*> 212*> \param[out] IWORK 213*> \verbatim 214*> IWORK is INTEGER array, dimension (N + 6) 215*> If JOB = 'E', IWORK is not referenced. 216*> \endverbatim 217*> 218*> \param[out] INFO 219*> \verbatim 220*> INFO is INTEGER 221*> =0: Successful exit 222*> <0: If INFO = -i, the i-th argument had an illegal value 223*> \endverbatim 224* 225* Authors: 226* ======== 227* 228*> \author Univ. of Tennessee 229*> \author Univ. of California Berkeley 230*> \author Univ. of Colorado Denver 231*> \author NAG Ltd. 232* 233*> \ingroup realOTHERcomputational 234* 235*> \par Further Details: 236* ===================== 237*> 238*> \verbatim 239*> 240*> The reciprocal of the condition number of a generalized eigenvalue 241*> w = (a, b) is defined as 242*> 243*> S(w) = (|u**TAv|**2 + |u**TBv|**2)**(1/2) / (norm(u)*norm(v)) 244*> 245*> where u and v are the left and right eigenvectors of (A, B) 246*> corresponding to w; |z| denotes the absolute value of the complex 247*> number, and norm(u) denotes the 2-norm of the vector u. 248*> The pair (a, b) corresponds to an eigenvalue w = a/b (= u**TAv/u**TBv) 249*> of the matrix pair (A, B). If both a and b equal zero, then (A B) is 250*> singular and S(I) = -1 is returned. 251*> 252*> An approximate error bound on the chordal distance between the i-th 253*> computed generalized eigenvalue w and the corresponding exact 254*> eigenvalue lambda is 255*> 256*> chord(w, lambda) <= EPS * norm(A, B) / S(I) 257*> 258*> where EPS is the machine precision. 259*> 260*> The reciprocal of the condition number DIF(i) of right eigenvector u 261*> and left eigenvector v corresponding to the generalized eigenvalue w 262*> is defined as follows: 263*> 264*> a) If the i-th eigenvalue w = (a,b) is real 265*> 266*> Suppose U and V are orthogonal transformations such that 267*> 268*> U**T*(A, B)*V = (S, T) = ( a * ) ( b * ) 1 269*> ( 0 S22 ),( 0 T22 ) n-1 270*> 1 n-1 1 n-1 271*> 272*> Then the reciprocal condition number DIF(i) is 273*> 274*> Difl((a, b), (S22, T22)) = sigma-min( Zl ), 275*> 276*> where sigma-min(Zl) denotes the smallest singular value of the 277*> 2(n-1)-by-2(n-1) matrix 278*> 279*> Zl = [ kron(a, In-1) -kron(1, S22) ] 280*> [ kron(b, In-1) -kron(1, T22) ] . 281*> 282*> Here In-1 is the identity matrix of size n-1. kron(X, Y) is the 283*> Kronecker product between the matrices X and Y. 284*> 285*> Note that if the default method for computing DIF(i) is wanted 286*> (see SLATDF), then the parameter DIFDRI (see below) should be 287*> changed from 3 to 4 (routine SLATDF(IJOB = 2 will be used)). 288*> See STGSYL for more details. 289*> 290*> b) If the i-th and (i+1)-th eigenvalues are complex conjugate pair, 291*> 292*> Suppose U and V are orthogonal transformations such that 293*> 294*> U**T*(A, B)*V = (S, T) = ( S11 * ) ( T11 * ) 2 295*> ( 0 S22 ),( 0 T22) n-2 296*> 2 n-2 2 n-2 297*> 298*> and (S11, T11) corresponds to the complex conjugate eigenvalue 299*> pair (w, conjg(w)). There exist unitary matrices U1 and V1 such 300*> that 301*> 302*> U1**T*S11*V1 = ( s11 s12 ) and U1**T*T11*V1 = ( t11 t12 ) 303*> ( 0 s22 ) ( 0 t22 ) 304*> 305*> where the generalized eigenvalues w = s11/t11 and 306*> conjg(w) = s22/t22. 307*> 308*> Then the reciprocal condition number DIF(i) is bounded by 309*> 310*> min( d1, max( 1, |real(s11)/real(s22)| )*d2 ) 311*> 312*> where, d1 = Difl((s11, t11), (s22, t22)) = sigma-min(Z1), where 313*> Z1 is the complex 2-by-2 matrix 314*> 315*> Z1 = [ s11 -s22 ] 316*> [ t11 -t22 ], 317*> 318*> This is done by computing (using real arithmetic) the 319*> roots of the characteristical polynomial det(Z1**T * Z1 - lambda I), 320*> where Z1**T denotes the transpose of Z1 and det(X) denotes 321*> the determinant of X. 322*> 323*> and d2 is an upper bound on Difl((S11, T11), (S22, T22)), i.e. an 324*> upper bound on sigma-min(Z2), where Z2 is (2n-2)-by-(2n-2) 325*> 326*> Z2 = [ kron(S11**T, In-2) -kron(I2, S22) ] 327*> [ kron(T11**T, In-2) -kron(I2, T22) ] 328*> 329*> Note that if the default method for computing DIF is wanted (see 330*> SLATDF), then the parameter DIFDRI (see below) should be changed 331*> from 3 to 4 (routine SLATDF(IJOB = 2 will be used)). See STGSYL 332*> for more details. 333*> 334*> For each eigenvalue/vector specified by SELECT, DIF stores a 335*> Frobenius norm-based estimate of Difl. 336*> 337*> An approximate error bound for the i-th computed eigenvector VL(i) or 338*> VR(i) is given by 339*> 340*> EPS * norm(A, B) / DIF(i). 341*> 342*> See ref. [2-3] for more details and further references. 343*> \endverbatim 344* 345*> \par Contributors: 346* ================== 347*> 348*> Bo Kagstrom and Peter Poromaa, Department of Computing Science, 349*> Umea University, S-901 87 Umea, Sweden. 350* 351*> \par References: 352* ================ 353*> 354*> \verbatim 355*> 356*> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the 357*> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in 358*> M.S. Moonen et al (eds), Linear Algebra for Large Scale and 359*> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218. 360*> 361*> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified 362*> Eigenvalues of a Regular Matrix Pair (A, B) and Condition 363*> Estimation: Theory, Algorithms and Software, 364*> Report UMINF - 94.04, Department of Computing Science, Umea 365*> University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working 366*> Note 87. To appear in Numerical Algorithms, 1996. 367*> 368*> [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software 369*> for Solving the Generalized Sylvester Equation and Estimating the 370*> Separation between Regular Matrix Pairs, Report UMINF - 93.23, 371*> Department of Computing Science, Umea University, S-901 87 Umea, 372*> Sweden, December 1993, Revised April 1994, Also as LAPACK Working 373*> Note 75. To appear in ACM Trans. on Math. Software, Vol 22, 374*> No 1, 1996. 375*> \endverbatim 376*> 377* ===================================================================== 378 SUBROUTINE STGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, 379 $ LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK, 380 $ IWORK, INFO ) 381* 382* -- LAPACK computational routine -- 383* -- LAPACK is a software package provided by Univ. of Tennessee, -- 384* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 385* 386* .. Scalar Arguments .. 387 CHARACTER HOWMNY, JOB 388 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N 389* .. 390* .. Array Arguments .. 391 LOGICAL SELECT( * ) 392 INTEGER IWORK( * ) 393 REAL A( LDA, * ), B( LDB, * ), DIF( * ), S( * ), 394 $ VL( LDVL, * ), VR( LDVR, * ), WORK( * ) 395* .. 396* 397* ===================================================================== 398* 399* .. Parameters .. 400 INTEGER DIFDRI 401 PARAMETER ( DIFDRI = 3 ) 402 REAL ZERO, ONE, TWO, FOUR 403 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0, 404 $ FOUR = 4.0E+0 ) 405* .. 406* .. Local Scalars .. 407 LOGICAL LQUERY, PAIR, SOMCON, WANTBH, WANTDF, WANTS 408 INTEGER I, IERR, IFST, ILST, IZ, K, KS, LWMIN, N1, N2 409 REAL ALPHAI, ALPHAR, ALPRQT, BETA, C1, C2, COND, 410 $ EPS, LNRM, RNRM, ROOT1, ROOT2, SCALE, SMLNUM, 411 $ TMPII, TMPIR, TMPRI, TMPRR, UHAV, UHAVI, UHBV, 412 $ UHBVI 413* .. 414* .. Local Arrays .. 415 REAL DUMMY( 1 ), DUMMY1( 1 ) 416* .. 417* .. External Functions .. 418 LOGICAL LSAME 419 REAL SDOT, SLAMCH, SLAPY2, SNRM2 420 EXTERNAL LSAME, SDOT, SLAMCH, SLAPY2, SNRM2 421* .. 422* .. External Subroutines .. 423 EXTERNAL SGEMV, SLACPY, SLAG2, STGEXC, STGSYL, XERBLA 424* .. 425* .. Intrinsic Functions .. 426 INTRINSIC MAX, MIN, SQRT 427* .. 428* .. Executable Statements .. 429* 430* Decode and test the input parameters 431* 432 WANTBH = LSAME( JOB, 'B' ) 433 WANTS = LSAME( JOB, 'E' ) .OR. WANTBH 434 WANTDF = LSAME( JOB, 'V' ) .OR. WANTBH 435* 436 SOMCON = LSAME( HOWMNY, 'S' ) 437* 438 INFO = 0 439 LQUERY = ( LWORK.EQ.-1 ) 440* 441 IF( .NOT.WANTS .AND. .NOT.WANTDF ) THEN 442 INFO = -1 443 ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN 444 INFO = -2 445 ELSE IF( N.LT.0 ) THEN 446 INFO = -4 447 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 448 INFO = -6 449 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 450 INFO = -8 451 ELSE IF( WANTS .AND. LDVL.LT.N ) THEN 452 INFO = -10 453 ELSE IF( WANTS .AND. LDVR.LT.N ) THEN 454 INFO = -12 455 ELSE 456* 457* Set M to the number of eigenpairs for which condition numbers 458* are required, and test MM. 459* 460 IF( SOMCON ) THEN 461 M = 0 462 PAIR = .FALSE. 463 DO 10 K = 1, N 464 IF( PAIR ) THEN 465 PAIR = .FALSE. 466 ELSE 467 IF( K.LT.N ) THEN 468 IF( A( K+1, K ).EQ.ZERO ) THEN 469 IF( SELECT( K ) ) 470 $ M = M + 1 471 ELSE 472 PAIR = .TRUE. 473 IF( SELECT( K ) .OR. SELECT( K+1 ) ) 474 $ M = M + 2 475 END IF 476 ELSE 477 IF( SELECT( N ) ) 478 $ M = M + 1 479 END IF 480 END IF 481 10 CONTINUE 482 ELSE 483 M = N 484 END IF 485* 486 IF( N.EQ.0 ) THEN 487 LWMIN = 1 488 ELSE IF( LSAME( JOB, 'V' ) .OR. LSAME( JOB, 'B' ) ) THEN 489 LWMIN = 2*N*( N + 2 ) + 16 490 ELSE 491 LWMIN = N 492 END IF 493 WORK( 1 ) = LWMIN 494* 495 IF( MM.LT.M ) THEN 496 INFO = -15 497 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 498 INFO = -18 499 END IF 500 END IF 501* 502 IF( INFO.NE.0 ) THEN 503 CALL XERBLA( 'STGSNA', -INFO ) 504 RETURN 505 ELSE IF( LQUERY ) THEN 506 RETURN 507 END IF 508* 509* Quick return if possible 510* 511 IF( N.EQ.0 ) 512 $ RETURN 513* 514* Get machine constants 515* 516 EPS = SLAMCH( 'P' ) 517 SMLNUM = SLAMCH( 'S' ) / EPS 518 KS = 0 519 PAIR = .FALSE. 520* 521 DO 20 K = 1, N 522* 523* Determine whether A(k,k) begins a 1-by-1 or 2-by-2 block. 524* 525 IF( PAIR ) THEN 526 PAIR = .FALSE. 527 GO TO 20 528 ELSE 529 IF( K.LT.N ) 530 $ PAIR = A( K+1, K ).NE.ZERO 531 END IF 532* 533* Determine whether condition numbers are required for the k-th 534* eigenpair. 535* 536 IF( SOMCON ) THEN 537 IF( PAIR ) THEN 538 IF( .NOT.SELECT( K ) .AND. .NOT.SELECT( K+1 ) ) 539 $ GO TO 20 540 ELSE 541 IF( .NOT.SELECT( K ) ) 542 $ GO TO 20 543 END IF 544 END IF 545* 546 KS = KS + 1 547* 548 IF( WANTS ) THEN 549* 550* Compute the reciprocal condition number of the k-th 551* eigenvalue. 552* 553 IF( PAIR ) THEN 554* 555* Complex eigenvalue pair. 556* 557 RNRM = SLAPY2( SNRM2( N, VR( 1, KS ), 1 ), 558 $ SNRM2( N, VR( 1, KS+1 ), 1 ) ) 559 LNRM = SLAPY2( SNRM2( N, VL( 1, KS ), 1 ), 560 $ SNRM2( N, VL( 1, KS+1 ), 1 ) ) 561 CALL SGEMV( 'N', N, N, ONE, A, LDA, VR( 1, KS ), 1, ZERO, 562 $ WORK, 1 ) 563 TMPRR = SDOT( N, WORK, 1, VL( 1, KS ), 1 ) 564 TMPRI = SDOT( N, WORK, 1, VL( 1, KS+1 ), 1 ) 565 CALL SGEMV( 'N', N, N, ONE, A, LDA, VR( 1, KS+1 ), 1, 566 $ ZERO, WORK, 1 ) 567 TMPII = SDOT( N, WORK, 1, VL( 1, KS+1 ), 1 ) 568 TMPIR = SDOT( N, WORK, 1, VL( 1, KS ), 1 ) 569 UHAV = TMPRR + TMPII 570 UHAVI = TMPIR - TMPRI 571 CALL SGEMV( 'N', N, N, ONE, B, LDB, VR( 1, KS ), 1, ZERO, 572 $ WORK, 1 ) 573 TMPRR = SDOT( N, WORK, 1, VL( 1, KS ), 1 ) 574 TMPRI = SDOT( N, WORK, 1, VL( 1, KS+1 ), 1 ) 575 CALL SGEMV( 'N', N, N, ONE, B, LDB, VR( 1, KS+1 ), 1, 576 $ ZERO, WORK, 1 ) 577 TMPII = SDOT( N, WORK, 1, VL( 1, KS+1 ), 1 ) 578 TMPIR = SDOT( N, WORK, 1, VL( 1, KS ), 1 ) 579 UHBV = TMPRR + TMPII 580 UHBVI = TMPIR - TMPRI 581 UHAV = SLAPY2( UHAV, UHAVI ) 582 UHBV = SLAPY2( UHBV, UHBVI ) 583 COND = SLAPY2( UHAV, UHBV ) 584 S( KS ) = COND / ( RNRM*LNRM ) 585 S( KS+1 ) = S( KS ) 586* 587 ELSE 588* 589* Real eigenvalue. 590* 591 RNRM = SNRM2( N, VR( 1, KS ), 1 ) 592 LNRM = SNRM2( N, VL( 1, KS ), 1 ) 593 CALL SGEMV( 'N', N, N, ONE, A, LDA, VR( 1, KS ), 1, ZERO, 594 $ WORK, 1 ) 595 UHAV = SDOT( N, WORK, 1, VL( 1, KS ), 1 ) 596 CALL SGEMV( 'N', N, N, ONE, B, LDB, VR( 1, KS ), 1, ZERO, 597 $ WORK, 1 ) 598 UHBV = SDOT( N, WORK, 1, VL( 1, KS ), 1 ) 599 COND = SLAPY2( UHAV, UHBV ) 600 IF( COND.EQ.ZERO ) THEN 601 S( KS ) = -ONE 602 ELSE 603 S( KS ) = COND / ( RNRM*LNRM ) 604 END IF 605 END IF 606 END IF 607* 608 IF( WANTDF ) THEN 609 IF( N.EQ.1 ) THEN 610 DIF( KS ) = SLAPY2( A( 1, 1 ), B( 1, 1 ) ) 611 GO TO 20 612 END IF 613* 614* Estimate the reciprocal condition number of the k-th 615* eigenvectors. 616 IF( PAIR ) THEN 617* 618* Copy the 2-by 2 pencil beginning at (A(k,k), B(k, k)). 619* Compute the eigenvalue(s) at position K. 620* 621 WORK( 1 ) = A( K, K ) 622 WORK( 2 ) = A( K+1, K ) 623 WORK( 3 ) = A( K, K+1 ) 624 WORK( 4 ) = A( K+1, K+1 ) 625 WORK( 5 ) = B( K, K ) 626 WORK( 6 ) = B( K+1, K ) 627 WORK( 7 ) = B( K, K+1 ) 628 WORK( 8 ) = B( K+1, K+1 ) 629 CALL SLAG2( WORK, 2, WORK( 5 ), 2, SMLNUM*EPS, BETA, 630 $ DUMMY1( 1 ), ALPHAR, DUMMY( 1 ), ALPHAI ) 631 ALPRQT = ONE 632 C1 = TWO*( ALPHAR*ALPHAR+ALPHAI*ALPHAI+BETA*BETA ) 633 C2 = FOUR*BETA*BETA*ALPHAI*ALPHAI 634 ROOT1 = C1 + SQRT( C1*C1-4.0*C2 ) 635 ROOT2 = C2 / ROOT1 636 ROOT1 = ROOT1 / TWO 637 COND = MIN( SQRT( ROOT1 ), SQRT( ROOT2 ) ) 638 END IF 639* 640* Copy the matrix (A, B) to the array WORK and swap the 641* diagonal block beginning at A(k,k) to the (1,1) position. 642* 643 CALL SLACPY( 'Full', N, N, A, LDA, WORK, N ) 644 CALL SLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N ) 645 IFST = K 646 ILST = 1 647* 648 CALL STGEXC( .FALSE., .FALSE., N, WORK, N, WORK( N*N+1 ), N, 649 $ DUMMY, 1, DUMMY1, 1, IFST, ILST, 650 $ WORK( N*N*2+1 ), LWORK-2*N*N, IERR ) 651* 652 IF( IERR.GT.0 ) THEN 653* 654* Ill-conditioned problem - swap rejected. 655* 656 DIF( KS ) = ZERO 657 ELSE 658* 659* Reordering successful, solve generalized Sylvester 660* equation for R and L, 661* A22 * R - L * A11 = A12 662* B22 * R - L * B11 = B12, 663* and compute estimate of Difl((A11,B11), (A22, B22)). 664* 665 N1 = 1 666 IF( WORK( 2 ).NE.ZERO ) 667 $ N1 = 2 668 N2 = N - N1 669 IF( N2.EQ.0 ) THEN 670 DIF( KS ) = COND 671 ELSE 672 I = N*N + 1 673 IZ = 2*N*N + 1 674 CALL STGSYL( 'N', DIFDRI, N2, N1, WORK( N*N1+N1+1 ), 675 $ N, WORK, N, WORK( N1+1 ), N, 676 $ WORK( N*N1+N1+I ), N, WORK( I ), N, 677 $ WORK( N1+I ), N, SCALE, DIF( KS ), 678 $ WORK( IZ+1 ), LWORK-2*N*N, IWORK, IERR ) 679* 680 IF( PAIR ) 681 $ DIF( KS ) = MIN( MAX( ONE, ALPRQT )*DIF( KS ), 682 $ COND ) 683 END IF 684 END IF 685 IF( PAIR ) 686 $ DIF( KS+1 ) = DIF( KS ) 687 END IF 688 IF( PAIR ) 689 $ KS = KS + 1 690* 691 20 CONTINUE 692 WORK( 1 ) = LWMIN 693 RETURN 694* 695* End of STGSNA 696* 697 END 698