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.
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 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: