LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 array, dimension (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 array, dimension (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 (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 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 complexGEauxiliary
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 ctgex2( 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 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
201  $ z( ldz, * )
202 * ..
203 *
204 * =====================================================================
205 *
206 * .. Parameters ..
207  COMPLEX CZERO, CONE
208  parameter( czero = ( 0.0e+0, 0.0e+0 ),
209  $ cone = ( 1.0e+0, 0.0e+0 ) )
210  REAL TWENTY
211  parameter( twenty = 2.0e+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  REAL CQ, CZ, EPS, SA, SB, SCALE, SMLNUM, SUM,
221  $ thresha, threshb
222  COMPLEX CDUM, F, G, SQ, SZ
223 * ..
224 * .. Local Arrays ..
225  COMPLEX S( LDST, LDST ), T( LDST, LDST ), WORK( 8 )
226 * ..
227 * .. External Functions ..
228  REAL SLAMCH
229  EXTERNAL slamch
230 * ..
231 * .. External Subroutines ..
232  EXTERNAL clacpy, clartg, classq, crot
233 * ..
234 * .. Intrinsic Functions ..
235  INTRINSIC abs, conjg, max, real, 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 clacpy( 'Full', m, m, a( j1, j1 ), lda, s, ldst )
253  CALL clacpy( 'Full', m, m, b( j1, j1 ), ldb, t, ldst )
254 *
255 * Compute the threshold for testing the acceptance of swapping.
256 *
257  eps = slamch( 'P' )
258  smlnum = slamch( 'S' ) / eps
259  scale = real( czero )
260  sum = real( cone )
261  CALL clacpy( 'Full', m, m, s, ldst, work, m )
262  CALL clacpy( 'Full', m, m, t, ldst, work( m*m+1 ), m )
263  CALL classq( m*m, work, 1, scale, sum )
264  sa = scale*sqrt( sum )
265  scale = dble( czero )
266  sum = dble( cone )
267  CALL classq( 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 clartg( g, f, cz, sz, cdum )
289  sz = -sz
290  CALL crot( 2, s( 1, 1 ), 1, s( 1, 2 ), 1, cz, conjg( sz ) )
291  CALL crot( 2, t( 1, 1 ), 1, t( 1, 2 ), 1, cz, conjg( sz ) )
292  IF( sa.GE.sb ) THEN
293  CALL clartg( s( 1, 1 ), s( 2, 1 ), cq, sq, cdum )
294  ELSE
295  CALL clartg( t( 1, 1 ), t( 2, 1 ), cq, sq, cdum )
296  END IF
297  CALL crot( 2, s( 1, 1 ), ldst, s( 2, 1 ), ldst, cq, sq )
298  CALL crot( 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, B-QL**H*T*QR)) <= O(EPS*F-norm((A, B)))
312 *
313  CALL clacpy( 'Full', m, m, s, ldst, work, m )
314  CALL clacpy( 'Full', m, m, t, ldst, work( m*m+1 ), m )
315  CALL crot( 2, work, 1, work( 3 ), 1, cz, -conjg( sz ) )
316  CALL crot( 2, work( 5 ), 1, work( 7 ), 1, cz, -conjg( sz ) )
317  CALL crot( 2, work, 2, work( 2 ), 2, cq, -sq )
318  CALL crot( 2, work( 5 ), 2, work( 6 ), 2, cq, -sq )
319  DO 10 i = 1, 2
320  work( i ) = work( i ) - a( j1+i-1, j1 )
321  work( i+2 ) = work( i+2 ) - a( j1+i-1, j1+1 )
322  work( i+4 ) = work( i+4 ) - b( j1+i-1, j1 )
323  work( i+6 ) = work( i+6 ) - b( j1+i-1, j1+1 )
324  10 CONTINUE
325  scale = dble( czero )
326  sum = dble( cone )
327  CALL classq( m*m, work, 1, scale, sum )
328  sa = scale*sqrt( sum )
329  scale = dble( czero )
330  sum = dble( cone )
331  CALL classq( m*m, work(m*m+1), 1, scale, sum )
332  sb = scale*sqrt( sum )
333  strong = sa.LE.thresha .AND. sb.LE.threshb
334  IF( .NOT.strong )
335  $ GO TO 20
336  END IF
337 *
338 * If the swap is accepted ("weakly" and "strongly"), apply the
339 * equivalence transformations to the original matrix pair (A,B)
340 *
341  CALL crot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, cz, conjg( sz ) )
342  CALL crot( j1+1, b( 1, j1 ), 1, b( 1, j1+1 ), 1, cz, conjg( sz ) )
343  CALL crot( n-j1+1, a( j1, j1 ), lda, a( j1+1, j1 ), lda, cq, sq )
344  CALL crot( n-j1+1, b( j1, j1 ), ldb, b( j1+1, j1 ), ldb, cq, sq )
345 *
346 * Set N1 by N2 (2,1) blocks to 0
347 *
348  a( j1+1, j1 ) = czero
349  b( j1+1, j1 ) = czero
350 *
351 * Accumulate transformations into Q and Z if requested.
352 *
353  IF( wantz )
354  $ CALL crot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, cz, conjg( sz ) )
355  IF( wantq )
356  $ CALL crot( n, q( 1, j1 ), 1, q( 1, j1+1 ), 1, cq, conjg( sq ) )
357 *
358 * Exit with INFO = 0 if swap was successfully performed.
359 *
360  RETURN
361 *
362 * Exit with INFO = 1 if swap was rejected.
363 *
364  20 CONTINUE
365  info = 1
366  RETURN
367 *
368 * End of CTGEX2
369 *
370  END
subroutine clartg(f, g, c, s, r)
CLARTG generates a plane rotation with real cosine and complex sine.
Definition: clartg.f90:118
subroutine classq(n, x, incx, scl, sumsq)
CLASSQ updates a sum of squares represented in scaled form.
Definition: classq.f90:126
subroutine ctgex2(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, J1, INFO)
CTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an unitary equiva...
Definition: ctgex2.f:190
subroutine crot(N, CX, INCX, CY, INCY, C, S)
CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
Definition: crot.f:103
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103