001: SUBROUTINE CTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB, 002: $ ALPHA, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF, 003: $ WORK, LWORK, IWORK, LIWORK, INFO ) 004: * 005: * -- LAPACK routine (version 3.2) -- 006: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 007: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 008: * January 2007 009: * 010: * Modified to call CLACN2 in place of CLACON, 10 Feb 03, SJH. 011: * 012: * .. Scalar Arguments .. 013: LOGICAL WANTQ, WANTZ 014: INTEGER IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK, 015: $ M, N 016: REAL PL, PR 017: * .. 018: * .. Array Arguments .. 019: LOGICAL SELECT( * ) 020: INTEGER IWORK( * ) 021: REAL DIF( * ) 022: COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ), 023: $ BETA( * ), Q( LDQ, * ), WORK( * ), Z( LDZ, * ) 024: * .. 025: * 026: * Purpose 027: * ======= 028: * 029: * CTGSEN reorders the generalized Schur decomposition of a complex 030: * matrix pair (A, B) (in terms of an unitary equivalence trans- 031: * formation Q' * (A, B) * Z), so that a selected cluster of eigenvalues 032: * appears in the leading diagonal blocks of the pair (A,B). The leading 033: * columns of Q and Z form unitary bases of the corresponding left and 034: * right eigenspaces (deflating subspaces). (A, B) must be in 035: * generalized Schur canonical form, that is, A and B are both upper 036: * triangular. 037: * 038: * CTGSEN also computes the generalized eigenvalues 039: * 040: * w(j)= ALPHA(j) / BETA(j) 041: * 042: * of the reordered matrix pair (A, B). 043: * 044: * Optionally, the routine computes estimates of reciprocal condition 045: * numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11), 046: * (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s) 047: * between the matrix pairs (A11, B11) and (A22,B22) that correspond to 048: * the selected cluster and the eigenvalues outside the cluster, resp., 049: * and norms of "projections" onto left and right eigenspaces w.r.t. 050: * the selected cluster in the (1,1)-block. 051: * 052: * 053: * Arguments 054: * ========= 055: * 056: * IJOB (input) integer 057: * Specifies whether condition numbers are required for the 058: * cluster of eigenvalues (PL and PR) or the deflating subspaces 059: * (Difu and Difl): 060: * =0: Only reorder w.r.t. SELECT. No extras. 061: * =1: Reciprocal of norms of "projections" onto left and right 062: * eigenspaces w.r.t. the selected cluster (PL and PR). 063: * =2: Upper bounds on Difu and Difl. F-norm-based estimate 064: * (DIF(1:2)). 065: * =3: Estimate of Difu and Difl. 1-norm-based estimate 066: * (DIF(1:2)). 067: * About 5 times as expensive as IJOB = 2. 068: * =4: Compute PL, PR and DIF (i.e. 0, 1 and 2 above): Economic 069: * version to get it all. 070: * =5: Compute PL, PR and DIF (i.e. 0, 1 and 3 above) 071: * 072: * WANTQ (input) LOGICAL 073: * .TRUE. : update the left transformation matrix Q; 074: * .FALSE.: do not update Q. 075: * 076: * WANTZ (input) LOGICAL 077: * .TRUE. : update the right transformation matrix Z; 078: * .FALSE.: do not update Z. 079: * 080: * SELECT (input) LOGICAL array, dimension (N) 081: * SELECT specifies the eigenvalues in the selected cluster. To 082: * select an eigenvalue w(j), SELECT(j) must be set to 083: * .TRUE.. 084: * 085: * N (input) INTEGER 086: * The order of the matrices A and B. N >= 0. 087: * 088: * A (input/output) COMPLEX array, dimension(LDA,N) 089: * On entry, the upper triangular matrix A, in generalized 090: * Schur canonical form. 091: * On exit, A is overwritten by the reordered matrix A. 092: * 093: * LDA (input) INTEGER 094: * The leading dimension of the array A. LDA >= max(1,N). 095: * 096: * B (input/output) COMPLEX array, dimension(LDB,N) 097: * On entry, the upper triangular matrix B, in generalized 098: * Schur canonical form. 099: * On exit, B is overwritten by the reordered matrix B. 100: * 101: * LDB (input) INTEGER 102: * The leading dimension of the array B. LDB >= max(1,N). 103: * 104: * ALPHA (output) COMPLEX array, dimension (N) 105: * BETA (output) COMPLEX array, dimension (N) 106: * The diagonal elements of A and B, respectively, 107: * when the pair (A,B) has been reduced to generalized Schur 108: * form. ALPHA(i)/BETA(i) i=1,...,N are the generalized 109: * eigenvalues. 110: * 111: * Q (input/output) COMPLEX array, dimension (LDQ,N) 112: * On entry, if WANTQ = .TRUE., Q is an N-by-N matrix. 113: * On exit, Q has been postmultiplied by the left unitary 114: * transformation matrix which reorder (A, B); The leading M 115: * columns of Q form orthonormal bases for the specified pair of 116: * left eigenspaces (deflating subspaces). 117: * If WANTQ = .FALSE., Q is not referenced. 118: * 119: * LDQ (input) INTEGER 120: * The leading dimension of the array Q. LDQ >= 1. 121: * If WANTQ = .TRUE., LDQ >= N. 122: * 123: * Z (input/output) COMPLEX array, dimension (LDZ,N) 124: * On entry, if WANTZ = .TRUE., Z is an N-by-N matrix. 125: * On exit, Z has been postmultiplied by the left unitary 126: * transformation matrix which reorder (A, B); The leading M 127: * columns of Z form orthonormal bases for the specified pair of 128: * left eigenspaces (deflating subspaces). 129: * If WANTZ = .FALSE., Z is not referenced. 130: * 131: * LDZ (input) INTEGER 132: * The leading dimension of the array Z. LDZ >= 1. 133: * If WANTZ = .TRUE., LDZ >= N. 134: * 135: * M (output) INTEGER 136: * The dimension of the specified pair of left and right 137: * eigenspaces, (deflating subspaces) 0 <= M <= N. 138: * 139: * PL (output) REAL 140: * PR (output) REAL 141: * If IJOB = 1, 4 or 5, PL, PR are lower bounds on the 142: * reciprocal of the norm of "projections" onto left and right 143: * eigenspace with respect to the selected cluster. 144: * 0 < PL, PR <= 1. 145: * If M = 0 or M = N, PL = PR = 1. 146: * If IJOB = 0, 2 or 3 PL, PR are not referenced. 147: * 148: * DIF (output) REAL array, dimension (2). 149: * If IJOB >= 2, DIF(1:2) store the estimates of Difu and Difl. 150: * If IJOB = 2 or 4, DIF(1:2) are F-norm-based upper bounds on 151: * Difu and Difl. If IJOB = 3 or 5, DIF(1:2) are 1-norm-based 152: * estimates of Difu and Difl, computed using reversed 153: * communication with CLACN2. 154: * If M = 0 or N, DIF(1:2) = F-norm([A, B]). 155: * If IJOB = 0 or 1, DIF is not referenced. 156: * 157: * WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) 158: * IF IJOB = 0, WORK is not referenced. Otherwise, 159: * on exit, if INFO = 0, WORK(1) returns the optimal LWORK. 160: * 161: * LWORK (input) INTEGER 162: * The dimension of the array WORK. LWORK >= 1 163: * If IJOB = 1, 2 or 4, LWORK >= 2*M*(N-M) 164: * If IJOB = 3 or 5, LWORK >= 4*M*(N-M) 165: * 166: * If LWORK = -1, then a workspace query is assumed; the routine 167: * only calculates the optimal size of the WORK array, returns 168: * this value as the first entry of the WORK array, and no error 169: * message related to LWORK is issued by XERBLA. 170: * 171: * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) 172: * IF IJOB = 0, IWORK is not referenced. Otherwise, 173: * on exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 174: * 175: * LIWORK (input) INTEGER 176: * The dimension of the array IWORK. LIWORK >= 1. 177: * If IJOB = 1, 2 or 4, LIWORK >= N+2; 178: * If IJOB = 3 or 5, LIWORK >= MAX(N+2, 2*M*(N-M)); 179: * 180: * If LIWORK = -1, then a workspace query is assumed; the 181: * routine only calculates the optimal size of the IWORK array, 182: * returns this value as the first entry of the IWORK array, and 183: * no error message related to LIWORK is issued by XERBLA. 184: * 185: * INFO (output) INTEGER 186: * =0: Successful exit. 187: * <0: If INFO = -i, the i-th argument had an illegal value. 188: * =1: Reordering of (A, B) failed because the transformed 189: * matrix pair (A, B) would be too far from generalized 190: * Schur form; the problem is very ill-conditioned. 191: * (A, B) may have been partially reordered. 192: * If requested, 0 is returned in DIF(*), PL and PR. 193: * 194: * 195: * Further Details 196: * =============== 197: * 198: * CTGSEN first collects the selected eigenvalues by computing unitary 199: * U and W that move them to the top left corner of (A, B). In other 200: * words, the selected eigenvalues are the eigenvalues of (A11, B11) in 201: * 202: * U'*(A, B)*W = (A11 A12) (B11 B12) n1 203: * ( 0 A22),( 0 B22) n2 204: * n1 n2 n1 n2 205: * 206: * where N = n1+n2 and U' means the conjugate transpose of U. The first 207: * n1 columns of U and W span the specified pair of left and right 208: * eigenspaces (deflating subspaces) of (A, B). 209: * 210: * If (A, B) has been obtained from the generalized real Schur 211: * decomposition of a matrix pair (C, D) = Q*(A, B)*Z', then the 212: * reordered generalized Schur form of (C, D) is given by 213: * 214: * (C, D) = (Q*U)*(U'*(A, B)*W)*(Z*W)', 215: * 216: * and the first n1 columns of Q*U and Z*W span the corresponding 217: * deflating subspaces of (C, D) (Q and Z store Q*U and Z*W, resp.). 218: * 219: * Note that if the selected eigenvalue is sufficiently ill-conditioned, 220: * then its value may differ significantly from its value before 221: * reordering. 222: * 223: * The reciprocal condition numbers of the left and right eigenspaces 224: * spanned by the first n1 columns of U and W (or Q*U and Z*W) may 225: * be returned in DIF(1:2), corresponding to Difu and Difl, resp. 226: * 227: * The Difu and Difl are defined as: 228: * 229: * Difu[(A11, B11), (A22, B22)] = sigma-min( Zu ) 230: * and 231: * Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11, B11)], 232: * 233: * where sigma-min(Zu) is the smallest singular value of the 234: * (2*n1*n2)-by-(2*n1*n2) matrix 235: * 236: * Zu = [ kron(In2, A11) -kron(A22', In1) ] 237: * [ kron(In2, B11) -kron(B22', In1) ]. 238: * 239: * Here, Inx is the identity matrix of size nx and A22' is the 240: * transpose of A22. kron(X, Y) is the Kronecker product between 241: * the matrices X and Y. 242: * 243: * When DIF(2) is small, small changes in (A, B) can cause large changes 244: * in the deflating subspace. An approximate (asymptotic) bound on the 245: * maximum angular error in the computed deflating subspaces is 246: * 247: * EPS * norm((A, B)) / DIF(2), 248: * 249: * where EPS is the machine precision. 250: * 251: * The reciprocal norm of the projectors on the left and right 252: * eigenspaces associated with (A11, B11) may be returned in PL and PR. 253: * They are computed as follows. First we compute L and R so that 254: * P*(A, B)*Q is block diagonal, where 255: * 256: * P = ( I -L ) n1 Q = ( I R ) n1 257: * ( 0 I ) n2 and ( 0 I ) n2 258: * n1 n2 n1 n2 259: * 260: * and (L, R) is the solution to the generalized Sylvester equation 261: * 262: * A11*R - L*A22 = -A12 263: * B11*R - L*B22 = -B12 264: * 265: * Then PL = (F-norm(L)**2+1)**(-1/2) and PR = (F-norm(R)**2+1)**(-1/2). 266: * An approximate (asymptotic) bound on the average absolute error of 267: * the selected eigenvalues is 268: * 269: * EPS * norm((A, B)) / PL. 270: * 271: * There are also global error bounds which valid for perturbations up 272: * to a certain restriction: A lower bound (x) on the smallest 273: * F-norm(E,F) for which an eigenvalue of (A11, B11) may move and 274: * coalesce with an eigenvalue of (A22, B22) under perturbation (E,F), 275: * (i.e. (A + E, B + F), is 276: * 277: * x = min(Difu,Difl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)). 278: * 279: * An approximate bound on x can be computed from DIF(1:2), PL and PR. 280: * 281: * If y = ( F-norm(E,F) / x) <= 1, the angles between the perturbed 282: * (L', R') and unperturbed (L, R) left and right deflating subspaces 283: * associated with the selected cluster in the (1,1)-blocks can be 284: * bounded as 285: * 286: * max-angle(L, L') <= arctan( y * PL / (1 - y * (1 - PL * PL)**(1/2)) 287: * max-angle(R, R') <= arctan( y * PR / (1 - y * (1 - PR * PR)**(1/2)) 288: * 289: * See LAPACK User's Guide section 4.11 or the following references 290: * for more information. 291: * 292: * Note that if the default method for computing the Frobenius-norm- 293: * based estimate DIF is not wanted (see CLATDF), then the parameter 294: * IDIFJB (see below) should be changed from 3 to 4 (routine CLATDF 295: * (IJOB = 2 will be used)). See CTGSYL for more details. 296: * 297: * Based on contributions by 298: * Bo Kagstrom and Peter Poromaa, Department of Computing Science, 299: * Umea University, S-901 87 Umea, Sweden. 300: * 301: * References 302: * ========== 303: * 304: * [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the 305: * Generalized Real Schur Form of a Regular Matrix Pair (A, B), in 306: * M.S. Moonen et al (eds), Linear Algebra for Large Scale and 307: * Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218. 308: * 309: * [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified 310: * Eigenvalues of a Regular Matrix Pair (A, B) and Condition 311: * Estimation: Theory, Algorithms and Software, Report 312: * UMINF - 94.04, Department of Computing Science, Umea University, 313: * S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87. 314: * To appear in Numerical Algorithms, 1996. 315: * 316: * [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software 317: * for Solving the Generalized Sylvester Equation and Estimating the 318: * Separation between Regular Matrix Pairs, Report UMINF - 93.23, 319: * Department of Computing Science, Umea University, S-901 87 Umea, 320: * Sweden, December 1993, Revised April 1994, Also as LAPACK working 321: * Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1, 322: * 1996. 323: * 324: * ===================================================================== 325: * 326: * .. Parameters .. 327: INTEGER IDIFJB 328: PARAMETER ( IDIFJB = 3 ) 329: REAL ZERO, ONE 330: PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 331: * .. 332: * .. Local Scalars .. 333: LOGICAL LQUERY, SWAP, WANTD, WANTD1, WANTD2, WANTP 334: INTEGER I, IERR, IJB, K, KASE, KS, LIWMIN, LWMIN, MN2, 335: $ N1, N2 336: REAL DSCALE, DSUM, RDSCAL, SAFMIN 337: COMPLEX TEMP1, TEMP2 338: * .. 339: * .. Local Arrays .. 340: INTEGER ISAVE( 3 ) 341: * .. 342: * .. External Subroutines .. 343: REAL SLAMCH 344: EXTERNAL CLACN2, CLACPY, CLASSQ, CSCAL, CTGEXC, CTGSYL, 345: $ SLAMCH, XERBLA 346: * .. 347: * .. Intrinsic Functions .. 348: INTRINSIC ABS, CMPLX, CONJG, MAX, SQRT 349: * .. 350: * .. Executable Statements .. 351: * 352: * Decode and test the input parameters 353: * 354: INFO = 0 355: LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 356: * 357: IF( IJOB.LT.0 .OR. IJOB.GT.5 ) THEN 358: INFO = -1 359: ELSE IF( N.LT.0 ) THEN 360: INFO = -5 361: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 362: INFO = -7 363: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 364: INFO = -9 365: ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 366: INFO = -13 367: ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 368: INFO = -15 369: END IF 370: * 371: IF( INFO.NE.0 ) THEN 372: CALL XERBLA( 'CTGSEN', -INFO ) 373: RETURN 374: END IF 375: * 376: IERR = 0 377: * 378: WANTP = IJOB.EQ.1 .OR. IJOB.GE.4 379: WANTD1 = IJOB.EQ.2 .OR. IJOB.EQ.4 380: WANTD2 = IJOB.EQ.3 .OR. IJOB.EQ.5 381: WANTD = WANTD1 .OR. WANTD2 382: * 383: * Set M to the dimension of the specified pair of deflating 384: * subspaces. 385: * 386: M = 0 387: DO 10 K = 1, N 388: ALPHA( K ) = A( K, K ) 389: BETA( K ) = B( K, K ) 390: IF( K.LT.N ) THEN 391: IF( SELECT( K ) ) 392: $ M = M + 1 393: ELSE 394: IF( SELECT( N ) ) 395: $ M = M + 1 396: END IF 397: 10 CONTINUE 398: * 399: IF( IJOB.EQ.1 .OR. IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN 400: LWMIN = MAX( 1, 2*M*(N-M) ) 401: LIWMIN = MAX( 1, N+2 ) 402: ELSE IF( IJOB.EQ.3 .OR. IJOB.EQ.5 ) THEN 403: LWMIN = MAX( 1, 4*M*(N-M) ) 404: LIWMIN = MAX( 1, 2*M*(N-M), N+2 ) 405: ELSE 406: LWMIN = 1 407: LIWMIN = 1 408: END IF 409: * 410: WORK( 1 ) = LWMIN 411: IWORK( 1 ) = LIWMIN 412: * 413: IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 414: INFO = -21 415: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 416: INFO = -23 417: END IF 418: * 419: IF( INFO.NE.0 ) THEN 420: CALL XERBLA( 'CTGSEN', -INFO ) 421: RETURN 422: ELSE IF( LQUERY ) THEN 423: RETURN 424: END IF 425: * 426: * Quick return if possible. 427: * 428: IF( M.EQ.N .OR. M.EQ.0 ) THEN 429: IF( WANTP ) THEN 430: PL = ONE 431: PR = ONE 432: END IF 433: IF( WANTD ) THEN 434: DSCALE = ZERO 435: DSUM = ONE 436: DO 20 I = 1, N 437: CALL CLASSQ( N, A( 1, I ), 1, DSCALE, DSUM ) 438: CALL CLASSQ( N, B( 1, I ), 1, DSCALE, DSUM ) 439: 20 CONTINUE 440: DIF( 1 ) = DSCALE*SQRT( DSUM ) 441: DIF( 2 ) = DIF( 1 ) 442: END IF 443: GO TO 70 444: END IF 445: * 446: * Get machine constant 447: * 448: SAFMIN = SLAMCH( 'S' ) 449: * 450: * Collect the selected blocks at the top-left corner of (A, B). 451: * 452: KS = 0 453: DO 30 K = 1, N 454: SWAP = SELECT( K ) 455: IF( SWAP ) THEN 456: KS = KS + 1 457: * 458: * Swap the K-th block to position KS. Compute unitary Q 459: * and Z that will swap adjacent diagonal blocks in (A, B). 460: * 461: IF( K.NE.KS ) 462: $ CALL CTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 463: $ LDZ, K, KS, IERR ) 464: * 465: IF( IERR.GT.0 ) THEN 466: * 467: * Swap is rejected: exit. 468: * 469: INFO = 1 470: IF( WANTP ) THEN 471: PL = ZERO 472: PR = ZERO 473: END IF 474: IF( WANTD ) THEN 475: DIF( 1 ) = ZERO 476: DIF( 2 ) = ZERO 477: END IF 478: GO TO 70 479: END IF 480: END IF 481: 30 CONTINUE 482: IF( WANTP ) THEN 483: * 484: * Solve generalized Sylvester equation for R and L: 485: * A11 * R - L * A22 = A12 486: * B11 * R - L * B22 = B12 487: * 488: N1 = M 489: N2 = N - M 490: I = N1 + 1 491: CALL CLACPY( 'Full', N1, N2, A( 1, I ), LDA, WORK, N1 ) 492: CALL CLACPY( 'Full', N1, N2, B( 1, I ), LDB, WORK( N1*N2+1 ), 493: $ N1 ) 494: IJB = 0 495: CALL CTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK, 496: $ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ), N1, 497: $ DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ), 498: $ LWORK-2*N1*N2, IWORK, IERR ) 499: * 500: * Estimate the reciprocal of norms of "projections" onto 501: * left and right eigenspaces 502: * 503: RDSCAL = ZERO 504: DSUM = ONE 505: CALL CLASSQ( N1*N2, WORK, 1, RDSCAL, DSUM ) 506: PL = RDSCAL*SQRT( DSUM ) 507: IF( PL.EQ.ZERO ) THEN 508: PL = ONE 509: ELSE 510: PL = DSCALE / ( SQRT( DSCALE*DSCALE / PL+PL )*SQRT( PL ) ) 511: END IF 512: RDSCAL = ZERO 513: DSUM = ONE 514: CALL CLASSQ( N1*N2, WORK( N1*N2+1 ), 1, RDSCAL, DSUM ) 515: PR = RDSCAL*SQRT( DSUM ) 516: IF( PR.EQ.ZERO ) THEN 517: PR = ONE 518: ELSE 519: PR = DSCALE / ( SQRT( DSCALE*DSCALE / PR+PR )*SQRT( PR ) ) 520: END IF 521: END IF 522: IF( WANTD ) THEN 523: * 524: * Compute estimates Difu and Difl. 525: * 526: IF( WANTD1 ) THEN 527: N1 = M 528: N2 = N - M 529: I = N1 + 1 530: IJB = IDIFJB 531: * 532: * Frobenius norm-based Difu estimate. 533: * 534: CALL CTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK, 535: $ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ), 536: $ N1, DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ), 537: $ LWORK-2*N1*N2, IWORK, IERR ) 538: * 539: * Frobenius norm-based Difl estimate. 540: * 541: CALL CTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA, WORK, 542: $ N2, B( I, I ), LDB, B, LDB, WORK( N1*N2+1 ), 543: $ N2, DSCALE, DIF( 2 ), WORK( N1*N2*2+1 ), 544: $ LWORK-2*N1*N2, IWORK, IERR ) 545: ELSE 546: * 547: * Compute 1-norm-based estimates of Difu and Difl using 548: * reversed communication with CLACN2. In each step a 549: * generalized Sylvester equation or a transposed variant 550: * is solved. 551: * 552: KASE = 0 553: N1 = M 554: N2 = N - M 555: I = N1 + 1 556: IJB = 0 557: MN2 = 2*N1*N2 558: * 559: * 1-norm-based estimate of Difu. 560: * 561: 40 CONTINUE 562: CALL CLACN2( MN2, WORK( MN2+1 ), WORK, DIF( 1 ), KASE, 563: $ ISAVE ) 564: IF( KASE.NE.0 ) THEN 565: IF( KASE.EQ.1 ) THEN 566: * 567: * Solve generalized Sylvester equation 568: * 569: CALL CTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, 570: $ WORK, N1, B, LDB, B( I, I ), LDB, 571: $ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ), 572: $ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK, 573: $ IERR ) 574: ELSE 575: * 576: * Solve the transposed variant. 577: * 578: CALL CTGSYL( 'C', IJB, N1, N2, A, LDA, A( I, I ), LDA, 579: $ WORK, N1, B, LDB, B( I, I ), LDB, 580: $ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ), 581: $ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK, 582: $ IERR ) 583: END IF 584: GO TO 40 585: END IF 586: DIF( 1 ) = DSCALE / DIF( 1 ) 587: * 588: * 1-norm-based estimate of Difl. 589: * 590: 50 CONTINUE 591: CALL CLACN2( MN2, WORK( MN2+1 ), WORK, DIF( 2 ), KASE, 592: $ ISAVE ) 593: IF( KASE.NE.0 ) THEN 594: IF( KASE.EQ.1 ) THEN 595: * 596: * Solve generalized Sylvester equation 597: * 598: CALL CTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA, 599: $ WORK, N2, B( I, I ), LDB, B, LDB, 600: $ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ), 601: $ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK, 602: $ IERR ) 603: ELSE 604: * 605: * Solve the transposed variant. 606: * 607: CALL CTGSYL( 'C', IJB, N2, N1, A( I, I ), LDA, A, LDA, 608: $ WORK, N2, B, LDB, B( I, I ), LDB, 609: $ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ), 610: $ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK, 611: $ IERR ) 612: END IF 613: GO TO 50 614: END IF 615: DIF( 2 ) = DSCALE / DIF( 2 ) 616: END IF 617: END IF 618: * 619: * If B(K,K) is complex, make it real and positive (normalization 620: * of the generalized Schur form) and Store the generalized 621: * eigenvalues of reordered pair (A, B) 622: * 623: DO 60 K = 1, N 624: DSCALE = ABS( B( K, K ) ) 625: IF( DSCALE.GT.SAFMIN ) THEN 626: TEMP1 = CONJG( B( K, K ) / DSCALE ) 627: TEMP2 = B( K, K ) / DSCALE 628: B( K, K ) = DSCALE 629: CALL CSCAL( N-K, TEMP1, B( K, K+1 ), LDB ) 630: CALL CSCAL( N-K+1, TEMP1, A( K, K ), LDA ) 631: IF( WANTQ ) 632: $ CALL CSCAL( N, TEMP2, Q( 1, K ), 1 ) 633: ELSE 634: B( K, K ) = CMPLX( ZERO, ZERO ) 635: END IF 636: * 637: ALPHA( K ) = A( K, K ) 638: BETA( K ) = B( K, K ) 639: * 640: 60 CONTINUE 641: * 642: 70 CONTINUE 643: * 644: WORK( 1 ) = LWMIN 645: IWORK( 1 ) = LIWMIN 646: * 647: RETURN 648: * 649: * End of CTGSEN 650: * 651: END 652: