LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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 186 of file ztgex2.f.

188*
189* -- LAPACK auxiliary routine --
190* -- LAPACK is a software package provided by Univ. of Tennessee, --
191* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
192*
193* .. Scalar Arguments ..
194 LOGICAL WANTQ, WANTZ
195 INTEGER INFO, J1, LDA, LDB, LDQ, LDZ, N
196* ..
197* .. Array Arguments ..
198 COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
199 $ Z( LDZ, * )
200* ..
201*
202* =====================================================================
203*
204* .. Parameters ..
205 COMPLEX*16 CZERO, CONE
206 parameter( czero = ( 0.0d+0, 0.0d+0 ),
207 $ cone = ( 1.0d+0, 0.0d+0 ) )
208 DOUBLE PRECISION TWENTY
209 parameter( twenty = 2.0d+1 )
210 INTEGER LDST
211 parameter( ldst = 2 )
212 LOGICAL WANDS
213 parameter( wands = .true. )
214* ..
215* .. Local Scalars ..
216 LOGICAL STRONG, WEAK
217 INTEGER I, M
218 DOUBLE PRECISION CQ, CZ, EPS, SA, SB, SCALE, SMLNUM, SUM,
219 $ THRESHA, THRESHB
220 COMPLEX*16 CDUM, F, G, SQ, SZ
221* ..
222* .. Local Arrays ..
223 COMPLEX*16 S( LDST, LDST ), T( LDST, LDST ), WORK( 8 )
224* ..
225* .. External Functions ..
226 DOUBLE PRECISION DLAMCH
227 EXTERNAL dlamch
228* ..
229* .. External Subroutines ..
230 EXTERNAL zlacpy, zlartg, zlassq, zrot
231* ..
232* .. Intrinsic Functions ..
233 INTRINSIC abs, dble, dconjg, max, sqrt
234* ..
235* .. Executable Statements ..
236*
237 info = 0
238*
239* Quick return if possible
240*
241 IF( n.LE.1 )
242 $ RETURN
243*
244 m = ldst
245 weak = .false.
246 strong = .false.
247*
248* Make a local copy of selected block in (A, B)
249*
250 CALL zlacpy( 'Full', m, m, a( j1, j1 ), lda, s, ldst )
251 CALL zlacpy( 'Full', m, m, b( j1, j1 ), ldb, t, ldst )
252*
253* Compute the threshold for testing the acceptance of swapping.
254*
255 eps = dlamch( 'P' )
256 smlnum = dlamch( 'S' ) / eps
257 scale = dble( czero )
258 sum = dble( cone )
259 CALL zlacpy( 'Full', m, m, s, ldst, work, m )
260 CALL zlacpy( 'Full', m, m, t, ldst, work( m*m+1 ), m )
261 CALL zlassq( m*m, work, 1, scale, sum )
262 sa = scale*sqrt( sum )
263 scale = dble( czero )
264 sum = dble( cone )
265 CALL zlassq( m*m, work(m*m+1), 1, scale, sum )
266 sb = scale*sqrt( sum )
267*
268* THRES has been changed from
269* THRESH = MAX( TEN*EPS*SA, SMLNUM )
270* to
271* THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
272* on 04/01/10.
273* "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
274* Jim Demmel and Guillaume Revy. See forum post 1783.
275*
276 thresha = max( twenty*eps*sa, smlnum )
277 threshb = max( twenty*eps*sb, 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 ) ) * abs( t( 1, 1 ) )
285 sb = abs( s( 1, 1 ) ) * 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| <= O(EPS F-norm((A)))
299* and |T21| <= O(EPS F-norm((B)))
300*
301 weak = abs( s( 2, 1 ) ).LE.thresha .AND.
302 $ abs( t( 2, 1 ) ).LE.threshb
303 IF( .NOT.weak )
304 $ GO TO 20
305*
306 IF( wands ) THEN
307*
308* Strong stability test:
309* F-norm((A-QL**H*S*QR)) <= O(EPS*F-norm((A)))
310* and
311* F-norm((B-QL**H*T*QR)) <= O(EPS*F-norm((B)))
312*
313 CALL zlacpy( 'Full', m, m, s, ldst, work, m )
314 CALL zlacpy( 'Full', m, m, t, ldst, work( m*m+1 ), m )
315 CALL zrot( 2, work, 1, work( 3 ), 1, cz, -dconjg( sz ) )
316 CALL zrot( 2, work( 5 ), 1, work( 7 ), 1, cz,
317 $ -dconjg( sz ) )
318 CALL zrot( 2, work, 2, work( 2 ), 2, cq, -sq )
319 CALL zrot( 2, work( 5 ), 2, work( 6 ), 2, cq, -sq )
320 DO 10 i = 1, 2
321 work( i ) = work( i ) - a( j1+i-1, j1 )
322 work( i+2 ) = work( i+2 ) - a( j1+i-1, j1+1 )
323 work( i+4 ) = work( i+4 ) - b( j1+i-1, j1 )
324 work( i+6 ) = work( i+6 ) - b( j1+i-1, j1+1 )
325 10 CONTINUE
326 scale = dble( czero )
327 sum = dble( cone )
328 CALL zlassq( m*m, work, 1, scale, sum )
329 sa = scale*sqrt( sum )
330 scale = dble( czero )
331 sum = dble( cone )
332 CALL zlassq( m*m, work(m*m+1), 1, scale, sum )
333 sb = scale*sqrt( sum )
334 strong = sa.LE.thresha .AND. sb.LE.threshb
335 IF( .NOT.strong )
336 $ GO TO 20
337 END IF
338*
339* If the swap is accepted ("weakly" and "strongly"), apply the
340* equivalence transformations to the original matrix pair (A,B)
341*
342 CALL zrot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, cz,
343 $ dconjg( sz ) )
344 CALL zrot( j1+1, b( 1, j1 ), 1, b( 1, j1+1 ), 1, cz,
345 $ dconjg( sz ) )
346 CALL zrot( n-j1+1, a( j1, j1 ), lda, a( j1+1, j1 ), lda, cq,
347 $ sq )
348 CALL zrot( n-j1+1, b( j1, j1 ), ldb, b( j1+1, j1 ), ldb, cq,
349 $ sq )
350*
351* Set N1 by N2 (2,1) blocks to 0
352*
353 a( j1+1, j1 ) = czero
354 b( j1+1, j1 ) = czero
355*
356* Accumulate transformations into Q and Z if requested.
357*
358 IF( wantz )
359 $ CALL zrot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, cz,
360 $ dconjg( sz ) )
361 IF( wantq )
362 $ CALL zrot( n, q( 1, j1 ), 1, q( 1, j1+1 ), 1, cq,
363 $ dconjg( sq ) )
364*
365* Exit with INFO = 0 if swap was successfully performed.
366*
367 RETURN
368*
369* Exit with INFO = 1 if swap was rejected.
370*
371 20 CONTINUE
372 info = 1
373 RETURN
374*
375* End of ZTGEX2
376*
zlacpy
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:101
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
zlartg
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition zlartg.f90:116
zlassq
subroutine zlassq(n, x, incx, scale, sumsq)
ZLASSQ updates a sum of squares represented in scaled form.
Definition zlassq.f90:122
zrot
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:101
Here is the call graph for this function:
Here is the caller graph for this function: