190 SUBROUTINE ctgex2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
200 INTEGER info, j1, lda, ldb, ldq, ldz, n
203 COMPLEX a( lda, * ),
b( ldb, * ), q( ldq, * ),
211 parameter( czero = ( 0.0e+0, 0.0e+0 ),
212 $ cone = ( 1.0e+0, 0.0e+0 ) )
214 parameter( twenty = 2.0e+1 )
216 parameter( ldst = 2 )
218 parameter( wands = .true. )
223 REAL cq, cz, eps, sa, sb, scale, smlnum, ss, sum,
225 COMPLEX cdum, f, g, sq, sz
228 COMPLEX s( ldst, ldst ), t( ldst, ldst ), work( 8 )
238 INTRINSIC abs, conjg, max,
REAL, sqrt
255 CALL
clacpy(
'Full', m, m, a( j1, j1 ), lda, s, ldst )
256 CALL
clacpy(
'Full', m, m,
b( j1, j1 ), ldb, t, ldst )
261 smlnum =
slamch(
'S' ) / eps
262 scale =
REAL( czero )
264 CALL
clacpy(
'Full', m, m, s, ldst, work, m )
265 CALL
clacpy(
'Full', m, m, t, ldst, work( m*m+1 ), m )
266 CALL
classq( 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
clartg( g, f, cz, sz, cdum )
288 CALL
crot( 2, s( 1, 1 ), 1, s( 1, 2 ), 1, cz, conjg( sz ) )
289 CALL
crot( 2, t( 1, 1 ), 1, t( 1, 2 ), 1, cz, conjg( sz ) )
291 CALL
clartg( s( 1, 1 ), s( 2, 1 ), cq, sq, cdum )
293 CALL
clartg( t( 1, 1 ), t( 2, 1 ), cq, sq, cdum )
295 CALL
crot( 2, s( 1, 1 ), ldst, s( 2, 1 ), ldst, cq, sq )
296 CALL
crot( 2, t( 1, 1 ), ldst, t( 2, 1 ), ldst, cq, sq )
300 ws = abs( s( 2, 1 ) ) + abs( t( 2, 1 ) )
310 CALL
clacpy(
'Full', m, m, s, ldst, work, m )
311 CALL
clacpy(
'Full', m, m, t, ldst, work( m*m+1 ), m )
312 CALL
crot( 2, work, 1, work( 3 ), 1, cz, -conjg( sz ) )
313 CALL
crot( 2, work( 5 ), 1, work( 7 ), 1, cz, -conjg( sz ) )
314 CALL
crot( 2, work, 2, work( 2 ), 2, cq, -sq )
315 CALL
crot( 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 =
REAL( czero )
324 CALL
classq( 2*m*m, work, 1, scale, sum )
325 ss = scale*sqrt( sum )
326 strong = ss.LE.thresh
334 CALL
crot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, cz, conjg( sz ) )
335 CALL
crot( j1+1,
b( 1, j1 ), 1,
b( 1, j1+1 ), 1, cz, conjg( sz ) )
336 CALL
crot( n-j1+1, a( j1, j1 ), lda, a( j1+1, j1 ), lda, cq, sq )
337 CALL
crot( n-j1+1,
b( j1, j1 ), ldb,
b( j1+1, j1 ), ldb, cq, sq )
341 a( j1+1, j1 ) = czero
342 b( j1+1, j1 ) = czero
347 $ CALL
crot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, cz, conjg( sz ) )
349 $ CALL
crot( n, q( 1, j1 ), 1, q( 1, j1+1 ), 1, cq, conjg( sq ) )
subroutine crot(N, CX, INCX, CY, INCY, C, S)
CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors...
REAL function slamch(CMACH)
SLAMCH
subroutine clartg(F, G, CS, SN, R)
CLARTG generates a plane rotation with real cosine and complex sine.
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine classq(N, X, INCX, SCALE, SUMSQ)
CLASSQ updates a sum of squares represented in scaled form.
subroutine ctgex2(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, J1, INFO)
CTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an unitary equiva...