190 SUBROUTINE ztgex2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
200 INTEGER info, j1, lda, ldb, ldq, ldz, n
203 COMPLEX*16 a( lda, * ),
b( ldb, * ), q( ldq, * ),
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 )
216 parameter( ldst = 2 )
218 parameter( wands = .true. )
223 DOUBLE PRECISION cq, cz, eps, sa, sb, scale, smlnum, ss, sum,
225 COMPLEX*16 cdum, f, g, sq, sz
228 COMPLEX*16 s( ldst, ldst ), t( ldst, ldst ), work( 8 )
238 INTRINSIC abs, dble, dconjg, max, sqrt
255 CALL
zlacpy(
'Full', m, m, a( j1, j1 ), lda, s, ldst )
256 CALL
zlacpy(
'Full', m, m,
b( j1, j1 ), ldb, t, ldst )
261 smlnum =
dlamch(
'S' ) / eps
262 scale = dble( czero )
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 )
277 thresh = max( twenty*eps*sa, smlnum )
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 )
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 ) )
291 CALL
zlartg( s( 1, 1 ), s( 2, 1 ), cq, sq, cdum )
293 CALL
zlartg( t( 1, 1 ), t( 2, 1 ), cq, sq, cdum )
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 )
300 ws = abs( s( 2, 1 ) ) + abs( t( 2, 1 ) )
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 )
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 )
322 scale = dble( czero )
324 CALL
zlassq( 2*m*m, work, 1, scale, sum )
325 ss = scale*sqrt( sum )
326 dtrong = ss.LE.thresh
334 CALL
zrot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, cz,
336 CALL
zrot( j1+1,
b( 1, j1 ), 1,
b( 1, j1+1 ), 1, cz,
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 )
343 a( j1+1, j1 ) = czero
344 b( j1+1, j1 ) = czero
349 $ CALL
zrot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, cz,
352 $ CALL
zrot( n, q( 1, j1 ), 1, q( 1, j1+1 ), 1, cq,
subroutine zlassq(N, X, INCX, SCALE, SUMSQ)
ZLASSQ updates a sum of squares represented in scaled form.
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zlartg(F, G, CS, SN, R)
ZLARTG generates a plane rotation with real cosine and complex sine.
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 equiva...
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
double precision function dlamch(CMACH)
DLAMCH
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...