LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
ztgex2.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine ztgex2 (WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, J1, INFO)
 ZTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an unitary equivalence transformation. More...
 

Function/Subroutine Documentation

subroutine ztgex2 ( logical  WANTQ,
logical  WANTZ,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( ldq, * )  Q,
integer  LDQ,
complex*16, dimension( ldz, * )  Z,
integer  LDZ,
integer  J1,
integer  INFO 
)

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

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

Purpose:
 ZTGEX2 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*16 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*16 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*16 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*16 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 ztgex2.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*16 a( lda, * ), b( ldb, * ), q( ldq, * ),
204  $ z( ldz, * )
205 * ..
206 *
207 * =====================================================================
208 *
209 * .. Parameters ..
210  COMPLEX*16 czero, cone
211  parameter( czero = ( 0.0d+0, 0.0d+0 ),
212  $ cone = ( 1.0d+0, 0.0d+0 ) )
213  DOUBLE PRECISION twenty
214  parameter( twenty = 2.0d+1 )
215  INTEGER ldst
216  parameter( ldst = 2 )
217  LOGICAL wands
218  parameter( wands = .true. )
219 * ..
220 * .. Local Scalars ..
221  LOGICAL dtrong, weak
222  INTEGER i, m
223  DOUBLE PRECISION cq, cz, eps, sa, sb, scale, smlnum, ss, sum,
224  $ thresh, ws
225  COMPLEX*16 cdum, f, g, sq, sz
226 * ..
227 * .. Local Arrays ..
228  COMPLEX*16 s( ldst, ldst ), t( ldst, ldst ), work( 8 )
229 * ..
230 * .. External Functions ..
231  DOUBLE PRECISION dlamch
232  EXTERNAL dlamch
233 * ..
234 * .. External Subroutines ..
235  EXTERNAL zlacpy, zlartg, zlassq, zrot
236 * ..
237 * .. Intrinsic Functions ..
238  INTRINSIC abs, dble, dconjg, max, 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  dtrong = .false.
252 *
253 * Make a local copy of selected block in (A, B)
254 *
255  CALL zlacpy( 'Full', m, m, a( j1, j1 ), lda, s, ldst )
256  CALL zlacpy( 'Full', m, m, b( j1, j1 ), ldb, t, ldst )
257 *
258 * Compute the threshold for testing the acceptance of swapping.
259 *
260  eps = dlamch( 'P' )
261  smlnum = dlamch( 'S' ) / eps
262  scale = dble( czero )
263  sum = dble( cone )
264  CALL zlacpy( 'Full', m, m, s, ldst, work, m )
265  CALL zlacpy( 'Full', m, m, t, ldst, work( m*m+1 ), m )
266  CALL zlassq( 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 zlartg( g, f, cz, sz, cdum )
287  sz = -sz
288  CALL zrot( 2, s( 1, 1 ), 1, s( 1, 2 ), 1, cz, dconjg( sz ) )
289  CALL zrot( 2, t( 1, 1 ), 1, t( 1, 2 ), 1, cz, dconjg( sz ) )
290  IF( sa.GE.sb ) THEN
291  CALL zlartg( s( 1, 1 ), s( 2, 1 ), cq, sq, cdum )
292  ELSE
293  CALL zlartg( t( 1, 1 ), t( 2, 1 ), cq, sq, cdum )
294  END IF
295  CALL zrot( 2, s( 1, 1 ), ldst, s( 2, 1 ), ldst, cq, sq )
296  CALL zrot( 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 zlacpy( 'Full', m, m, s, ldst, work, m )
311  CALL zlacpy( 'Full', m, m, t, ldst, work( m*m+1 ), m )
312  CALL zrot( 2, work, 1, work( 3 ), 1, cz, -dconjg( sz ) )
313  CALL zrot( 2, work( 5 ), 1, work( 7 ), 1, cz, -dconjg( sz ) )
314  CALL zrot( 2, work, 2, work( 2 ), 2, cq, -sq )
315  CALL zrot( 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 = dble( czero )
323  sum = dble( cone )
324  CALL zlassq( 2*m*m, work, 1, scale, sum )
325  ss = scale*sqrt( sum )
326  dtrong = ss.LE.thresh
327  IF( .NOT.dtrong )
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 zrot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, cz,
335  $ dconjg( sz ) )
336  CALL zrot( j1+1, b( 1, j1 ), 1, b( 1, j1+1 ), 1, cz,
337  $ dconjg( sz ) )
338  CALL zrot( n-j1+1, a( j1, j1 ), lda, a( j1+1, j1 ), lda, cq, sq )
339  CALL zrot( n-j1+1, b( j1, j1 ), ldb, b( j1+1, j1 ), ldb, cq, sq )
340 *
341 * Set N1 by N2 (2,1) blocks to 0
342 *
343  a( j1+1, j1 ) = czero
344  b( j1+1, j1 ) = czero
345 *
346 * Accumulate transformations into Q and Z if requested.
347 *
348  IF( wantz )
349  $ CALL zrot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, cz,
350  $ dconjg( sz ) )
351  IF( wantq )
352  $ CALL zrot( n, q( 1, j1 ), 1, q( 1, j1+1 ), 1, cq,
353  $ dconjg( sq ) )
354 *
355 * Exit with INFO = 0 if swap was successfully performed.
356 *
357  RETURN
358 *
359 * Exit with INFO = 1 if swap was rejected.
360 *
361  20 CONTINUE
362  info = 1
363  RETURN
364 *
365 * End of ZTGEX2
366 *
subroutine zlassq(N, X, INCX, SCALE, SUMSQ)
ZLASSQ updates a sum of squares represented in scaled form.
Definition: zlassq.f:108
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine zlartg(F, G, CS, SN, R)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition: zlartg.f:105
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zrot(N, CX, INCX, CY, INCY, C, S)
ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors...
Definition: zrot.f:105

Here is the call graph for this function:

Here is the caller graph for this function: