LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ ctgex2()

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 array, dimension (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 array, dimension (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 (LDQ,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.
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 188 of file ctgex2.f.

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 *
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 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
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: