LAPACK  3.7.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
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 December 2016
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.7.0) --
194 * -- LAPACK is a software package provided by Univ. of Tennessee, --
195 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
196 * December 2016
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
subroutine classq(N, X, INCX, SCALE, SUMSQ)
CLASSQ updates a sum of squares represented in scaled form.
Definition: classq.f:108
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:105
subroutine clartg(F, G, CS, SN, R)
CLARTG generates a plane rotation with real cosine and complex sine.
Definition: clartg.f:105
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:192
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105