![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI, 2: $ M, S, SEP, WORK, LWORK, IWORK, LIWORK, 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: * .. Scalar Arguments .. 10: CHARACTER COMPQ, JOB 11: INTEGER INFO, LDQ, LDT, LIWORK, LWORK, M, N 12: DOUBLE PRECISION S, SEP 13: * .. 14: * .. Array Arguments .. 15: LOGICAL SELECT( * ) 16: INTEGER IWORK( * ) 17: DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( * ), 18: $ WR( * ) 19: * .. 20: * 21: * Purpose 22: * ======= 23: * 24: * DTRSEN reorders the real Schur factorization of a real matrix 25: * A = Q*T*Q**T, so that a selected cluster of eigenvalues appears in 26: * the leading diagonal blocks of the upper quasi-triangular matrix T, 27: * 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: * T must be in Schur canonical form (as returned by DHSEQR), that is, 34: * block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each 35: * 2-by-2 diagonal block has its diagonal elemnts equal and its 36: * off-diagonal elements of opposite sign. 37: * 38: * Arguments 39: * ========= 40: * 41: * JOB (input) CHARACTER*1 42: * Specifies whether condition numbers are required for the 43: * cluster of eigenvalues (S) or the invariant subspace (SEP): 44: * = 'N': none; 45: * = 'E': for eigenvalues only (S); 46: * = 'V': for invariant subspace only (SEP); 47: * = 'B': for both eigenvalues and invariant subspace (S and 48: * SEP). 49: * 50: * COMPQ (input) CHARACTER*1 51: * = 'V': update the matrix Q of Schur vectors; 52: * = 'N': do not update Q. 53: * 54: * SELECT (input) LOGICAL array, dimension (N) 55: * SELECT specifies the eigenvalues in the selected cluster. To 56: * select a real eigenvalue w(j), SELECT(j) must be set to 57: * .TRUE.. To select a complex conjugate pair of eigenvalues 58: * w(j) and w(j+1), corresponding to a 2-by-2 diagonal block, 59: * either SELECT(j) or SELECT(j+1) or both must be set to 60: * .TRUE.; a complex conjugate pair of eigenvalues must be 61: * either both included in the cluster or both excluded. 62: * 63: * N (input) INTEGER 64: * The order of the matrix T. N >= 0. 65: * 66: * T (input/output) DOUBLE PRECISION array, dimension (LDT,N) 67: * On entry, the upper quasi-triangular matrix T, in Schur 68: * canonical form. 69: * On exit, T is overwritten by the reordered matrix T, again in 70: * Schur canonical form, with the selected eigenvalues in the 71: * leading diagonal blocks. 72: * 73: * LDT (input) INTEGER 74: * The leading dimension of the array T. LDT >= max(1,N). 75: * 76: * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N) 77: * On entry, if COMPQ = 'V', the matrix Q of Schur vectors. 78: * On exit, if COMPQ = 'V', Q has been postmultiplied by the 79: * orthogonal transformation matrix which reorders T; the 80: * leading M columns of Q form an orthonormal basis for the 81: * specified invariant subspace. 82: * If COMPQ = 'N', Q is not referenced. 83: * 84: * LDQ (input) INTEGER 85: * The leading dimension of the array Q. 86: * LDQ >= 1; and if COMPQ = 'V', LDQ >= N. 87: * 88: * WR (output) DOUBLE PRECISION array, dimension (N) 89: * WI (output) DOUBLE PRECISION array, dimension (N) 90: * The real and imaginary parts, respectively, of the reordered 91: * eigenvalues of T. The eigenvalues are stored in the same 92: * order as on the diagonal of T, with WR(i) = T(i,i) and, if 93: * T(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) > 0 and 94: * WI(i+1) = -WI(i). Note that if a complex eigenvalue is 95: * sufficiently ill-conditioned, then its value may differ 96: * significantly from its value before reordering. 97: * 98: * M (output) INTEGER 99: * The dimension of the specified invariant subspace. 100: * 0 < = M <= N. 101: * 102: * S (output) DOUBLE PRECISION 103: * If JOB = 'E' or 'B', S is a lower bound on the reciprocal 104: * condition number for the selected cluster of eigenvalues. 105: * S cannot underestimate the true reciprocal condition number 106: * by more than a factor of sqrt(N). If M = 0 or N, S = 1. 107: * If JOB = 'N' or 'V', S is not referenced. 108: * 109: * SEP (output) DOUBLE PRECISION 110: * If JOB = 'V' or 'B', SEP is the estimated reciprocal 111: * condition number of the specified invariant subspace. If 112: * M = 0 or N, SEP = norm(T). 113: * If JOB = 'N' or 'E', SEP is not referenced. 114: * 115: * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 116: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 117: * 118: * LWORK (input) INTEGER 119: * The dimension of the array WORK. 120: * If JOB = 'N', LWORK >= max(1,N); 121: * if JOB = 'E', LWORK >= max(1,M*(N-M)); 122: * if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)). 123: * 124: * If LWORK = -1, then a workspace query is assumed; the routine 125: * only calculates the optimal size of the WORK array, returns 126: * this value as the first entry of the WORK array, and no error 127: * message related to LWORK is issued by XERBLA. 128: * 129: * IWORK (workspace) INTEGER array, dimension (MAX(1,LIWORK)) 130: * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 131: * 132: * LIWORK (input) INTEGER 133: * The dimension of the array IWORK. 134: * If JOB = 'N' or 'E', LIWORK >= 1; 135: * if JOB = 'V' or 'B', LIWORK >= max(1,M*(N-M)). 136: * 137: * If LIWORK = -1, then a workspace query is assumed; the 138: * routine only calculates the optimal size of the IWORK array, 139: * returns this value as the first entry of the IWORK array, and 140: * no error message related to LIWORK is issued by XERBLA. 141: * 142: * INFO (output) INTEGER 143: * = 0: successful exit 144: * < 0: if INFO = -i, the i-th argument had an illegal value 145: * = 1: reordering of T failed because some eigenvalues are too 146: * close to separate (the problem is very ill-conditioned); 147: * T may have been partially reordered, and WR and WI 148: * contain the eigenvalues in the same order as in T; S and 149: * SEP (if requested) are set to zero. 150: * 151: * Further Details 152: * =============== 153: * 154: * DTRSEN first collects the selected eigenvalues by computing an 155: * orthogonal transformation Z to move them to the top left corner of T. 156: * In other words, the selected eigenvalues are the eigenvalues of T11 157: * in: 158: * 159: * Z'*T*Z = ( T11 T12 ) n1 160: * ( 0 T22 ) n2 161: * n1 n2 162: * 163: * where N = n1+n2 and Z' means the transpose of Z. The first n1 columns 164: * of Z span the specified invariant subspace of T. 165: * 166: * If T has been obtained from the real Schur factorization of a matrix 167: * A = Q*T*Q', then the reordered real Schur factorization of A is given 168: * by A = (Q*Z)*(Z'*T*Z)*(Q*Z)', and the first n1 columns of Q*Z span 169: * the corresponding invariant subspace of A. 170: * 171: * The reciprocal condition number of the average of the eigenvalues of 172: * T11 may be returned in S. S lies between 0 (very badly conditioned) 173: * and 1 (very well conditioned). It is computed as follows. First we 174: * compute R so that 175: * 176: * P = ( I R ) n1 177: * ( 0 0 ) n2 178: * n1 n2 179: * 180: * is the projector on the invariant subspace associated with T11. 181: * R is the solution of the Sylvester equation: 182: * 183: * T11*R - R*T22 = T12. 184: * 185: * Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote 186: * the two-norm of M. Then S is computed as the lower bound 187: * 188: * (1 + F-norm(R)**2)**(-1/2) 189: * 190: * on the reciprocal of 2-norm(P), the true reciprocal condition number. 191: * S cannot underestimate 1 / 2-norm(P) by more than a factor of 192: * sqrt(N). 193: * 194: * An approximate error bound for the computed average of the 195: * eigenvalues of T11 is 196: * 197: * EPS * norm(T) / S 198: * 199: * where EPS is the machine precision. 200: * 201: * The reciprocal condition number of the right invariant subspace 202: * spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP. 203: * SEP is defined as the separation of T11 and T22: 204: * 205: * sep( T11, T22 ) = sigma-min( C ) 206: * 207: * where sigma-min(C) is the smallest singular value of the 208: * n1*n2-by-n1*n2 matrix 209: * 210: * C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) ) 211: * 212: * I(m) is an m by m identity matrix, and kprod denotes the Kronecker 213: * product. We estimate sigma-min(C) by the reciprocal of an estimate of 214: * the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C) 215: * cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2). 216: * 217: * When SEP is small, small changes in T can cause large changes in 218: * the invariant subspace. An approximate bound on the maximum angular 219: * error in the computed right invariant subspace is 220: * 221: * EPS * norm(T) / SEP 222: * 223: * ===================================================================== 224: * 225: * .. Parameters .. 226: DOUBLE PRECISION ZERO, ONE 227: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 228: * .. 229: * .. Local Scalars .. 230: LOGICAL LQUERY, PAIR, SWAP, WANTBH, WANTQ, WANTS, 231: $ WANTSP 232: INTEGER IERR, K, KASE, KK, KS, LIWMIN, LWMIN, N1, N2, 233: $ NN 234: DOUBLE PRECISION EST, RNORM, SCALE 235: * .. 236: * .. Local Arrays .. 237: INTEGER ISAVE( 3 ) 238: * .. 239: * .. External Functions .. 240: LOGICAL LSAME 241: DOUBLE PRECISION DLANGE 242: EXTERNAL LSAME, DLANGE 243: * .. 244: * .. External Subroutines .. 245: EXTERNAL DLACN2, DLACPY, DTREXC, DTRSYL, XERBLA 246: * .. 247: * .. Intrinsic Functions .. 248: INTRINSIC ABS, MAX, SQRT 249: * .. 250: * .. Executable Statements .. 251: * 252: * Decode and test the input parameters 253: * 254: WANTBH = LSAME( JOB, 'B' ) 255: WANTS = LSAME( JOB, 'E' ) .OR. WANTBH 256: WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH 257: WANTQ = LSAME( COMPQ, 'V' ) 258: * 259: INFO = 0 260: LQUERY = ( LWORK.EQ.-1 ) 261: IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.WANTS .AND. .NOT.WANTSP ) 262: $ THEN 263: INFO = -1 264: ELSE IF( .NOT.LSAME( COMPQ, 'N' ) .AND. .NOT.WANTQ ) THEN 265: INFO = -2 266: ELSE IF( N.LT.0 ) THEN 267: INFO = -4 268: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN 269: INFO = -6 270: ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 271: INFO = -8 272: ELSE 273: * 274: * Set M to the dimension of the specified invariant subspace, 275: * and test LWORK and LIWORK. 276: * 277: M = 0 278: PAIR = .FALSE. 279: DO 10 K = 1, N 280: IF( PAIR ) THEN 281: PAIR = .FALSE. 282: ELSE 283: IF( K.LT.N ) THEN 284: IF( T( K+1, K ).EQ.ZERO ) THEN 285: IF( SELECT( K ) ) 286: $ M = M + 1 287: ELSE 288: PAIR = .TRUE. 289: IF( SELECT( K ) .OR. SELECT( K+1 ) ) 290: $ M = M + 2 291: END IF 292: ELSE 293: IF( SELECT( N ) ) 294: $ M = M + 1 295: END IF 296: END IF 297: 10 CONTINUE 298: * 299: N1 = M 300: N2 = N - M 301: NN = N1*N2 302: * 303: IF( WANTSP ) THEN 304: LWMIN = MAX( 1, 2*NN ) 305: LIWMIN = MAX( 1, NN ) 306: ELSE IF( LSAME( JOB, 'N' ) ) THEN 307: LWMIN = MAX( 1, N ) 308: LIWMIN = 1 309: ELSE IF( LSAME( JOB, 'E' ) ) THEN 310: LWMIN = MAX( 1, NN ) 311: LIWMIN = 1 312: END IF 313: * 314: IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 315: INFO = -15 316: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 317: INFO = -17 318: END IF 319: END IF 320: * 321: IF( INFO.EQ.0 ) THEN 322: WORK( 1 ) = LWMIN 323: IWORK( 1 ) = LIWMIN 324: END IF 325: * 326: IF( INFO.NE.0 ) THEN 327: CALL XERBLA( 'DTRSEN', -INFO ) 328: RETURN 329: ELSE IF( LQUERY ) THEN 330: RETURN 331: END IF 332: * 333: * Quick return if possible. 334: * 335: IF( M.EQ.N .OR. M.EQ.0 ) THEN 336: IF( WANTS ) 337: $ S = ONE 338: IF( WANTSP ) 339: $ SEP = DLANGE( '1', N, N, T, LDT, WORK ) 340: GO TO 40 341: END IF 342: * 343: * Collect the selected blocks at the top-left corner of T. 344: * 345: KS = 0 346: PAIR = .FALSE. 347: DO 20 K = 1, N 348: IF( PAIR ) THEN 349: PAIR = .FALSE. 350: ELSE 351: SWAP = SELECT( K ) 352: IF( K.LT.N ) THEN 353: IF( T( K+1, K ).NE.ZERO ) THEN 354: PAIR = .TRUE. 355: SWAP = SWAP .OR. SELECT( K+1 ) 356: END IF 357: END IF 358: IF( SWAP ) THEN 359: KS = KS + 1 360: * 361: * Swap the K-th block to position KS. 362: * 363: IERR = 0 364: KK = K 365: IF( K.NE.KS ) 366: $ CALL DTREXC( COMPQ, N, T, LDT, Q, LDQ, KK, KS, WORK, 367: $ IERR ) 368: IF( IERR.EQ.1 .OR. IERR.EQ.2 ) THEN 369: * 370: * Blocks too close to swap: exit. 371: * 372: INFO = 1 373: IF( WANTS ) 374: $ S = ZERO 375: IF( WANTSP ) 376: $ SEP = ZERO 377: GO TO 40 378: END IF 379: IF( PAIR ) 380: $ KS = KS + 1 381: END IF 382: END IF 383: 20 CONTINUE 384: * 385: IF( WANTS ) THEN 386: * 387: * Solve Sylvester equation for R: 388: * 389: * T11*R - R*T22 = scale*T12 390: * 391: CALL DLACPY( 'F', N1, N2, T( 1, N1+1 ), LDT, WORK, N1 ) 392: CALL DTRSYL( 'N', 'N', -1, N1, N2, T, LDT, T( N1+1, N1+1 ), 393: $ LDT, WORK, N1, SCALE, IERR ) 394: * 395: * Estimate the reciprocal of the condition number of the cluster 396: * of eigenvalues. 397: * 398: RNORM = DLANGE( 'F', N1, N2, WORK, N1, WORK ) 399: IF( RNORM.EQ.ZERO ) THEN 400: S = ONE 401: ELSE 402: S = SCALE / ( SQRT( SCALE*SCALE / RNORM+RNORM )* 403: $ SQRT( RNORM ) ) 404: END IF 405: END IF 406: * 407: IF( WANTSP ) THEN 408: * 409: * Estimate sep(T11,T22). 410: * 411: EST = ZERO 412: KASE = 0 413: 30 CONTINUE 414: CALL DLACN2( NN, WORK( NN+1 ), WORK, IWORK, EST, KASE, ISAVE ) 415: IF( KASE.NE.0 ) THEN 416: IF( KASE.EQ.1 ) THEN 417: * 418: * Solve T11*R - R*T22 = scale*X. 419: * 420: CALL DTRSYL( 'N', 'N', -1, N1, N2, T, LDT, 421: $ T( N1+1, N1+1 ), LDT, WORK, N1, SCALE, 422: $ IERR ) 423: ELSE 424: * 425: * Solve T11'*R - R*T22' = scale*X. 426: * 427: CALL DTRSYL( 'T', 'T', -1, N1, N2, T, LDT, 428: $ T( N1+1, N1+1 ), LDT, WORK, N1, SCALE, 429: $ IERR ) 430: END IF 431: GO TO 30 432: END IF 433: * 434: SEP = SCALE / EST 435: END IF 436: * 437: 40 CONTINUE 438: * 439: * Store the output eigenvalues in WR and WI. 440: * 441: DO 50 K = 1, N 442: WR( K ) = T( K, K ) 443: WI( K ) = ZERO 444: 50 CONTINUE 445: DO 60 K = 1, N - 1 446: IF( T( K+1, K ).NE.ZERO ) THEN 447: WI( K ) = SQRT( ABS( T( K, K+1 ) ) )* 448: $ SQRT( ABS( T( K+1, K ) ) ) 449: WI( K+1 ) = -WI( K ) 450: END IF 451: 60 CONTINUE 452: * 453: WORK( 1 ) = LWMIN 454: IWORK( 1 ) = LIWMIN 455: * 456: RETURN 457: * 458: * End of DTRSEN 459: * 460: END