1*> \brief \b STGSEN 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download STGSEN + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stgsen.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stgsen.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stgsen.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE STGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB, 22* ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, M, PL, 23* PR, DIF, WORK, LWORK, IWORK, LIWORK, INFO ) 24* 25* .. Scalar Arguments .. 26* LOGICAL WANTQ, WANTZ 27* INTEGER IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK, 28* $ M, N 29* REAL PL, PR 30* .. 31* .. Array Arguments .. 32* LOGICAL SELECT( * ) 33* INTEGER IWORK( * ) 34* REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ), 35* $ B( LDB, * ), BETA( * ), DIF( * ), Q( LDQ, * ), 36* $ WORK( * ), Z( LDZ, * ) 37* .. 38* 39* 40*> \par Purpose: 41* ============= 42*> 43*> \verbatim 44*> 45*> STGSEN reorders the generalized real Schur decomposition of a real 46*> matrix pair (A, B) (in terms of an orthonormal equivalence trans- 47*> formation Q**T * (A, B) * Z), so that a selected cluster of eigenvalues 48*> appears in the leading diagonal blocks of the upper quasi-triangular 49*> matrix A and the upper triangular B. The leading columns of Q and 50*> Z form orthonormal bases of the corresponding left and right eigen- 51*> spaces (deflating subspaces). (A, B) must be in generalized real 52*> Schur canonical form (as returned by SGGES), i.e. A is block upper 53*> triangular with 1-by-1 and 2-by-2 diagonal blocks. B is upper 54*> triangular. 55*> 56*> STGSEN also computes the generalized eigenvalues 57*> 58*> w(j) = (ALPHAR(j) + i*ALPHAI(j))/BETA(j) 59*> 60*> of the reordered matrix pair (A, B). 61*> 62*> Optionally, STGSEN computes the estimates of reciprocal condition 63*> numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11), 64*> (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s) 65*> between the matrix pairs (A11, B11) and (A22,B22) that correspond to 66*> the selected cluster and the eigenvalues outside the cluster, resp., 67*> and norms of "projections" onto left and right eigenspaces w.r.t. 68*> the selected cluster in the (1,1)-block. 69*> \endverbatim 70* 71* Arguments: 72* ========== 73* 74*> \param[in] IJOB 75*> \verbatim 76*> IJOB is INTEGER 77*> Specifies whether condition numbers are required for the 78*> cluster of eigenvalues (PL and PR) or the deflating subspaces 79*> (Difu and Difl): 80*> =0: Only reorder w.r.t. SELECT. No extras. 81*> =1: Reciprocal of norms of "projections" onto left and right 82*> eigenspaces w.r.t. the selected cluster (PL and PR). 83*> =2: Upper bounds on Difu and Difl. F-norm-based estimate 84*> (DIF(1:2)). 85*> =3: Estimate of Difu and Difl. 1-norm-based estimate 86*> (DIF(1:2)). 87*> About 5 times as expensive as IJOB = 2. 88*> =4: Compute PL, PR and DIF (i.e. 0, 1 and 2 above): Economic 89*> version to get it all. 90*> =5: Compute PL, PR and DIF (i.e. 0, 1 and 3 above) 91*> \endverbatim 92*> 93*> \param[in] WANTQ 94*> \verbatim 95*> WANTQ is LOGICAL 96*> .TRUE. : update the left transformation matrix Q; 97*> .FALSE.: do not update Q. 98*> \endverbatim 99*> 100*> \param[in] WANTZ 101*> \verbatim 102*> WANTZ is LOGICAL 103*> .TRUE. : update the right transformation matrix Z; 104*> .FALSE.: do not update Z. 105*> \endverbatim 106*> 107*> \param[in] SELECT 108*> \verbatim 109*> SELECT is LOGICAL array, dimension (N) 110*> SELECT specifies the eigenvalues in the selected cluster. 111*> To select a real eigenvalue w(j), SELECT(j) must be set to 112*> .TRUE.. To select a complex conjugate pair of eigenvalues 113*> w(j) and w(j+1), corresponding to a 2-by-2 diagonal block, 114*> either SELECT(j) or SELECT(j+1) or both must be set to 115*> .TRUE.; a complex conjugate pair of eigenvalues must be 116*> either both included in the cluster or both excluded. 117*> \endverbatim 118*> 119*> \param[in] N 120*> \verbatim 121*> N is INTEGER 122*> The order of the matrices A and B. N >= 0. 123*> \endverbatim 124*> 125*> \param[in,out] A 126*> \verbatim 127*> A is REAL array, dimension(LDA,N) 128*> On entry, the upper quasi-triangular matrix A, with (A, B) in 129*> generalized real Schur canonical form. 130*> On exit, A is overwritten by the reordered matrix A. 131*> \endverbatim 132*> 133*> \param[in] LDA 134*> \verbatim 135*> LDA is INTEGER 136*> The leading dimension of the array A. LDA >= max(1,N). 137*> \endverbatim 138*> 139*> \param[in,out] B 140*> \verbatim 141*> B is REAL array, dimension(LDB,N) 142*> On entry, the upper triangular matrix B, with (A, B) in 143*> generalized real Schur canonical form. 144*> On exit, B is overwritten by the reordered matrix B. 145*> \endverbatim 146*> 147*> \param[in] LDB 148*> \verbatim 149*> LDB is INTEGER 150*> The leading dimension of the array B. LDB >= max(1,N). 151*> \endverbatim 152*> 153*> \param[out] ALPHAR 154*> \verbatim 155*> ALPHAR is REAL array, dimension (N) 156*> \endverbatim 157*> 158*> \param[out] ALPHAI 159*> \verbatim 160*> ALPHAI is REAL array, dimension (N) 161*> \endverbatim 162*> 163*> \param[out] BETA 164*> \verbatim 165*> BETA is REAL array, dimension (N) 166*> 167*> On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will 168*> be the generalized eigenvalues. ALPHAR(j) + ALPHAI(j)*i 169*> and BETA(j),j=1,...,N are the diagonals of the complex Schur 170*> form (S,T) that would result if the 2-by-2 diagonal blocks of 171*> the real generalized Schur form of (A,B) were further reduced 172*> to triangular form using complex unitary transformations. 173*> If ALPHAI(j) is zero, then the j-th eigenvalue is real; if 174*> positive, then the j-th and (j+1)-st eigenvalues are a 175*> complex conjugate pair, with ALPHAI(j+1) negative. 176*> \endverbatim 177*> 178*> \param[in,out] Q 179*> \verbatim 180*> Q is REAL array, dimension (LDQ,N) 181*> On entry, if WANTQ = .TRUE., Q is an N-by-N matrix. 182*> On exit, Q has been postmultiplied by the left orthogonal 183*> transformation matrix which reorder (A, B); The leading M 184*> columns of Q form orthonormal bases for the specified pair of 185*> left eigenspaces (deflating subspaces). 186*> If WANTQ = .FALSE., Q is not referenced. 187*> \endverbatim 188*> 189*> \param[in] LDQ 190*> \verbatim 191*> LDQ is INTEGER 192*> The leading dimension of the array Q. LDQ >= 1; 193*> and if WANTQ = .TRUE., LDQ >= N. 194*> \endverbatim 195*> 196*> \param[in,out] Z 197*> \verbatim 198*> Z is REAL array, dimension (LDZ,N) 199*> On entry, if WANTZ = .TRUE., Z is an N-by-N matrix. 200*> On exit, Z has been postmultiplied by the left orthogonal 201*> transformation matrix which reorder (A, B); The leading M 202*> columns of Z form orthonormal bases for the specified pair of 203*> left eigenspaces (deflating subspaces). 204*> If WANTZ = .FALSE., Z is not referenced. 205*> \endverbatim 206*> 207*> \param[in] LDZ 208*> \verbatim 209*> LDZ is INTEGER 210*> The leading dimension of the array Z. LDZ >= 1; 211*> If WANTZ = .TRUE., LDZ >= N. 212*> \endverbatim 213*> 214*> \param[out] M 215*> \verbatim 216*> M is INTEGER 217*> The dimension of the specified pair of left and right eigen- 218*> spaces (deflating subspaces). 0 <= M <= N. 219*> \endverbatim 220*> 221*> \param[out] PL 222*> \verbatim 223*> PL is REAL 224*> \endverbatim 225*> 226*> \param[out] PR 227*> \verbatim 228*> PR is REAL 229*> 230*> If IJOB = 1, 4 or 5, PL, PR are lower bounds on the 231*> reciprocal of the norm of "projections" onto left and right 232*> eigenspaces with respect to the selected cluster. 233*> 0 < PL, PR <= 1. 234*> If M = 0 or M = N, PL = PR = 1. 235*> If IJOB = 0, 2 or 3, PL and PR are not referenced. 236*> \endverbatim 237*> 238*> \param[out] DIF 239*> \verbatim 240*> DIF is REAL array, dimension (2). 241*> If IJOB >= 2, DIF(1:2) store the estimates of Difu and Difl. 242*> If IJOB = 2 or 4, DIF(1:2) are F-norm-based upper bounds on 243*> Difu and Difl. If IJOB = 3 or 5, DIF(1:2) are 1-norm-based 244*> estimates of Difu and Difl. 245*> If M = 0 or N, DIF(1:2) = F-norm([A, B]). 246*> If IJOB = 0 or 1, DIF is not referenced. 247*> \endverbatim 248*> 249*> \param[out] WORK 250*> \verbatim 251*> WORK is REAL array, dimension (MAX(1,LWORK)) 252*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 253*> \endverbatim 254*> 255*> \param[in] LWORK 256*> \verbatim 257*> LWORK is INTEGER 258*> The dimension of the array WORK. LWORK >= 4*N+16. 259*> If IJOB = 1, 2 or 4, LWORK >= MAX(4*N+16, 2*M*(N-M)). 260*> If IJOB = 3 or 5, LWORK >= MAX(4*N+16, 4*M*(N-M)). 261*> 262*> If LWORK = -1, then a workspace query is assumed; the routine 263*> only calculates the optimal size of the WORK array, returns 264*> this value as the first entry of the WORK array, and no error 265*> message related to LWORK is issued by XERBLA. 266*> \endverbatim 267*> 268*> \param[out] IWORK 269*> \verbatim 270*> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) 271*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 272*> \endverbatim 273*> 274*> \param[in] LIWORK 275*> \verbatim 276*> LIWORK is INTEGER 277*> The dimension of the array IWORK. LIWORK >= 1. 278*> If IJOB = 1, 2 or 4, LIWORK >= N+6. 279*> If IJOB = 3 or 5, LIWORK >= MAX(2*M*(N-M), N+6). 280*> 281*> If LIWORK = -1, then a workspace query is assumed; the 282*> routine only calculates the optimal size of the IWORK array, 283*> returns this value as the first entry of the IWORK array, and 284*> no error message related to LIWORK is issued by XERBLA. 285*> \endverbatim 286*> 287*> \param[out] INFO 288*> \verbatim 289*> INFO is INTEGER 290*> =0: Successful exit. 291*> <0: If INFO = -i, the i-th argument had an illegal value. 292*> =1: Reordering of (A, B) failed because the transformed 293*> matrix pair (A, B) would be too far from generalized 294*> Schur form; the problem is very ill-conditioned. 295*> (A, B) may have been partially reordered. 296*> If requested, 0 is returned in DIF(*), PL and PR. 297*> \endverbatim 298* 299* Authors: 300* ======== 301* 302*> \author Univ. of Tennessee 303*> \author Univ. of California Berkeley 304*> \author Univ. of Colorado Denver 305*> \author NAG Ltd. 306* 307*> \ingroup realOTHERcomputational 308* 309*> \par Further Details: 310* ===================== 311*> 312*> \verbatim 313*> 314*> STGSEN first collects the selected eigenvalues by computing 315*> orthogonal U and W that move them to the top left corner of (A, B). 316*> In other words, the selected eigenvalues are the eigenvalues of 317*> (A11, B11) in: 318*> 319*> U**T*(A, B)*W = (A11 A12) (B11 B12) n1 320*> ( 0 A22),( 0 B22) n2 321*> n1 n2 n1 n2 322*> 323*> where N = n1+n2 and U**T means the transpose of U. The first n1 columns 324*> of U and W span the specified pair of left and right eigenspaces 325*> (deflating subspaces) of (A, B). 326*> 327*> If (A, B) has been obtained from the generalized real Schur 328*> decomposition of a matrix pair (C, D) = Q*(A, B)*Z**T, then the 329*> reordered generalized real Schur form of (C, D) is given by 330*> 331*> (C, D) = (Q*U)*(U**T*(A, B)*W)*(Z*W)**T, 332*> 333*> and the first n1 columns of Q*U and Z*W span the corresponding 334*> deflating subspaces of (C, D) (Q and Z store Q*U and Z*W, resp.). 335*> 336*> Note that if the selected eigenvalue is sufficiently ill-conditioned, 337*> then its value may differ significantly from its value before 338*> reordering. 339*> 340*> The reciprocal condition numbers of the left and right eigenspaces 341*> spanned by the first n1 columns of U and W (or Q*U and Z*W) may 342*> be returned in DIF(1:2), corresponding to Difu and Difl, resp. 343*> 344*> The Difu and Difl are defined as: 345*> 346*> Difu[(A11, B11), (A22, B22)] = sigma-min( Zu ) 347*> and 348*> Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11, B11)], 349*> 350*> where sigma-min(Zu) is the smallest singular value of the 351*> (2*n1*n2)-by-(2*n1*n2) matrix 352*> 353*> Zu = [ kron(In2, A11) -kron(A22**T, In1) ] 354*> [ kron(In2, B11) -kron(B22**T, In1) ]. 355*> 356*> Here, Inx is the identity matrix of size nx and A22**T is the 357*> transpose of A22. kron(X, Y) is the Kronecker product between 358*> the matrices X and Y. 359*> 360*> When DIF(2) is small, small changes in (A, B) can cause large changes 361*> in the deflating subspace. An approximate (asymptotic) bound on the 362*> maximum angular error in the computed deflating subspaces is 363*> 364*> EPS * norm((A, B)) / DIF(2), 365*> 366*> where EPS is the machine precision. 367*> 368*> The reciprocal norm of the projectors on the left and right 369*> eigenspaces associated with (A11, B11) may be returned in PL and PR. 370*> They are computed as follows. First we compute L and R so that 371*> P*(A, B)*Q is block diagonal, where 372*> 373*> P = ( I -L ) n1 Q = ( I R ) n1 374*> ( 0 I ) n2 and ( 0 I ) n2 375*> n1 n2 n1 n2 376*> 377*> and (L, R) is the solution to the generalized Sylvester equation 378*> 379*> A11*R - L*A22 = -A12 380*> B11*R - L*B22 = -B12 381*> 382*> Then PL = (F-norm(L)**2+1)**(-1/2) and PR = (F-norm(R)**2+1)**(-1/2). 383*> An approximate (asymptotic) bound on the average absolute error of 384*> the selected eigenvalues is 385*> 386*> EPS * norm((A, B)) / PL. 387*> 388*> There are also global error bounds which valid for perturbations up 389*> to a certain restriction: A lower bound (x) on the smallest 390*> F-norm(E,F) for which an eigenvalue of (A11, B11) may move and 391*> coalesce with an eigenvalue of (A22, B22) under perturbation (E,F), 392*> (i.e. (A + E, B + F), is 393*> 394*> x = min(Difu,Difl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)). 395*> 396*> An approximate bound on x can be computed from DIF(1:2), PL and PR. 397*> 398*> If y = ( F-norm(E,F) / x) <= 1, the angles between the perturbed 399*> (L', R') and unperturbed (L, R) left and right deflating subspaces 400*> associated with the selected cluster in the (1,1)-blocks can be 401*> bounded as 402*> 403*> max-angle(L, L') <= arctan( y * PL / (1 - y * (1 - PL * PL)**(1/2)) 404*> max-angle(R, R') <= arctan( y * PR / (1 - y * (1 - PR * PR)**(1/2)) 405*> 406*> See LAPACK User's Guide section 4.11 or the following references 407*> for more information. 408*> 409*> Note that if the default method for computing the Frobenius-norm- 410*> based estimate DIF is not wanted (see SLATDF), then the parameter 411*> IDIFJB (see below) should be changed from 3 to 4 (routine SLATDF 412*> (IJOB = 2 will be used)). See STGSYL for more details. 413*> \endverbatim 414* 415*> \par Contributors: 416* ================== 417*> 418*> Bo Kagstrom and Peter Poromaa, Department of Computing Science, 419*> Umea University, S-901 87 Umea, Sweden. 420* 421*> \par References: 422* ================ 423*> 424*> \verbatim 425*> 426*> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the 427*> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in 428*> M.S. Moonen et al (eds), Linear Algebra for Large Scale and 429*> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218. 430*> 431*> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified 432*> Eigenvalues of a Regular Matrix Pair (A, B) and Condition 433*> Estimation: Theory, Algorithms and Software, 434*> Report UMINF - 94.04, Department of Computing Science, Umea 435*> University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working 436*> Note 87. To appear in Numerical Algorithms, 1996. 437*> 438*> [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software 439*> for Solving the Generalized Sylvester Equation and Estimating the 440*> Separation between Regular Matrix Pairs, Report UMINF - 93.23, 441*> Department of Computing Science, Umea University, S-901 87 Umea, 442*> Sweden, December 1993, Revised April 1994, Also as LAPACK Working 443*> Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1, 444*> 1996. 445*> \endverbatim 446*> 447* ===================================================================== 448 SUBROUTINE STGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB, 449 $ ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, M, PL, 450 $ PR, DIF, WORK, LWORK, IWORK, LIWORK, INFO ) 451* 452* -- LAPACK computational routine -- 453* -- LAPACK is a software package provided by Univ. of Tennessee, -- 454* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 455* 456* .. Scalar Arguments .. 457 LOGICAL WANTQ, WANTZ 458 INTEGER IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK, 459 $ M, N 460 REAL PL, PR 461* .. 462* .. Array Arguments .. 463 LOGICAL SELECT( * ) 464 INTEGER IWORK( * ) 465 REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ), 466 $ B( LDB, * ), BETA( * ), DIF( * ), Q( LDQ, * ), 467 $ WORK( * ), Z( LDZ, * ) 468* .. 469* 470* ===================================================================== 471* 472* .. Parameters .. 473 INTEGER IDIFJB 474 PARAMETER ( IDIFJB = 3 ) 475 REAL ZERO, ONE 476 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 477* .. 478* .. Local Scalars .. 479 LOGICAL LQUERY, PAIR, SWAP, WANTD, WANTD1, WANTD2, 480 $ WANTP 481 INTEGER I, IERR, IJB, K, KASE, KK, KS, LIWMIN, LWMIN, 482 $ MN2, N1, N2 483 REAL DSCALE, DSUM, EPS, RDSCAL, SMLNUM 484* .. 485* .. Local Arrays .. 486 INTEGER ISAVE( 3 ) 487* .. 488* .. External Subroutines .. 489 EXTERNAL SLACN2, SLACPY, SLAG2, SLASSQ, STGEXC, STGSYL, 490 $ XERBLA 491* .. 492* .. External Functions .. 493 REAL SLAMCH 494 EXTERNAL SLAMCH 495* .. 496* .. Intrinsic Functions .. 497 INTRINSIC MAX, SIGN, SQRT 498* .. 499* .. Executable Statements .. 500* 501* Decode and test the input parameters 502* 503 INFO = 0 504 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 505* 506 IF( IJOB.LT.0 .OR. IJOB.GT.5 ) THEN 507 INFO = -1 508 ELSE IF( N.LT.0 ) THEN 509 INFO = -5 510 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 511 INFO = -7 512 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 513 INFO = -9 514 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 515 INFO = -14 516 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 517 INFO = -16 518 END IF 519* 520 IF( INFO.NE.0 ) THEN 521 CALL XERBLA( 'STGSEN', -INFO ) 522 RETURN 523 END IF 524* 525* Get machine constants 526* 527 EPS = SLAMCH( 'P' ) 528 SMLNUM = SLAMCH( 'S' ) / EPS 529 IERR = 0 530* 531 WANTP = IJOB.EQ.1 .OR. IJOB.GE.4 532 WANTD1 = IJOB.EQ.2 .OR. IJOB.EQ.4 533 WANTD2 = IJOB.EQ.3 .OR. IJOB.EQ.5 534 WANTD = WANTD1 .OR. WANTD2 535* 536* Set M to the dimension of the specified pair of deflating 537* subspaces. 538* 539 M = 0 540 PAIR = .FALSE. 541 IF( .NOT.LQUERY .OR. IJOB.NE.0 ) THEN 542 DO 10 K = 1, N 543 IF( PAIR ) THEN 544 PAIR = .FALSE. 545 ELSE 546 IF( K.LT.N ) THEN 547 IF( A( K+1, K ).EQ.ZERO ) THEN 548 IF( SELECT( K ) ) 549 $ M = M + 1 550 ELSE 551 PAIR = .TRUE. 552 IF( SELECT( K ) .OR. SELECT( K+1 ) ) 553 $ M = M + 2 554 END IF 555 ELSE 556 IF( SELECT( N ) ) 557 $ M = M + 1 558 END IF 559 END IF 560 10 CONTINUE 561 END IF 562* 563 IF( IJOB.EQ.1 .OR. IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN 564 LWMIN = MAX( 1, 4*N+16, 2*M*(N-M) ) 565 LIWMIN = MAX( 1, N+6 ) 566 ELSE IF( IJOB.EQ.3 .OR. IJOB.EQ.5 ) THEN 567 LWMIN = MAX( 1, 4*N+16, 4*M*(N-M) ) 568 LIWMIN = MAX( 1, 2*M*(N-M), N+6 ) 569 ELSE 570 LWMIN = MAX( 1, 4*N+16 ) 571 LIWMIN = 1 572 END IF 573* 574 WORK( 1 ) = LWMIN 575 IWORK( 1 ) = LIWMIN 576* 577 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 578 INFO = -22 579 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 580 INFO = -24 581 END IF 582* 583 IF( INFO.NE.0 ) THEN 584 CALL XERBLA( 'STGSEN', -INFO ) 585 RETURN 586 ELSE IF( LQUERY ) THEN 587 RETURN 588 END IF 589* 590* Quick return if possible. 591* 592 IF( M.EQ.N .OR. M.EQ.0 ) THEN 593 IF( WANTP ) THEN 594 PL = ONE 595 PR = ONE 596 END IF 597 IF( WANTD ) THEN 598 DSCALE = ZERO 599 DSUM = ONE 600 DO 20 I = 1, N 601 CALL SLASSQ( N, A( 1, I ), 1, DSCALE, DSUM ) 602 CALL SLASSQ( N, B( 1, I ), 1, DSCALE, DSUM ) 603 20 CONTINUE 604 DIF( 1 ) = DSCALE*SQRT( DSUM ) 605 DIF( 2 ) = DIF( 1 ) 606 END IF 607 GO TO 60 608 END IF 609* 610* Collect the selected blocks at the top-left corner of (A, B). 611* 612 KS = 0 613 PAIR = .FALSE. 614 DO 30 K = 1, N 615 IF( PAIR ) THEN 616 PAIR = .FALSE. 617 ELSE 618* 619 SWAP = SELECT( K ) 620 IF( K.LT.N ) THEN 621 IF( A( K+1, K ).NE.ZERO ) THEN 622 PAIR = .TRUE. 623 SWAP = SWAP .OR. SELECT( K+1 ) 624 END IF 625 END IF 626* 627 IF( SWAP ) THEN 628 KS = KS + 1 629* 630* Swap the K-th block to position KS. 631* Perform the reordering of diagonal blocks in (A, B) 632* by orthogonal transformation matrices and update 633* Q and Z accordingly (if requested): 634* 635 KK = K 636 IF( K.NE.KS ) 637 $ CALL STGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, 638 $ Z, LDZ, KK, KS, WORK, LWORK, IERR ) 639* 640 IF( IERR.GT.0 ) THEN 641* 642* Swap is rejected: exit. 643* 644 INFO = 1 645 IF( WANTP ) THEN 646 PL = ZERO 647 PR = ZERO 648 END IF 649 IF( WANTD ) THEN 650 DIF( 1 ) = ZERO 651 DIF( 2 ) = ZERO 652 END IF 653 GO TO 60 654 END IF 655* 656 IF( PAIR ) 657 $ KS = KS + 1 658 END IF 659 END IF 660 30 CONTINUE 661 IF( WANTP ) THEN 662* 663* Solve generalized Sylvester equation for R and L 664* and compute PL and PR. 665* 666 N1 = M 667 N2 = N - M 668 I = N1 + 1 669 IJB = 0 670 CALL SLACPY( 'Full', N1, N2, A( 1, I ), LDA, WORK, N1 ) 671 CALL SLACPY( 'Full', N1, N2, B( 1, I ), LDB, WORK( N1*N2+1 ), 672 $ N1 ) 673 CALL STGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK, 674 $ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ), N1, 675 $ DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ), 676 $ LWORK-2*N1*N2, IWORK, IERR ) 677* 678* Estimate the reciprocal of norms of "projections" onto left 679* and right eigenspaces. 680* 681 RDSCAL = ZERO 682 DSUM = ONE 683 CALL SLASSQ( N1*N2, WORK, 1, RDSCAL, DSUM ) 684 PL = RDSCAL*SQRT( DSUM ) 685 IF( PL.EQ.ZERO ) THEN 686 PL = ONE 687 ELSE 688 PL = DSCALE / ( SQRT( DSCALE*DSCALE / PL+PL )*SQRT( PL ) ) 689 END IF 690 RDSCAL = ZERO 691 DSUM = ONE 692 CALL SLASSQ( N1*N2, WORK( N1*N2+1 ), 1, RDSCAL, DSUM ) 693 PR = RDSCAL*SQRT( DSUM ) 694 IF( PR.EQ.ZERO ) THEN 695 PR = ONE 696 ELSE 697 PR = DSCALE / ( SQRT( DSCALE*DSCALE / PR+PR )*SQRT( PR ) ) 698 END IF 699 END IF 700* 701 IF( WANTD ) THEN 702* 703* Compute estimates of Difu and Difl. 704* 705 IF( WANTD1 ) THEN 706 N1 = M 707 N2 = N - M 708 I = N1 + 1 709 IJB = IDIFJB 710* 711* Frobenius norm-based Difu-estimate. 712* 713 CALL STGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK, 714 $ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ), 715 $ N1, DSCALE, DIF( 1 ), WORK( 2*N1*N2+1 ), 716 $ LWORK-2*N1*N2, IWORK, IERR ) 717* 718* Frobenius norm-based Difl-estimate. 719* 720 CALL STGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA, WORK, 721 $ N2, B( I, I ), LDB, B, LDB, WORK( N1*N2+1 ), 722 $ N2, DSCALE, DIF( 2 ), WORK( 2*N1*N2+1 ), 723 $ LWORK-2*N1*N2, IWORK, IERR ) 724 ELSE 725* 726* 727* Compute 1-norm-based estimates of Difu and Difl using 728* reversed communication with SLACN2. In each step a 729* generalized Sylvester equation or a transposed variant 730* is solved. 731* 732 KASE = 0 733 N1 = M 734 N2 = N - M 735 I = N1 + 1 736 IJB = 0 737 MN2 = 2*N1*N2 738* 739* 1-norm-based estimate of Difu. 740* 741 40 CONTINUE 742 CALL SLACN2( MN2, WORK( MN2+1 ), WORK, IWORK, DIF( 1 ), 743 $ KASE, ISAVE ) 744 IF( KASE.NE.0 ) THEN 745 IF( KASE.EQ.1 ) THEN 746* 747* Solve generalized Sylvester equation. 748* 749 CALL STGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, 750 $ WORK, N1, B, LDB, B( I, I ), LDB, 751 $ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ), 752 $ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK, 753 $ IERR ) 754 ELSE 755* 756* Solve the transposed variant. 757* 758 CALL STGSYL( 'T', IJB, N1, N2, A, LDA, A( I, I ), LDA, 759 $ WORK, N1, B, LDB, B( I, I ), LDB, 760 $ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ), 761 $ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK, 762 $ IERR ) 763 END IF 764 GO TO 40 765 END IF 766 DIF( 1 ) = DSCALE / DIF( 1 ) 767* 768* 1-norm-based estimate of Difl. 769* 770 50 CONTINUE 771 CALL SLACN2( MN2, WORK( MN2+1 ), WORK, IWORK, DIF( 2 ), 772 $ KASE, ISAVE ) 773 IF( KASE.NE.0 ) THEN 774 IF( KASE.EQ.1 ) THEN 775* 776* Solve generalized Sylvester equation. 777* 778 CALL STGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA, 779 $ WORK, N2, B( I, I ), LDB, B, LDB, 780 $ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ), 781 $ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK, 782 $ IERR ) 783 ELSE 784* 785* Solve the transposed variant. 786* 787 CALL STGSYL( 'T', IJB, N2, N1, A( I, I ), LDA, A, LDA, 788 $ WORK, N2, B( I, I ), LDB, B, LDB, 789 $ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ), 790 $ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK, 791 $ IERR ) 792 END IF 793 GO TO 50 794 END IF 795 DIF( 2 ) = DSCALE / DIF( 2 ) 796* 797 END IF 798 END IF 799* 800 60 CONTINUE 801* 802* Compute generalized eigenvalues of reordered pair (A, B) and 803* normalize the generalized Schur form. 804* 805 PAIR = .FALSE. 806 DO 70 K = 1, N 807 IF( PAIR ) THEN 808 PAIR = .FALSE. 809 ELSE 810* 811 IF( K.LT.N ) THEN 812 IF( A( K+1, K ).NE.ZERO ) THEN 813 PAIR = .TRUE. 814 END IF 815 END IF 816* 817 IF( PAIR ) THEN 818* 819* Compute the eigenvalue(s) at position K. 820* 821 WORK( 1 ) = A( K, K ) 822 WORK( 2 ) = A( K+1, K ) 823 WORK( 3 ) = A( K, K+1 ) 824 WORK( 4 ) = A( K+1, K+1 ) 825 WORK( 5 ) = B( K, K ) 826 WORK( 6 ) = B( K+1, K ) 827 WORK( 7 ) = B( K, K+1 ) 828 WORK( 8 ) = B( K+1, K+1 ) 829 CALL SLAG2( WORK, 2, WORK( 5 ), 2, SMLNUM*EPS, BETA( K ), 830 $ BETA( K+1 ), ALPHAR( K ), ALPHAR( K+1 ), 831 $ ALPHAI( K ) ) 832 ALPHAI( K+1 ) = -ALPHAI( K ) 833* 834 ELSE 835* 836 IF( SIGN( ONE, B( K, K ) ).LT.ZERO ) THEN 837* 838* If B(K,K) is negative, make it positive 839* 840 DO 80 I = 1, N 841 A( K, I ) = -A( K, I ) 842 B( K, I ) = -B( K, I ) 843 IF( WANTQ ) Q( I, K ) = -Q( I, K ) 844 80 CONTINUE 845 END IF 846* 847 ALPHAR( K ) = A( K, K ) 848 ALPHAI( K ) = ZERO 849 BETA( K ) = B( K, K ) 850* 851 END IF 852 END IF 853 70 CONTINUE 854* 855 WORK( 1 ) = LWMIN 856 IWORK( 1 ) = LIWMIN 857* 858 RETURN 859* 860* End of STGSEN 861* 862 END 863