LAPACK 3.12.0
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 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*
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:103
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:124
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:103
Here is the call graph for this function:
Here is the caller graph for this function: