 LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ ztgex2()

 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 array, 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 array, 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 (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*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.```
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.  for details.
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.
References:
 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.
 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 ztgex2.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*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
201  \$ Z( LDZ, * )
202 * ..
203 *
204 * =====================================================================
205 *
206 * .. Parameters ..
207  COMPLEX*16 CZERO, CONE
208  parameter( czero = ( 0.0d+0, 0.0d+0 ),
209  \$ cone = ( 1.0d+0, 0.0d+0 ) )
210  DOUBLE PRECISION TWENTY
211  parameter( twenty = 2.0d+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  DOUBLE PRECISION CQ, CZ, EPS, SA, SB, SCALE, SMLNUM, SUM,
221  \$ THRESHA, THRESHB
222  COMPLEX*16 CDUM, F, G, SQ, SZ
223 * ..
224 * .. Local Arrays ..
225  COMPLEX*16 S( LDST, LDST ), T( LDST, LDST ), WORK( 8 )
226 * ..
227 * .. External Functions ..
228  DOUBLE PRECISION DLAMCH
229  EXTERNAL dlamch
230 * ..
231 * .. External Subroutines ..
232  EXTERNAL zlacpy, zlartg, zlassq, zrot
233 * ..
234 * .. Intrinsic Functions ..
235  INTRINSIC abs, dble, dconjg, max, 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 zlacpy( 'Full', m, m, a( j1, j1 ), lda, s, ldst )
253  CALL zlacpy( 'Full', m, m, b( j1, j1 ), ldb, t, ldst )
254 *
255 * Compute the threshold for testing the acceptance of swapping.
256 *
257  eps = dlamch( 'P' )
258  smlnum = dlamch( 'S' ) / eps
259  scale = dble( czero )
260  sum = dble( cone )
261  CALL zlacpy( 'Full', m, m, s, ldst, work, m )
262  CALL zlacpy( 'Full', m, m, t, ldst, work( m*m+1 ), m )
263  CALL zlassq( m*m, work, 1, scale, sum )
264  sa = scale*sqrt( sum )
265  scale = dble( czero )
266  sum = dble( cone )
267  CALL zlassq( 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 zlartg( g, f, cz, sz, cdum )
289  sz = -sz
290  CALL zrot( 2, s( 1, 1 ), 1, s( 1, 2 ), 1, cz, dconjg( sz ) )
291  CALL zrot( 2, t( 1, 1 ), 1, t( 1, 2 ), 1, cz, dconjg( sz ) )
292  IF( sa.GE.sb ) THEN
293  CALL zlartg( s( 1, 1 ), s( 2, 1 ), cq, sq, cdum )
294  ELSE
295  CALL zlartg( t( 1, 1 ), t( 2, 1 ), cq, sq, cdum )
296  END IF
297  CALL zrot( 2, s( 1, 1 ), ldst, s( 2, 1 ), ldst, cq, sq )
298  CALL zrot( 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)) <= O(EPS*F-norm((A)))
312 * and
313 * F-norm((B-QL**H*T*QR)) <= O(EPS*F-norm((B)))
314 *
315  CALL zlacpy( 'Full', m, m, s, ldst, work, m )
316  CALL zlacpy( 'Full', m, m, t, ldst, work( m*m+1 ), m )
317  CALL zrot( 2, work, 1, work( 3 ), 1, cz, -dconjg( sz ) )
318  CALL zrot( 2, work( 5 ), 1, work( 7 ), 1, cz, -dconjg( sz ) )
319  CALL zrot( 2, work, 2, work( 2 ), 2, cq, -sq )
320  CALL zrot( 2, work( 5 ), 2, work( 6 ), 2, cq, -sq )
321  DO 10 i = 1, 2
322  work( i ) = work( i ) - a( j1+i-1, j1 )
323  work( i+2 ) = work( i+2 ) - a( j1+i-1, j1+1 )
324  work( i+4 ) = work( i+4 ) - b( j1+i-1, j1 )
325  work( i+6 ) = work( i+6 ) - b( j1+i-1, j1+1 )
326  10 CONTINUE
327  scale = dble( czero )
328  sum = dble( cone )
329  CALL zlassq( m*m, work, 1, scale, sum )
330  sa = scale*sqrt( sum )
331  scale = dble( czero )
332  sum = dble( cone )
333  CALL zlassq( m*m, work(m*m+1), 1, scale, sum )
334  sb = scale*sqrt( sum )
335  strong = sa.LE.thresha .AND. sb.LE.threshb
336  IF( .NOT.strong )
337  \$ GO TO 20
338  END IF
339 *
340 * If the swap is accepted ("weakly" and "strongly"), apply the
341 * equivalence transformations to the original matrix pair (A,B)
342 *
343  CALL zrot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, cz,
344  \$ dconjg( sz ) )
345  CALL zrot( j1+1, b( 1, j1 ), 1, b( 1, j1+1 ), 1, cz,
346  \$ dconjg( sz ) )
347  CALL zrot( n-j1+1, a( j1, j1 ), lda, a( j1+1, j1 ), lda, cq, sq )
348  CALL zrot( n-j1+1, b( j1, j1 ), ldb, b( j1+1, j1 ), ldb, cq, sq )
349 *
350 * Set N1 by N2 (2,1) blocks to 0
351 *
352  a( j1+1, j1 ) = czero
353  b( j1+1, j1 ) = czero
354 *
355 * Accumulate transformations into Q and Z if requested.
356 *
357  IF( wantz )
358  \$ CALL zrot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, cz,
359  \$ dconjg( sz ) )
360  IF( wantq )
361  \$ CALL zrot( n, q( 1, j1 ), 1, q( 1, j1+1 ), 1, cq,
362  \$ dconjg( sq ) )
363 *
364 * Exit with INFO = 0 if swap was successfully performed.
365 *
366  RETURN
367 *
368 * Exit with INFO = 1 if swap was rejected.
369 *
370  20 CONTINUE
371  info = 1
372  RETURN
373 *
374 * End of ZTGEX2
375 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition: zlartg.f90:118
subroutine zlassq(n, x, incx, scl, sumsq)
ZLASSQ updates a sum of squares represented in scaled form.
Definition: zlassq.f90:137
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:103
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
Here is the call graph for this function:
Here is the caller graph for this function: