![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 2: $ LDZ, J1, INFO ) 3: * 4: * -- LAPACK auxiliary routine (version 3.2.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: * June 2010 8: * 9: * .. Scalar Arguments .. 10: LOGICAL WANTQ, WANTZ 11: INTEGER INFO, J1, LDA, LDB, LDQ, LDZ, N 12: * .. 13: * .. Array Arguments .. 14: COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 15: $ Z( LDZ, * ) 16: * .. 17: * 18: * Purpose 19: * ======= 20: * 21: * ZTGEX2 swaps adjacent diagonal 1 by 1 blocks (A11,B11) and (A22,B22) 22: * in an upper triangular matrix pair (A, B) by an unitary equivalence 23: * transformation. 24: * 25: * (A, B) must be in generalized Schur canonical form, that is, A and 26: * B are both upper triangular. 27: * 28: * Optionally, the matrices Q and Z of generalized Schur vectors are 29: * updated. 30: * 31: * Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)' 32: * Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)' 33: * 34: * 35: * Arguments 36: * ========= 37: * 38: * WANTQ (input) LOGICAL 39: * .TRUE. : update the left transformation matrix Q; 40: * .FALSE.: do not update Q. 41: * 42: * WANTZ (input) LOGICAL 43: * .TRUE. : update the right transformation matrix Z; 44: * .FALSE.: do not update Z. 45: * 46: * N (input) INTEGER 47: * The order of the matrices A and B. N >= 0. 48: * 49: * A (input/output) COMPLEX*16 arrays, dimensions (LDA,N) 50: * On entry, the matrix A in the pair (A, B). 51: * On exit, the updated matrix A. 52: * 53: * LDA (input) INTEGER 54: * The leading dimension of the array A. LDA >= max(1,N). 55: * 56: * B (input/output) COMPLEX*16 arrays, dimensions (LDB,N) 57: * On entry, the matrix B in the pair (A, B). 58: * On exit, the updated matrix B. 59: * 60: * LDB (input) INTEGER 61: * The leading dimension of the array B. LDB >= max(1,N). 62: * 63: * Q (input/output) COMPLEX*16 array, dimension (LDZ,N) 64: * If WANTQ = .TRUE, on entry, the unitary matrix Q. On exit, 65: * the updated matrix Q. 66: * Not referenced if WANTQ = .FALSE.. 67: * 68: * LDQ (input) INTEGER 69: * The leading dimension of the array Q. LDQ >= 1; 70: * If WANTQ = .TRUE., LDQ >= N. 71: * 72: * Z (input/output) COMPLEX*16 array, dimension (LDZ,N) 73: * If WANTZ = .TRUE, on entry, the unitary matrix Z. On exit, 74: * the updated matrix Z. 75: * Not referenced if WANTZ = .FALSE.. 76: * 77: * LDZ (input) INTEGER 78: * The leading dimension of the array Z. LDZ >= 1; 79: * If WANTZ = .TRUE., LDZ >= N. 80: * 81: * J1 (input) INTEGER 82: * The index to the first block (A11, B11). 83: * 84: * INFO (output) INTEGER 85: * =0: Successful exit. 86: * =1: The transformed matrix pair (A, B) would be too far 87: * from generalized Schur form; the problem is ill- 88: * conditioned. 89: * 90: * 91: * Further Details 92: * =============== 93: * 94: * Based on contributions by 95: * Bo Kagstrom and Peter Poromaa, Department of Computing Science, 96: * Umea University, S-901 87 Umea, Sweden. 97: * 98: * In the current code both weak and strong stability tests are 99: * performed. The user can omit the strong stability test by changing 100: * the internal logical parameter WANDS to .FALSE.. See ref. [2] for 101: * details. 102: * 103: * [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the 104: * Generalized Real Schur Form of a Regular Matrix Pair (A, B), in 105: * M.S. Moonen et al (eds), Linear Algebra for Large Scale and 106: * Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218. 107: * 108: * [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified 109: * Eigenvalues of a Regular Matrix Pair (A, B) and Condition 110: * Estimation: Theory, Algorithms and Software, Report UMINF-94.04, 111: * Department of Computing Science, Umea University, S-901 87 Umea, 112: * Sweden, 1994. Also as LAPACK Working Note 87. To appear in 113: * Numerical Algorithms, 1996. 114: * 115: * ===================================================================== 116: * 117: * .. Parameters .. 118: COMPLEX*16 CZERO, CONE 119: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 120: $ CONE = ( 1.0D+0, 0.0D+0 ) ) 121: DOUBLE PRECISION TWENTY 122: PARAMETER ( TWENTY = 2.0D+1 ) 123: INTEGER LDST 124: PARAMETER ( LDST = 2 ) 125: LOGICAL WANDS 126: PARAMETER ( WANDS = .TRUE. ) 127: * .. 128: * .. Local Scalars .. 129: LOGICAL DTRONG, WEAK 130: INTEGER I, M 131: DOUBLE PRECISION CQ, CZ, EPS, SA, SB, SCALE, SMLNUM, SS, SUM, 132: $ THRESH, WS 133: COMPLEX*16 CDUM, F, G, SQ, SZ 134: * .. 135: * .. Local Arrays .. 136: COMPLEX*16 S( LDST, LDST ), T( LDST, LDST ), WORK( 8 ) 137: * .. 138: * .. External Functions .. 139: DOUBLE PRECISION DLAMCH 140: EXTERNAL DLAMCH 141: * .. 142: * .. External Subroutines .. 143: EXTERNAL ZLACPY, ZLARTG, ZLASSQ, ZROT 144: * .. 145: * .. Intrinsic Functions .. 146: INTRINSIC ABS, DBLE, DCONJG, MAX, SQRT 147: * .. 148: * .. Executable Statements .. 149: * 150: INFO = 0 151: * 152: * Quick return if possible 153: * 154: IF( N.LE.1 ) 155: $ RETURN 156: * 157: M = LDST 158: WEAK = .FALSE. 159: DTRONG = .FALSE. 160: * 161: * Make a local copy of selected block in (A, B) 162: * 163: CALL ZLACPY( 'Full', M, M, A( J1, J1 ), LDA, S, LDST ) 164: CALL ZLACPY( 'Full', M, M, B( J1, J1 ), LDB, T, LDST ) 165: * 166: * Compute the threshold for testing the acceptance of swapping. 167: * 168: EPS = DLAMCH( 'P' ) 169: SMLNUM = DLAMCH( 'S' ) / EPS 170: SCALE = DBLE( CZERO ) 171: SUM = DBLE( CONE ) 172: CALL ZLACPY( 'Full', M, M, S, LDST, WORK, M ) 173: CALL ZLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M ) 174: CALL ZLASSQ( 2*M*M, WORK, 1, SCALE, SUM ) 175: SA = SCALE*SQRT( SUM ) 176: * 177: * THRES has been changed from 178: * THRESH = MAX( TEN*EPS*SA, SMLNUM ) 179: * to 180: * THRESH = MAX( TWENTY*EPS*SA, SMLNUM ) 181: * on 04/01/10. 182: * "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by 183: * Jim Demmel and Guillaume Revy. See forum post 1783. 184: * 185: THRESH = MAX( TWENTY*EPS*SA, SMLNUM ) 186: * 187: * Compute unitary QL and RQ that swap 1-by-1 and 1-by-1 blocks 188: * using Givens rotations and perform the swap tentatively. 189: * 190: F = S( 2, 2 )*T( 1, 1 ) - T( 2, 2 )*S( 1, 1 ) 191: G = S( 2, 2 )*T( 1, 2 ) - T( 2, 2 )*S( 1, 2 ) 192: SA = ABS( S( 2, 2 ) ) 193: SB = ABS( T( 2, 2 ) ) 194: CALL ZLARTG( G, F, CZ, SZ, CDUM ) 195: SZ = -SZ 196: CALL ZROT( 2, S( 1, 1 ), 1, S( 1, 2 ), 1, CZ, DCONJG( SZ ) ) 197: CALL ZROT( 2, T( 1, 1 ), 1, T( 1, 2 ), 1, CZ, DCONJG( SZ ) ) 198: IF( SA.GE.SB ) THEN 199: CALL ZLARTG( S( 1, 1 ), S( 2, 1 ), CQ, SQ, CDUM ) 200: ELSE 201: CALL ZLARTG( T( 1, 1 ), T( 2, 1 ), CQ, SQ, CDUM ) 202: END IF 203: CALL ZROT( 2, S( 1, 1 ), LDST, S( 2, 1 ), LDST, CQ, SQ ) 204: CALL ZROT( 2, T( 1, 1 ), LDST, T( 2, 1 ), LDST, CQ, SQ ) 205: * 206: * Weak stability test: |S21| + |T21| <= O(EPS F-norm((S, T))) 207: * 208: WS = ABS( S( 2, 1 ) ) + ABS( T( 2, 1 ) ) 209: WEAK = WS.LE.THRESH 210: IF( .NOT.WEAK ) 211: $ GO TO 20 212: * 213: IF( WANDS ) THEN 214: * 215: * Strong stability test: 216: * F-norm((A-QL'*S*QR, B-QL'*T*QR)) <= O(EPS*F-norm((A, B))) 217: * 218: CALL ZLACPY( 'Full', M, M, S, LDST, WORK, M ) 219: CALL ZLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M ) 220: CALL ZROT( 2, WORK, 1, WORK( 3 ), 1, CZ, -DCONJG( SZ ) ) 221: CALL ZROT( 2, WORK( 5 ), 1, WORK( 7 ), 1, CZ, -DCONJG( SZ ) ) 222: CALL ZROT( 2, WORK, 2, WORK( 2 ), 2, CQ, -SQ ) 223: CALL ZROT( 2, WORK( 5 ), 2, WORK( 6 ), 2, CQ, -SQ ) 224: DO 10 I = 1, 2 225: WORK( I ) = WORK( I ) - A( J1+I-1, J1 ) 226: WORK( I+2 ) = WORK( I+2 ) - A( J1+I-1, J1+1 ) 227: WORK( I+4 ) = WORK( I+4 ) - B( J1+I-1, J1 ) 228: WORK( I+6 ) = WORK( I+6 ) - B( J1+I-1, J1+1 ) 229: 10 CONTINUE 230: SCALE = DBLE( CZERO ) 231: SUM = DBLE( CONE ) 232: CALL ZLASSQ( 2*M*M, WORK, 1, SCALE, SUM ) 233: SS = SCALE*SQRT( SUM ) 234: DTRONG = SS.LE.THRESH 235: IF( .NOT.DTRONG ) 236: $ GO TO 20 237: END IF 238: * 239: * If the swap is accepted ("weakly" and "strongly"), apply the 240: * equivalence transformations to the original matrix pair (A,B) 241: * 242: CALL ZROT( J1+1, A( 1, J1 ), 1, A( 1, J1+1 ), 1, CZ, 243: $ DCONJG( SZ ) ) 244: CALL ZROT( J1+1, B( 1, J1 ), 1, B( 1, J1+1 ), 1, CZ, 245: $ DCONJG( SZ ) ) 246: CALL ZROT( N-J1+1, A( J1, J1 ), LDA, A( J1+1, J1 ), LDA, CQ, SQ ) 247: CALL ZROT( N-J1+1, B( J1, J1 ), LDB, B( J1+1, J1 ), LDB, CQ, SQ ) 248: * 249: * Set N1 by N2 (2,1) blocks to 0 250: * 251: A( J1+1, J1 ) = CZERO 252: B( J1+1, J1 ) = CZERO 253: * 254: * Accumulate transformations into Q and Z if requested. 255: * 256: IF( WANTZ ) 257: $ CALL ZROT( N, Z( 1, J1 ), 1, Z( 1, J1+1 ), 1, CZ, 258: $ DCONJG( SZ ) ) 259: IF( WANTQ ) 260: $ CALL ZROT( N, Q( 1, J1 ), 1, Q( 1, J1+1 ), 1, CQ, 261: $ DCONJG( SQ ) ) 262: * 263: * Exit with INFO = 0 if swap was successfully performed. 264: * 265: RETURN 266: * 267: * Exit with INFO = 1 if swap was rejected. 268: * 269: 20 CONTINUE 270: INFO = 1 271: RETURN 272: * 273: * End of ZTGEX2 274: * 275: END