![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S, 2: $ SEP, WORK, LWORK, INFO ) 3: * 4: * -- LAPACK routine (version 3.2) -- 5: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 7: * November 2006 8: * 9: * Modified to call ZLACN2 in place of ZLACON, 10 Feb 03, SJH. 10: * 11: * .. Scalar Arguments .. 12: CHARACTER COMPQ, JOB 13: INTEGER INFO, LDQ, LDT, LWORK, M, N 14: DOUBLE PRECISION S, SEP 15: * .. 16: * .. Array Arguments .. 17: LOGICAL SELECT( * ) 18: COMPLEX*16 Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * ) 19: * .. 20: * 21: * Purpose 22: * ======= 23: * 24: * ZTRSEN reorders the Schur factorization of a complex matrix 25: * A = Q*T*Q**H, so that a selected cluster of eigenvalues appears in 26: * the leading positions on the diagonal of the upper triangular matrix 27: * T, and the leading columns of Q form an orthonormal basis of the 28: * corresponding right invariant subspace. 29: * 30: * Optionally the routine computes the reciprocal condition numbers of 31: * the cluster of eigenvalues and/or the invariant subspace. 32: * 33: * Arguments 34: * ========= 35: * 36: * JOB (input) CHARACTER*1 37: * Specifies whether condition numbers are required for the 38: * cluster of eigenvalues (S) or the invariant subspace (SEP): 39: * = 'N': none; 40: * = 'E': for eigenvalues only (S); 41: * = 'V': for invariant subspace only (SEP); 42: * = 'B': for both eigenvalues and invariant subspace (S and 43: * SEP). 44: * 45: * COMPQ (input) CHARACTER*1 46: * = 'V': update the matrix Q of Schur vectors; 47: * = 'N': do not update Q. 48: * 49: * SELECT (input) LOGICAL array, dimension (N) 50: * SELECT specifies the eigenvalues in the selected cluster. To 51: * select the j-th eigenvalue, SELECT(j) must be set to .TRUE.. 52: * 53: * N (input) INTEGER 54: * The order of the matrix T. N >= 0. 55: * 56: * T (input/output) COMPLEX*16 array, dimension (LDT,N) 57: * On entry, the upper triangular matrix T. 58: * On exit, T is overwritten by the reordered matrix T, with the 59: * selected eigenvalues as the leading diagonal elements. 60: * 61: * LDT (input) INTEGER 62: * The leading dimension of the array T. LDT >= max(1,N). 63: * 64: * Q (input/output) COMPLEX*16 array, dimension (LDQ,N) 65: * On entry, if COMPQ = 'V', the matrix Q of Schur vectors. 66: * On exit, if COMPQ = 'V', Q has been postmultiplied by the 67: * unitary transformation matrix which reorders T; the leading M 68: * columns of Q form an orthonormal basis for the specified 69: * invariant subspace. 70: * If COMPQ = 'N', Q is not referenced. 71: * 72: * LDQ (input) INTEGER 73: * The leading dimension of the array Q. 74: * LDQ >= 1; and if COMPQ = 'V', LDQ >= N. 75: * 76: * W (output) COMPLEX*16 array, dimension (N) 77: * The reordered eigenvalues of T, in the same order as they 78: * appear on the diagonal of T. 79: * 80: * M (output) INTEGER 81: * The dimension of the specified invariant subspace. 82: * 0 <= M <= N. 83: * 84: * S (output) DOUBLE PRECISION 85: * If JOB = 'E' or 'B', S is a lower bound on the reciprocal 86: * condition number for the selected cluster of eigenvalues. 87: * S cannot underestimate the true reciprocal condition number 88: * by more than a factor of sqrt(N). If M = 0 or N, S = 1. 89: * If JOB = 'N' or 'V', S is not referenced. 90: * 91: * SEP (output) DOUBLE PRECISION 92: * If JOB = 'V' or 'B', SEP is the estimated reciprocal 93: * condition number of the specified invariant subspace. If 94: * M = 0 or N, SEP = norm(T). 95: * If JOB = 'N' or 'E', SEP is not referenced. 96: * 97: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) 98: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 99: * 100: * LWORK (input) INTEGER 101: * The dimension of the array WORK. 102: * If JOB = 'N', LWORK >= 1; 103: * if JOB = 'E', LWORK = max(1,M*(N-M)); 104: * if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)). 105: * 106: * If LWORK = -1, then a workspace query is assumed; the routine 107: * only calculates the optimal size of the WORK array, returns 108: * this value as the first entry of the WORK array, and no error 109: * message related to LWORK is issued by XERBLA. 110: * 111: * INFO (output) INTEGER 112: * = 0: successful exit 113: * < 0: if INFO = -i, the i-th argument had an illegal value 114: * 115: * Further Details 116: * =============== 117: * 118: * ZTRSEN first collects the selected eigenvalues by computing a unitary 119: * transformation Z to move them to the top left corner of T. In other 120: * words, the selected eigenvalues are the eigenvalues of T11 in: 121: * 122: * Z'*T*Z = ( T11 T12 ) n1 123: * ( 0 T22 ) n2 124: * n1 n2 125: * 126: * where N = n1+n2 and Z' means the conjugate transpose of Z. The first 127: * n1 columns of Z span the specified invariant subspace of T. 128: * 129: * If T has been obtained from the Schur factorization of a matrix 130: * A = Q*T*Q', then the reordered Schur factorization of A is given by 131: * A = (Q*Z)*(Z'*T*Z)*(Q*Z)', and the first n1 columns of Q*Z span the 132: * corresponding invariant subspace of A. 133: * 134: * The reciprocal condition number of the average of the eigenvalues of 135: * T11 may be returned in S. S lies between 0 (very badly conditioned) 136: * and 1 (very well conditioned). It is computed as follows. First we 137: * compute R so that 138: * 139: * P = ( I R ) n1 140: * ( 0 0 ) n2 141: * n1 n2 142: * 143: * is the projector on the invariant subspace associated with T11. 144: * R is the solution of the Sylvester equation: 145: * 146: * T11*R - R*T22 = T12. 147: * 148: * Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote 149: * the two-norm of M. Then S is computed as the lower bound 150: * 151: * (1 + F-norm(R)**2)**(-1/2) 152: * 153: * on the reciprocal of 2-norm(P), the true reciprocal condition number. 154: * S cannot underestimate 1 / 2-norm(P) by more than a factor of 155: * sqrt(N). 156: * 157: * An approximate error bound for the computed average of the 158: * eigenvalues of T11 is 159: * 160: * EPS * norm(T) / S 161: * 162: * where EPS is the machine precision. 163: * 164: * The reciprocal condition number of the right invariant subspace 165: * spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP. 166: * SEP is defined as the separation of T11 and T22: 167: * 168: * sep( T11, T22 ) = sigma-min( C ) 169: * 170: * where sigma-min(C) is the smallest singular value of the 171: * n1*n2-by-n1*n2 matrix 172: * 173: * C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) ) 174: * 175: * I(m) is an m by m identity matrix, and kprod denotes the Kronecker 176: * product. We estimate sigma-min(C) by the reciprocal of an estimate of 177: * the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C) 178: * cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2). 179: * 180: * When SEP is small, small changes in T can cause large changes in 181: * the invariant subspace. An approximate bound on the maximum angular 182: * error in the computed right invariant subspace is 183: * 184: * EPS * norm(T) / SEP 185: * 186: * ===================================================================== 187: * 188: * .. Parameters .. 189: DOUBLE PRECISION ZERO, ONE 190: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 191: * .. 192: * .. Local Scalars .. 193: LOGICAL LQUERY, WANTBH, WANTQ, WANTS, WANTSP 194: INTEGER IERR, K, KASE, KS, LWMIN, N1, N2, NN 195: DOUBLE PRECISION EST, RNORM, SCALE 196: * .. 197: * .. Local Arrays .. 198: INTEGER ISAVE( 3 ) 199: DOUBLE PRECISION RWORK( 1 ) 200: * .. 201: * .. External Functions .. 202: LOGICAL LSAME 203: DOUBLE PRECISION ZLANGE 204: EXTERNAL LSAME, ZLANGE 205: * .. 206: * .. External Subroutines .. 207: EXTERNAL XERBLA, ZLACN2, ZLACPY, ZTREXC, ZTRSYL 208: * .. 209: * .. Intrinsic Functions .. 210: INTRINSIC MAX, SQRT 211: * .. 212: * .. Executable Statements .. 213: * 214: * Decode and test the input parameters. 215: * 216: WANTBH = LSAME( JOB, 'B' ) 217: WANTS = LSAME( JOB, 'E' ) .OR. WANTBH 218: WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH 219: WANTQ = LSAME( COMPQ, 'V' ) 220: * 221: * Set M to the number of selected eigenvalues. 222: * 223: M = 0 224: DO 10 K = 1, N 225: IF( SELECT( K ) ) 226: $ M = M + 1 227: 10 CONTINUE 228: * 229: N1 = M 230: N2 = N - M 231: NN = N1*N2 232: * 233: INFO = 0 234: LQUERY = ( LWORK.EQ.-1 ) 235: * 236: IF( WANTSP ) THEN 237: LWMIN = MAX( 1, 2*NN ) 238: ELSE IF( LSAME( JOB, 'N' ) ) THEN 239: LWMIN = 1 240: ELSE IF( LSAME( JOB, 'E' ) ) THEN 241: LWMIN = MAX( 1, NN ) 242: END IF 243: * 244: IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.WANTS .AND. .NOT.WANTSP ) 245: $ THEN 246: INFO = -1 247: ELSE IF( .NOT.LSAME( COMPQ, 'N' ) .AND. .NOT.WANTQ ) THEN 248: INFO = -2 249: ELSE IF( N.LT.0 ) THEN 250: INFO = -4 251: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN 252: INFO = -6 253: ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 254: INFO = -8 255: ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 256: INFO = -14 257: END IF 258: * 259: IF( INFO.EQ.0 ) THEN 260: WORK( 1 ) = LWMIN 261: END IF 262: * 263: IF( INFO.NE.0 ) THEN 264: CALL XERBLA( 'ZTRSEN', -INFO ) 265: RETURN 266: ELSE IF( LQUERY ) THEN 267: RETURN 268: END IF 269: * 270: * Quick return if possible 271: * 272: IF( M.EQ.N .OR. M.EQ.0 ) THEN 273: IF( WANTS ) 274: $ S = ONE 275: IF( WANTSP ) 276: $ SEP = ZLANGE( '1', N, N, T, LDT, RWORK ) 277: GO TO 40 278: END IF 279: * 280: * Collect the selected eigenvalues at the top left corner of T. 281: * 282: KS = 0 283: DO 20 K = 1, N 284: IF( SELECT( K ) ) THEN 285: KS = KS + 1 286: * 287: * Swap the K-th eigenvalue to position KS. 288: * 289: IF( K.NE.KS ) 290: $ CALL ZTREXC( COMPQ, N, T, LDT, Q, LDQ, K, KS, IERR ) 291: END IF 292: 20 CONTINUE 293: * 294: IF( WANTS ) THEN 295: * 296: * Solve the Sylvester equation for R: 297: * 298: * T11*R - R*T22 = scale*T12 299: * 300: CALL ZLACPY( 'F', N1, N2, T( 1, N1+1 ), LDT, WORK, N1 ) 301: CALL ZTRSYL( 'N', 'N', -1, N1, N2, T, LDT, T( N1+1, N1+1 ), 302: $ LDT, WORK, N1, SCALE, IERR ) 303: * 304: * Estimate the reciprocal of the condition number of the cluster 305: * of eigenvalues. 306: * 307: RNORM = ZLANGE( 'F', N1, N2, WORK, N1, RWORK ) 308: IF( RNORM.EQ.ZERO ) THEN 309: S = ONE 310: ELSE 311: S = SCALE / ( SQRT( SCALE*SCALE / RNORM+RNORM )* 312: $ SQRT( RNORM ) ) 313: END IF 314: END IF 315: * 316: IF( WANTSP ) THEN 317: * 318: * Estimate sep(T11,T22). 319: * 320: EST = ZERO 321: KASE = 0 322: 30 CONTINUE 323: CALL ZLACN2( NN, WORK( NN+1 ), WORK, EST, KASE, ISAVE ) 324: IF( KASE.NE.0 ) THEN 325: IF( KASE.EQ.1 ) THEN 326: * 327: * Solve T11*R - R*T22 = scale*X. 328: * 329: CALL ZTRSYL( 'N', 'N', -1, N1, N2, T, LDT, 330: $ T( N1+1, N1+1 ), LDT, WORK, N1, SCALE, 331: $ IERR ) 332: ELSE 333: * 334: * Solve T11'*R - R*T22' = scale*X. 335: * 336: CALL ZTRSYL( 'C', 'C', -1, N1, N2, T, LDT, 337: $ T( N1+1, N1+1 ), LDT, WORK, N1, SCALE, 338: $ IERR ) 339: END IF 340: GO TO 30 341: END IF 342: * 343: SEP = SCALE / EST 344: END IF 345: * 346: 40 CONTINUE 347: * 348: * Copy reordered eigenvalues to W. 349: * 350: DO 50 K = 1, N 351: W( K ) = T( K, K ) 352: 50 CONTINUE 353: * 354: WORK( 1 ) = LWMIN 355: * 356: RETURN 357: * 358: * End of ZTRSEN 359: * 360: END