Annotation of rpl/lapack/lapack/ztgex2.f, revision 1.8
1.1 bertrand 1: SUBROUTINE ZTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
2: $ LDZ, J1, INFO )
3: *
1.5 bertrand 4: * -- LAPACK auxiliary routine (version 3.2.2) --
1.1 bertrand 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.5 bertrand 7: * June 2010
1.1 bertrand 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 ) )
1.5 bertrand 121: DOUBLE PRECISION TWENTY
122: PARAMETER ( TWENTY = 2.0D+1 )
1.1 bertrand 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 )
1.5 bertrand 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 )
1.1 bertrand 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztgex2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack