LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
ctgex2.f
Go to the documentation of this file.
1 *> \brief \b CTGEX2 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 CTGEX2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctgex2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctgex2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctgex2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CTGEX2( 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 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
30 * $ Z( LDZ, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> CTGEX2 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 arrays, 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 arrays, 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 array, dimension (LDZ,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 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 *> \date September 2012
157 *
158 *> \ingroup complexGEauxiliary
159 *
160 *> \par Further Details:
161 * =====================
162 *>
163 *> In the current code both weak and strong stability tests are
164 *> performed. The user can omit the strong stability test by changing
165 *> the internal logical parameter WANDS to .FALSE.. See ref. [2] for
166 *> details.
167 *
168 *> \par Contributors:
169 * ==================
170 *>
171 *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
172 *> Umea University, S-901 87 Umea, Sweden.
173 *
174 *> \par References:
175 * ================
176 *>
177 *> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
178 *> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
179 *> M.S. Moonen et al (eds), Linear Algebra for Large Scale and
180 *> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
181 *> \n
182 *> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
183 *> Eigenvalues of a Regular Matrix Pair (A, B) and Condition
184 *> Estimation: Theory, Algorithms and Software, Report UMINF-94.04,
185 *> Department of Computing Science, Umea University, S-901 87 Umea,
186 *> Sweden, 1994. Also as LAPACK Working Note 87. To appear in
187 *> Numerical Algorithms, 1996.
188 *>
189 * =====================================================================
190  SUBROUTINE ctgex2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
191  $ ldz, j1, info )
192 *
193 * -- LAPACK auxiliary routine (version 3.4.2) --
194 * -- LAPACK is a software package provided by Univ. of Tennessee, --
195 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
196 * September 2012
197 *
198 * .. Scalar Arguments ..
199  LOGICAL wantq, wantz
200  INTEGER info, j1, lda, ldb, ldq, ldz, n
201 * ..
202 * .. Array Arguments ..
203  COMPLEX a( lda, * ), b( ldb, * ), q( ldq, * ),
204  $ z( ldz, * )
205 * ..
206 *
207 * =====================================================================
208 *
209 * .. Parameters ..
210  COMPLEX czero, cone
211  parameter( czero = ( 0.0e+0, 0.0e+0 ),
212  $ cone = ( 1.0e+0, 0.0e+0 ) )
213  REAL twenty
214  parameter( twenty = 2.0e+1 )
215  INTEGER ldst
216  parameter( ldst = 2 )
217  LOGICAL wands
218  parameter( wands = .true. )
219 * ..
220 * .. Local Scalars ..
221  LOGICAL strong, weak
222  INTEGER i, m
223  REAL cq, cz, eps, sa, sb, scale, smlnum, ss, sum,
224  $ thresh, ws
225  COMPLEX cdum, f, g, sq, sz
226 * ..
227 * .. Local Arrays ..
228  COMPLEX s( ldst, ldst ), t( ldst, ldst ), work( 8 )
229 * ..
230 * .. External Functions ..
231  REAL slamch
232  EXTERNAL slamch
233 * ..
234 * .. External Subroutines ..
235  EXTERNAL clacpy, clartg, classq, crot
236 * ..
237 * .. Intrinsic Functions ..
238  INTRINSIC abs, conjg, max, REAL, sqrt
239 * ..
240 * .. Executable Statements ..
241 *
242  info = 0
243 *
244 * Quick return if possible
245 *
246  IF( n.LE.1 )
247  $ RETURN
248 *
249  m = ldst
250  weak = .false.
251  strong = .false.
252 *
253 * Make a local copy of selected block in (A, B)
254 *
255  CALL clacpy( 'Full', m, m, a( j1, j1 ), lda, s, ldst )
256  CALL clacpy( 'Full', m, m, b( j1, j1 ), ldb, t, ldst )
257 *
258 * Compute the threshold for testing the acceptance of swapping.
259 *
260  eps = slamch( 'P' )
261  smlnum = slamch( 'S' ) / eps
262  scale = REAL( czero )
263  sum = REAL( cone )
264  CALL clacpy( 'Full', m, m, s, ldst, work, m )
265  CALL clacpy( 'Full', m, m, t, ldst, work( m*m+1 ), m )
266  CALL classq( 2*m*m, work, 1, scale, sum )
267  sa = scale*sqrt( sum )
268 *
269 * THRES has been changed from
270 * THRESH = MAX( TEN*EPS*SA, SMLNUM )
271 * to
272 * THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
273 * on 04/01/10.
274 * "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
275 * Jim Demmel and Guillaume Revy. See forum post 1783.
276 *
277  thresh = max( twenty*eps*sa, smlnum )
278 *
279 * Compute unitary QL and RQ that swap 1-by-1 and 1-by-1 blocks
280 * using Givens rotations and perform the swap tentatively.
281 *
282  f = s( 2, 2 )*t( 1, 1 ) - t( 2, 2 )*s( 1, 1 )
283  g = s( 2, 2 )*t( 1, 2 ) - t( 2, 2 )*s( 1, 2 )
284  sa = abs( s( 2, 2 ) )
285  sb = abs( t( 2, 2 ) )
286  CALL clartg( g, f, cz, sz, cdum )
287  sz = -sz
288  CALL crot( 2, s( 1, 1 ), 1, s( 1, 2 ), 1, cz, conjg( sz ) )
289  CALL crot( 2, t( 1, 1 ), 1, t( 1, 2 ), 1, cz, conjg( sz ) )
290  IF( sa.GE.sb ) THEN
291  CALL clartg( s( 1, 1 ), s( 2, 1 ), cq, sq, cdum )
292  ELSE
293  CALL clartg( t( 1, 1 ), t( 2, 1 ), cq, sq, cdum )
294  END IF
295  CALL crot( 2, s( 1, 1 ), ldst, s( 2, 1 ), ldst, cq, sq )
296  CALL crot( 2, t( 1, 1 ), ldst, t( 2, 1 ), ldst, cq, sq )
297 *
298 * Weak stability test: |S21| + |T21| <= O(EPS F-norm((S, T)))
299 *
300  ws = abs( s( 2, 1 ) ) + abs( t( 2, 1 ) )
301  weak = ws.LE.thresh
302  IF( .NOT.weak )
303  $ go to 20
304 *
305  IF( wands ) THEN
306 *
307 * Strong stability test:
308 * F-norm((A-QL**H*S*QR, B-QL**H*T*QR)) <= O(EPS*F-norm((A, B)))
309 *
310  CALL clacpy( 'Full', m, m, s, ldst, work, m )
311  CALL clacpy( 'Full', m, m, t, ldst, work( m*m+1 ), m )
312  CALL crot( 2, work, 1, work( 3 ), 1, cz, -conjg( sz ) )
313  CALL crot( 2, work( 5 ), 1, work( 7 ), 1, cz, -conjg( sz ) )
314  CALL crot( 2, work, 2, work( 2 ), 2, cq, -sq )
315  CALL crot( 2, work( 5 ), 2, work( 6 ), 2, cq, -sq )
316  DO 10 i = 1, 2
317  work( i ) = work( i ) - a( j1+i-1, j1 )
318  work( i+2 ) = work( i+2 ) - a( j1+i-1, j1+1 )
319  work( i+4 ) = work( i+4 ) - b( j1+i-1, j1 )
320  work( i+6 ) = work( i+6 ) - b( j1+i-1, j1+1 )
321  10 CONTINUE
322  scale = REAL( czero )
323  sum = REAL( cone )
324  CALL classq( 2*m*m, work, 1, scale, sum )
325  ss = scale*sqrt( sum )
326  strong = ss.LE.thresh
327  IF( .NOT.strong )
328  $ go to 20
329  END IF
330 *
331 * If the swap is accepted ("weakly" and "strongly"), apply the
332 * equivalence transformations to the original matrix pair (A,B)
333 *
334  CALL crot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, cz, conjg( sz ) )
335  CALL crot( j1+1, b( 1, j1 ), 1, b( 1, j1+1 ), 1, cz, conjg( sz ) )
336  CALL crot( n-j1+1, a( j1, j1 ), lda, a( j1+1, j1 ), lda, cq, sq )
337  CALL crot( n-j1+1, b( j1, j1 ), ldb, b( j1+1, j1 ), ldb, cq, sq )
338 *
339 * Set N1 by N2 (2,1) blocks to 0
340 *
341  a( j1+1, j1 ) = czero
342  b( j1+1, j1 ) = czero
343 *
344 * Accumulate transformations into Q and Z if requested.
345 *
346  IF( wantz )
347  $ CALL crot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, cz, conjg( sz ) )
348  IF( wantq )
349  $ CALL crot( n, q( 1, j1 ), 1, q( 1, j1+1 ), 1, cq, conjg( sq ) )
350 *
351 * Exit with INFO = 0 if swap was successfully performed.
352 *
353  RETURN
354 *
355 * Exit with INFO = 1 if swap was rejected.
356 *
357  20 CONTINUE
358  info = 1
359  RETURN
360 *
361 * End of CTGEX2
362 *
363  END