1: *> \brief \b ZTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an unitary equivalence transformation.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZTGEX2 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgex2.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgex2.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgex2.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
22: * LDZ, J1, INFO )
23: *
24: * .. Scalar Arguments ..
25: * LOGICAL WANTQ, WANTZ
26: * INTEGER INFO, J1, LDA, LDB, LDQ, LDZ, N
27: * ..
28: * .. Array Arguments ..
29: * COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
30: * $ Z( LDZ, * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> ZTGEX2 swaps adjacent diagonal 1 by 1 blocks (A11,B11) and (A22,B22)
40: *> in an upper triangular matrix pair (A, B) by an unitary equivalence
41: *> transformation.
42: *>
43: *> (A, B) must be in generalized Schur canonical form, that is, A and
44: *> B are both upper triangular.
45: *>
46: *> Optionally, the matrices Q and Z of generalized Schur vectors are
47: *> updated.
48: *>
49: *> Q(in) * A(in) * Z(in)**H = Q(out) * A(out) * Z(out)**H
50: *> Q(in) * B(in) * Z(in)**H = Q(out) * B(out) * Z(out)**H
51: *>
52: *> \endverbatim
53: *
54: * Arguments:
55: * ==========
56: *
57: *> \param[in] WANTQ
58: *> \verbatim
59: *> WANTQ is LOGICAL
60: *> .TRUE. : update the left transformation matrix Q;
61: *> .FALSE.: do not update Q.
62: *> \endverbatim
63: *>
64: *> \param[in] WANTZ
65: *> \verbatim
66: *> WANTZ is LOGICAL
67: *> .TRUE. : update the right transformation matrix Z;
68: *> .FALSE.: do not update Z.
69: *> \endverbatim
70: *>
71: *> \param[in] N
72: *> \verbatim
73: *> N is INTEGER
74: *> The order of the matrices A and B. N >= 0.
75: *> \endverbatim
76: *>
77: *> \param[in,out] A
78: *> \verbatim
79: *> A is COMPLEX*16 array, dimensions (LDA,N)
80: *> On entry, the matrix A in the pair (A, B).
81: *> On exit, the updated matrix A.
82: *> \endverbatim
83: *>
84: *> \param[in] LDA
85: *> \verbatim
86: *> LDA is INTEGER
87: *> The leading dimension of the array A. LDA >= max(1,N).
88: *> \endverbatim
89: *>
90: *> \param[in,out] B
91: *> \verbatim
92: *> B is COMPLEX*16 array, dimensions (LDB,N)
93: *> On entry, the matrix B in the pair (A, B).
94: *> On exit, the updated matrix B.
95: *> \endverbatim
96: *>
97: *> \param[in] LDB
98: *> \verbatim
99: *> LDB is INTEGER
100: *> The leading dimension of the array B. LDB >= max(1,N).
101: *> \endverbatim
102: *>
103: *> \param[in,out] Q
104: *> \verbatim
105: *> Q is COMPLEX*16 array, dimension (LDQ,N)
106: *> If WANTQ = .TRUE, on entry, the unitary matrix Q. On exit,
107: *> the updated matrix Q.
108: *> Not referenced if WANTQ = .FALSE..
109: *> \endverbatim
110: *>
111: *> \param[in] LDQ
112: *> \verbatim
113: *> LDQ is INTEGER
114: *> The leading dimension of the array Q. LDQ >= 1;
115: *> If WANTQ = .TRUE., LDQ >= N.
116: *> \endverbatim
117: *>
118: *> \param[in,out] Z
119: *> \verbatim
120: *> Z is COMPLEX*16 array, dimension (LDZ,N)
121: *> If WANTZ = .TRUE, on entry, the unitary matrix Z. On exit,
122: *> the updated matrix Z.
123: *> Not referenced if WANTZ = .FALSE..
124: *> \endverbatim
125: *>
126: *> \param[in] LDZ
127: *> \verbatim
128: *> LDZ is INTEGER
129: *> The leading dimension of the array Z. LDZ >= 1;
130: *> If WANTZ = .TRUE., LDZ >= N.
131: *> \endverbatim
132: *>
133: *> \param[in] J1
134: *> \verbatim
135: *> J1 is INTEGER
136: *> The index to the first block (A11, B11).
137: *> \endverbatim
138: *>
139: *> \param[out] INFO
140: *> \verbatim
141: *> INFO is INTEGER
142: *> =0: Successful exit.
143: *> =1: The transformed matrix pair (A, B) would be too far
144: *> from generalized Schur form; the problem is ill-
145: *> conditioned.
146: *> \endverbatim
147: *
148: * Authors:
149: * ========
150: *
151: *> \author Univ. of Tennessee
152: *> \author Univ. of California Berkeley
153: *> \author Univ. of Colorado Denver
154: *> \author NAG Ltd.
155: *
156: *> \ingroup complex16GEauxiliary
157: *
158: *> \par Further Details:
159: * =====================
160: *>
161: *> In the current code both weak and strong stability tests are
162: *> performed. The user can omit the strong stability test by changing
163: *> the internal logical parameter WANDS to .FALSE.. See ref. [2] for
164: *> details.
165: *
166: *> \par Contributors:
167: * ==================
168: *>
169: *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
170: *> Umea University, S-901 87 Umea, Sweden.
171: *
172: *> \par References:
173: * ================
174: *>
175: *> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
176: *> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
177: *> M.S. Moonen et al (eds), Linear Algebra for Large Scale and
178: *> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
179: *> \n
180: *> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
181: *> Eigenvalues of a Regular Matrix Pair (A, B) and Condition
182: *> Estimation: Theory, Algorithms and Software, Report UMINF-94.04,
183: *> Department of Computing Science, Umea University, S-901 87 Umea,
184: *> Sweden, 1994. Also as LAPACK Working Note 87. To appear in
185: *> Numerical Algorithms, 1996.
186: *>
187: * =====================================================================
188: SUBROUTINE ZTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
189: $ LDZ, J1, INFO )
190: *
191: * -- LAPACK auxiliary routine --
192: * -- LAPACK is a software package provided by Univ. of Tennessee, --
193: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
194: *
195: * .. Scalar Arguments ..
196: LOGICAL WANTQ, WANTZ
197: INTEGER INFO, J1, LDA, LDB, LDQ, LDZ, N
198: * ..
199: * .. Array Arguments ..
200: COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
201: $ Z( LDZ, * )
202: * ..
203: *
204: * =====================================================================
205: *
206: * .. Parameters ..
207: COMPLEX*16 CZERO, CONE
208: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
209: $ CONE = ( 1.0D+0, 0.0D+0 ) )
210: DOUBLE PRECISION TWENTY
211: PARAMETER ( TWENTY = 2.0D+1 )
212: INTEGER LDST
213: PARAMETER ( LDST = 2 )
214: LOGICAL WANDS
215: PARAMETER ( WANDS = .TRUE. )
216: * ..
217: * .. Local Scalars ..
218: LOGICAL STRONG, WEAK
219: INTEGER I, M
220: DOUBLE PRECISION CQ, CZ, EPS, SA, SB, SCALE, SMLNUM, SUM,
221: $ THRESHA, THRESHB
222: COMPLEX*16 CDUM, F, G, SQ, SZ
223: * ..
224: * .. Local Arrays ..
225: COMPLEX*16 S( LDST, LDST ), T( LDST, LDST ), WORK( 8 )
226: * ..
227: * .. External Functions ..
228: DOUBLE PRECISION DLAMCH
229: EXTERNAL DLAMCH
230: * ..
231: * .. External Subroutines ..
232: EXTERNAL ZLACPY, ZLARTG, ZLASSQ, ZROT
233: * ..
234: * .. Intrinsic Functions ..
235: INTRINSIC ABS, DBLE, DCONJG, MAX, SQRT
236: * ..
237: * .. Executable Statements ..
238: *
239: INFO = 0
240: *
241: * Quick return if possible
242: *
243: IF( N.LE.1 )
244: $ RETURN
245: *
246: M = LDST
247: WEAK = .FALSE.
248: STRONG = .FALSE.
249: *
250: * Make a local copy of selected block in (A, B)
251: *
252: CALL ZLACPY( 'Full', M, M, A( J1, J1 ), LDA, S, LDST )
253: CALL ZLACPY( 'Full', M, M, B( J1, J1 ), LDB, T, LDST )
254: *
255: * Compute the threshold for testing the acceptance of swapping.
256: *
257: EPS = DLAMCH( 'P' )
258: SMLNUM = DLAMCH( 'S' ) / EPS
259: SCALE = DBLE( CZERO )
260: SUM = DBLE( CONE )
261: CALL ZLACPY( 'Full', M, M, S, LDST, WORK, M )
262: CALL ZLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M )
263: CALL ZLASSQ( M*M, WORK, 1, SCALE, SUM )
264: SA = SCALE*SQRT( SUM )
265: SCALE = DBLE( CZERO )
266: SUM = DBLE( CONE )
267: CALL ZLASSQ( M*M, WORK(M*M+1), 1, SCALE, SUM )
268: SB = SCALE*SQRT( SUM )
269: *
270: * THRES has been changed from
271: * THRESH = MAX( TEN*EPS*SA, SMLNUM )
272: * to
273: * THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
274: * on 04/01/10.
275: * "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
276: * Jim Demmel and Guillaume Revy. See forum post 1783.
277: *
278: THRESHA = MAX( TWENTY*EPS*SA, SMLNUM )
279: THRESHB = MAX( TWENTY*EPS*SB, SMLNUM )
280: *
281: * Compute unitary QL and RQ that swap 1-by-1 and 1-by-1 blocks
282: * using Givens rotations and perform the swap tentatively.
283: *
284: F = S( 2, 2 )*T( 1, 1 ) - T( 2, 2 )*S( 1, 1 )
285: G = S( 2, 2 )*T( 1, 2 ) - T( 2, 2 )*S( 1, 2 )
286: SA = ABS( S( 2, 2 ) ) * ABS( T( 1, 1 ) )
287: SB = ABS( S( 1, 1 ) ) * ABS( T( 2, 2 ) )
288: CALL ZLARTG( G, F, CZ, SZ, CDUM )
289: SZ = -SZ
290: CALL ZROT( 2, S( 1, 1 ), 1, S( 1, 2 ), 1, CZ, DCONJG( SZ ) )
291: CALL ZROT( 2, T( 1, 1 ), 1, T( 1, 2 ), 1, CZ, DCONJG( SZ ) )
292: IF( SA.GE.SB ) THEN
293: CALL ZLARTG( S( 1, 1 ), S( 2, 1 ), CQ, SQ, CDUM )
294: ELSE
295: CALL ZLARTG( T( 1, 1 ), T( 2, 1 ), CQ, SQ, CDUM )
296: END IF
297: CALL ZROT( 2, S( 1, 1 ), LDST, S( 2, 1 ), LDST, CQ, SQ )
298: CALL ZROT( 2, T( 1, 1 ), LDST, T( 2, 1 ), LDST, CQ, SQ )
299: *
300: * Weak stability test: |S21| <= O(EPS F-norm((A)))
301: * and |T21| <= O(EPS F-norm((B)))
302: *
303: WEAK = ABS( S( 2, 1 ) ).LE.THRESHA .AND.
304: $ ABS( T( 2, 1 ) ).LE.THRESHB
305: IF( .NOT.WEAK )
306: $ GO TO 20
307: *
308: IF( WANDS ) THEN
309: *
310: * Strong stability test:
311: * F-norm((A-QL**H*S*QR)) <= O(EPS*F-norm((A)))
312: * and
313: * F-norm((B-QL**H*T*QR)) <= O(EPS*F-norm((B)))
314: *
315: CALL ZLACPY( 'Full', M, M, S, LDST, WORK, M )
316: CALL ZLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M )
317: CALL ZROT( 2, WORK, 1, WORK( 3 ), 1, CZ, -DCONJG( SZ ) )
318: CALL ZROT( 2, WORK( 5 ), 1, WORK( 7 ), 1, CZ, -DCONJG( SZ ) )
319: CALL ZROT( 2, WORK, 2, WORK( 2 ), 2, CQ, -SQ )
320: CALL ZROT( 2, WORK( 5 ), 2, WORK( 6 ), 2, CQ, -SQ )
321: DO 10 I = 1, 2
322: WORK( I ) = WORK( I ) - A( J1+I-1, J1 )
323: WORK( I+2 ) = WORK( I+2 ) - A( J1+I-1, J1+1 )
324: WORK( I+4 ) = WORK( I+4 ) - B( J1+I-1, J1 )
325: WORK( I+6 ) = WORK( I+6 ) - B( J1+I-1, J1+1 )
326: 10 CONTINUE
327: SCALE = DBLE( CZERO )
328: SUM = DBLE( CONE )
329: CALL ZLASSQ( M*M, WORK, 1, SCALE, SUM )
330: SA = SCALE*SQRT( SUM )
331: SCALE = DBLE( CZERO )
332: SUM = DBLE( CONE )
333: CALL ZLASSQ( M*M, WORK(M*M+1), 1, SCALE, SUM )
334: SB = SCALE*SQRT( SUM )
335: STRONG = SA.LE.THRESHA .AND. SB.LE.THRESHB
336: IF( .NOT.STRONG )
337: $ GO TO 20
338: END IF
339: *
340: * If the swap is accepted ("weakly" and "strongly"), apply the
341: * equivalence transformations to the original matrix pair (A,B)
342: *
343: CALL ZROT( J1+1, A( 1, J1 ), 1, A( 1, J1+1 ), 1, CZ,
344: $ DCONJG( SZ ) )
345: CALL ZROT( J1+1, B( 1, J1 ), 1, B( 1, J1+1 ), 1, CZ,
346: $ DCONJG( SZ ) )
347: CALL ZROT( N-J1+1, A( J1, J1 ), LDA, A( J1+1, J1 ), LDA, CQ, SQ )
348: CALL ZROT( N-J1+1, B( J1, J1 ), LDB, B( J1+1, J1 ), LDB, CQ, SQ )
349: *
350: * Set N1 by N2 (2,1) blocks to 0
351: *
352: A( J1+1, J1 ) = CZERO
353: B( J1+1, J1 ) = CZERO
354: *
355: * Accumulate transformations into Q and Z if requested.
356: *
357: IF( WANTZ )
358: $ CALL ZROT( N, Z( 1, J1 ), 1, Z( 1, J1+1 ), 1, CZ,
359: $ DCONJG( SZ ) )
360: IF( WANTQ )
361: $ CALL ZROT( N, Q( 1, J1 ), 1, Q( 1, J1+1 ), 1, CQ,
362: $ DCONJG( SQ ) )
363: *
364: * Exit with INFO = 0 if swap was successfully performed.
365: *
366: RETURN
367: *
368: * Exit with INFO = 1 if swap was rejected.
369: *
370: 20 CONTINUE
371: INFO = 1
372: RETURN
373: *
374: * End of ZTGEX2
375: *
376: 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