LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine ctgex2 ( logical  WANTQ,
logical  WANTZ,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
complex, dimension( ldb, * )  B,
integer  LDB,
complex, dimension( ldq, * )  Q,
integer  LDQ,
complex, dimension( ldz, * )  Z,
integer  LDZ,
integer  J1,
integer  INFO 
)

CTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an unitary equivalence transformation.

Download CTGEX2 + dependencies [TGZ] [ZIP] [TXT]

Purpose: 
 CTGEX2 swaps adjacent diagonal 1 by 1 blocks (A11,B11) and (A22,B22)
 in an upper triangular matrix pair (A, B) by an unitary equivalence
 transformation.

 (A, B) must be in generalized Schur canonical form, that is, A and
 B are both upper triangular.

 Optionally, the matrices Q and Z of generalized Schur vectors are
 updated.

        Q(in) * A(in) * Z(in)**H = Q(out) * A(out) * Z(out)**H
        Q(in) * B(in) * Z(in)**H = Q(out) * B(out) * Z(out)**H
Parameters
[in]WANTQ
          WANTQ is LOGICAL
          .TRUE. : update the left transformation matrix Q;
          .FALSE.: do not update Q.
[in]WANTZ
          WANTZ is LOGICAL
          .TRUE. : update the right transformation matrix Z;
          .FALSE.: do not update Z.
[in]N
          N is INTEGER
          The order of the matrices A and B. N >= 0.
[in,out]A
          A is COMPLEX arrays, dimensions (LDA,N)
          On entry, the matrix A in the pair (A, B).
          On exit, the updated matrix A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A. LDA >= max(1,N).
[in,out]B
          B is COMPLEX arrays, dimensions (LDB,N)
          On entry, the matrix B in the pair (A, B).
          On exit, the updated matrix B.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= max(1,N).
[in,out]Q
          Q is COMPLEX array, dimension (LDZ,N)
          If WANTQ = .TRUE, on entry, the unitary matrix Q. On exit,
          the updated matrix Q.
          Not referenced if WANTQ = .FALSE..
[in]LDQ
          LDQ is INTEGER
          The leading dimension of the array Q. LDQ >= 1;
          If WANTQ = .TRUE., LDQ >= N.
[in,out]Z
          Z is COMPLEX array, dimension (LDZ,N)
          If WANTZ = .TRUE, on entry, the unitary matrix Z. On exit,
          the updated matrix Z.
          Not referenced if WANTZ = .FALSE..
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z. LDZ >= 1;
          If WANTZ = .TRUE., LDZ >= N.
[in]J1
          J1 is INTEGER
          The index to the first block (A11, B11).
[out]INFO
          INFO is INTEGER
           =0:  Successful exit.
           =1:  The transformed matrix pair (A, B) would be too far
                from generalized Schur form; the problem is ill-
                conditioned.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Further Details:
In the current code both weak and strong stability tests are performed. The user can omit the strong stability test by changing the internal logical parameter WANDS to .FALSE.. See ref. [2] for details.
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.
References:
[1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the Generalized Real Schur Form of a Regular Matrix Pair (A, B), in M.S. Moonen et al (eds), Linear Algebra for Large Scale and Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
[2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified Eigenvalues of a Regular Matrix Pair (A, B) and Condition Estimation: Theory, Algorithms and Software, Report UMINF-94.04, Department of Computing Science, Umea University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87. To appear in Numerical Algorithms, 1996.

Definition at line 192 of file ctgex2.f.

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 *
clartg
subroutine clartg(F, G, CS, SN, R)
CLARTG generates a plane rotation with real cosine and complex sine.
Definition: clartg.f:105
classq
subroutine classq(N, X, INCX, SCALE, SUMSQ)
CLASSQ updates a sum of squares represented in scaled form.
Definition: classq.f:108
clacpy
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
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
crot
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

Here is the call graph for this function:

Here is the caller graph for this function: