LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
ztgex2.f
Go to the documentation of this file.
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
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 tgex2
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
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition zlartg.f90:116
subroutine zlassq(n, x, incx, scale, sumsq)
ZLASSQ updates a sum of squares represented in scaled form.
Definition zlassq.f90:124
subroutine zrot(n, cx, incx, cy, incy, c, s)
ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
Definition zrot.f:103
subroutine ztgex2(wantq, wantz, n, a, lda, b, ldb, q, ldq, z, ldz, j1, info)
ZTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an unitary equiva...
Definition ztgex2.f:190