221 SUBROUTINE dtgex2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
222 $ ldz, j1, n1, n2, work, lwork, info )
231 INTEGER INFO, J1, LDA, LDB, LDQ, LDZ, LWORK, N, N1, N2
234 DOUBLE PRECISION A( lda, * ), B( ldb, * ), Q( ldq, * ),
235 $ work( * ), z( ldz, * )
243 DOUBLE PRECISION ZERO, ONE
244 parameter( zero = 0.0d+0, one = 1.0d+0 )
245 DOUBLE PRECISION TWENTY
246 parameter( twenty = 2.0d+01 )
248 parameter( ldst = 4 )
250 parameter( wands = .true. )
254 INTEGER I, IDUM, LINFO, M
255 DOUBLE PRECISION BQRA21, BRQA21, DDUM, DNORM, DSCALE, DSUM, EPS,
256 $ f, g, sa, sb, scale, smlnum, ss, thresh, ws
259 INTEGER IWORK( ldst )
260 DOUBLE PRECISION AI( 2 ), AR( 2 ), BE( 2 ), IR( ldst, ldst ),
261 $ ircop( ldst, ldst ), li( ldst, ldst ),
262 $ licop( ldst, ldst ), s( ldst, ldst ),
263 $ scpy( ldst, ldst ), t( ldst, ldst ),
264 $ taul( ldst ), taur( ldst ), tcpy( ldst, ldst )
267 DOUBLE PRECISION DLAMCH
276 INTRINSIC abs, max, sqrt
284 IF( n.LE.1 .OR. n1.LE.0 .OR. n2.LE.0 )
286 IF( n1.GT.n .OR. ( j1+n1 ).GT.n )
289 IF( lwork.LT.max( 1, n*m, m*m*2 ) )
THEN
291 work( 1 ) = max( 1, n*m, m*m*2 )
300 CALL
dlaset(
'Full', ldst, ldst, zero, zero, li, ldst )
301 CALL
dlaset(
'Full', ldst, ldst, zero, zero, ir, ldst )
302 CALL
dlacpy(
'Full', m, m, a( j1, j1 ), lda, s, ldst )
303 CALL
dlacpy(
'Full', m, m, b( j1, j1 ), ldb, t, ldst )
308 smlnum = dlamch(
'S' ) / eps
311 CALL
dlacpy(
'Full', m, m, s, ldst, work, m )
312 CALL
dlassq( m*m, work, 1, dscale, dsum )
313 CALL
dlacpy(
'Full', m, m, t, ldst, work, m )
314 CALL
dlassq( m*m, work, 1, dscale, dsum )
315 dnorm = dscale*sqrt( dsum )
325 thresh = max( twenty*eps*dnorm, smlnum )
334 f = s( 2, 2 )*t( 1, 1 ) - t( 2, 2 )*s( 1, 1 )
335 g = s( 2, 2 )*t( 1, 2 ) - t( 2, 2 )*s( 1, 2 )
336 sb = abs( t( 2, 2 ) )
337 sa = abs( s( 2, 2 ) )
338 CALL
dlartg( f, g, ir( 1, 2 ), ir( 1, 1 ), ddum )
339 ir( 2, 1 ) = -ir( 1, 2 )
340 ir( 2, 2 ) = ir( 1, 1 )
341 CALL
drot( 2, s( 1, 1 ), 1, s( 1, 2 ), 1, ir( 1, 1 ),
343 CALL
drot( 2, t( 1, 1 ), 1, t( 1, 2 ), 1, ir( 1, 1 ),
346 CALL
dlartg( s( 1, 1 ), s( 2, 1 ), li( 1, 1 ), li( 2, 1 ),
349 CALL
dlartg( t( 1, 1 ), t( 2, 1 ), li( 1, 1 ), li( 2, 1 ),
352 CALL
drot( 2, s( 1, 1 ), ldst, s( 2, 1 ), ldst, li( 1, 1 ),
354 CALL
drot( 2, t( 1, 1 ), ldst, t( 2, 1 ), ldst, li( 1, 1 ),
356 li( 2, 2 ) = li( 1, 1 )
357 li( 1, 2 ) = -li( 2, 1 )
362 ws = abs( s( 2, 1 ) ) + abs( t( 2, 1 ) )
372 CALL
dlacpy(
'Full', m, m, a( j1, j1 ), lda, work( m*m+1 ),
374 CALL
dgemm(
'N',
'N', m, m, m, one, li, ldst, s, ldst, zero,
376 CALL
dgemm(
'N',
'T', m, m, m, -one, work, m, ir, ldst, one,
380 CALL
dlassq( m*m, work( m*m+1 ), 1, dscale, dsum )
382 CALL
dlacpy(
'Full', m, m, b( j1, j1 ), ldb, work( m*m+1 ),
384 CALL
dgemm(
'N',
'N', m, m, m, one, li, ldst, t, ldst, zero,
386 CALL
dgemm(
'N',
'T', m, m, m, -one, work, m, ir, ldst, one,
388 CALL
dlassq( m*m, work( m*m+1 ), 1, dscale, dsum )
389 ss = dscale*sqrt( dsum )
390 dtrong = ss.LE.thresh
398 CALL
drot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, ir( 1, 1 ),
400 CALL
drot( j1+1, b( 1, j1 ), 1, b( 1, j1+1 ), 1, ir( 1, 1 ),
402 CALL
drot( n-j1+1, a( j1, j1 ), lda, a( j1+1, j1 ), lda,
403 $ li( 1, 1 ), li( 2, 1 ) )
404 CALL
drot( n-j1+1, b( j1, j1 ), ldb, b( j1+1, j1 ), ldb,
405 $ li( 1, 1 ), li( 2, 1 ) )
415 $ CALL
drot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, ir( 1, 1 ),
418 $ CALL
drot( n, q( 1, j1 ), 1, q( 1, j1+1 ), 1, li( 1, 1 ),
435 CALL
dlacpy(
'Full', n1, n2, t( 1, n1+1 ), ldst, li, ldst )
436 CALL
dlacpy(
'Full', n1, n2, s( 1, n1+1 ), ldst,
437 $ ir( n2+1, n1+1 ), ldst )
438 CALL
dtgsy2(
'N', 0, n1, n2, s, ldst, s( n1+1, n1+1 ), ldst,
439 $ ir( n2+1, n1+1 ), ldst, t, ldst, t( n1+1, n1+1 ),
440 $ ldst, li, ldst, scale, dsum, dscale, iwork, idum,
452 CALL
dscal( n1, -one, li( 1, i ), 1 )
453 li( n1+i, i ) = scale
455 CALL
dgeqr2( m, n2, li, ldst, taul, work, linfo )
458 CALL
dorg2r( m, m, n2, li, ldst, taul, work, linfo )
469 ir( n2+i, i ) = scale
471 CALL
dgerq2( n1, m, ir( n2+1, 1 ), ldst, taur, work, linfo )
474 CALL
dorgr2( m, m, n1, ir, ldst, taur, work, linfo )
480 CALL
dgemm(
'T',
'N', m, m, m, one, li, ldst, s, ldst, zero,
482 CALL
dgemm(
'N',
'T', m, m, m, one, work, m, ir, ldst, zero, s,
484 CALL
dgemm(
'T',
'N', m, m, m, one, li, ldst, t, ldst, zero,
486 CALL
dgemm(
'N',
'T', m, m, m, one, work, m, ir, ldst, zero, t,
488 CALL
dlacpy(
'F', m, m, s, ldst, scpy, ldst )
489 CALL
dlacpy(
'F', m, m, t, ldst, tcpy, ldst )
490 CALL
dlacpy(
'F', m, m, ir, ldst, ircop, ldst )
491 CALL
dlacpy(
'F', m, m, li, ldst, licop, ldst )
496 CALL
dgerq2( m, m, t, ldst, taur, work, linfo )
499 CALL
dormr2(
'R',
'T', m, m, m, t, ldst, taur, s, ldst, work,
503 CALL
dormr2(
'L',
'N', m, m, m, t, ldst, taur, ir, ldst, work,
513 CALL
dlassq( n1, s( n2+1, i ), 1, dscale, dsum )
515 brqa21 = dscale*sqrt( dsum )
520 CALL
dgeqr2( m, m, tcpy, ldst, taul, work, linfo )
523 CALL
dorm2r(
'L',
'T', m, m, m, tcpy, ldst, taul, scpy, ldst,
525 CALL
dorm2r(
'R',
'N', m, m, m, tcpy, ldst, taul, licop, ldst,
535 CALL
dlassq( n1, scpy( n2+1, i ), 1, dscale, dsum )
537 bqra21 = dscale*sqrt( dsum )
543 IF( bqra21.LE.brqa21 .AND. bqra21.LE.thresh )
THEN
544 CALL
dlacpy(
'F', m, m, scpy, ldst, s, ldst )
545 CALL
dlacpy(
'F', m, m, tcpy, ldst, t, ldst )
546 CALL
dlacpy(
'F', m, m, ircop, ldst, ir, ldst )
547 CALL
dlacpy(
'F', m, m, licop, ldst, li, ldst )
548 ELSE IF( brqa21.GE.thresh )
THEN
554 CALL
dlaset(
'Lower', m-1, m-1, zero, zero, t(2,1), ldst )
561 CALL
dlacpy(
'Full', m, m, a( j1, j1 ), lda, work( m*m+1 ),
563 CALL
dgemm(
'N',
'N', m, m, m, one, li, ldst, s, ldst, zero,
565 CALL
dgemm(
'N',
'N', m, m, m, -one, work, m, ir, ldst, one,
569 CALL
dlassq( m*m, work( m*m+1 ), 1, dscale, dsum )
571 CALL
dlacpy(
'Full', m, m, b( j1, j1 ), ldb, work( m*m+1 ),
573 CALL
dgemm(
'N',
'N', m, m, m, one, li, ldst, t, ldst, zero,
575 CALL
dgemm(
'N',
'N', m, m, m, -one, work, m, ir, ldst, one,
577 CALL
dlassq( m*m, work( m*m+1 ), 1, dscale, dsum )
578 ss = dscale*sqrt( dsum )
579 dtrong = ( ss.LE.thresh )
588 CALL
dlaset(
'Full', n1, n2, zero, zero, s(n2+1,1), ldst )
592 CALL
dlacpy(
'F', m, m, s, ldst, a( j1, j1 ), lda )
593 CALL
dlacpy(
'F', m, m, t, ldst, b( j1, j1 ), ldb )
594 CALL
dlaset(
'Full', ldst, ldst, zero, zero, t, ldst )
603 idum = lwork - m*m - 2
605 CALL
dlagv2( a( j1, j1 ), lda, b( j1, j1 ), ldb, ar, ai, be,
606 $ work( 1 ), work( 2 ), t( 1, 1 ), t( 2, 1 ) )
607 work( m+1 ) = -work( 2 )
608 work( m+2 ) = work( 1 )
609 t( n2, n2 ) = t( 1, 1 )
610 t( 1, 2 ) = -t( 2, 1 )
616 CALL
dlagv2( a( j1+n2, j1+n2 ), lda, b( j1+n2, j1+n2 ), ldb,
617 $ taur, taul, work( m*m+1 ), work( n2*m+n2+1 ),
618 $ work( n2*m+n2+2 ), t( n2+1, n2+1 ),
620 work( m*m ) = work( n2*m+n2+1 )
621 work( m*m-1 ) = -work( n2*m+n2+2 )
622 t( m, m ) = t( n2+1, n2+1 )
623 t( m-1, m ) = -t( m, m-1 )
625 CALL
dgemm(
'T',
'N', n2, n1, n2, one, work, m, a( j1, j1+n2 ),
626 $ lda, zero, work( m*m+1 ), n2 )
627 CALL
dlacpy(
'Full', n2, n1, work( m*m+1 ), n2, a( j1, j1+n2 ),
629 CALL
dgemm(
'T',
'N', n2, n1, n2, one, work, m, b( j1, j1+n2 ),
630 $ ldb, zero, work( m*m+1 ), n2 )
631 CALL
dlacpy(
'Full', n2, n1, work( m*m+1 ), n2, b( j1, j1+n2 ),
633 CALL
dgemm(
'N',
'N', m, m, m, one, li, ldst, work, m, zero,
635 CALL
dlacpy(
'Full', m, m, work( m*m+1 ), m, li, ldst )
636 CALL
dgemm(
'N',
'N', n2, n1, n1, one, a( j1, j1+n2 ), lda,
637 $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
638 CALL
dlacpy(
'Full', n2, n1, work, n2, a( j1, j1+n2 ), lda )
639 CALL
dgemm(
'N',
'N', n2, n1, n1, one, b( j1, j1+n2 ), ldb,
640 $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
641 CALL
dlacpy(
'Full', n2, n1, work, n2, b( j1, j1+n2 ), ldb )
642 CALL
dgemm(
'T',
'N', m, m, m, one, ir, ldst, t, ldst, zero,
644 CALL
dlacpy(
'Full', m, m, work, m, ir, ldst )
649 CALL
dgemm(
'N',
'N', n, m, m, one, q( 1, j1 ), ldq, li,
650 $ ldst, zero, work, n )
651 CALL
dlacpy(
'Full', n, m, work, n, q( 1, j1 ), ldq )
656 CALL
dgemm(
'N',
'N', n, m, m, one, z( 1, j1 ), ldz, ir,
657 $ ldst, zero, work, n )
658 CALL
dlacpy(
'Full', n, m, work, n, z( 1, j1 ), ldz )
667 CALL
dgemm(
'T',
'N', m, n-i+1, m, one, li, ldst,
668 $ a( j1, i ), lda, zero, work, m )
669 CALL
dlacpy(
'Full', m, n-i+1, work, m, a( j1, i ), lda )
670 CALL
dgemm(
'T',
'N', m, n-i+1, m, one, li, ldst,
671 $ b( j1, i ), lda, zero, work, m )
672 CALL
dlacpy(
'Full', m, n-i+1, work, m, b( j1, i ), ldb )
676 CALL
dgemm(
'N',
'N', i, m, m, one, a( 1, j1 ), lda, ir,
677 $ ldst, zero, work, i )
678 CALL
dlacpy(
'Full', i, m, work, i, a( 1, j1 ), lda )
679 CALL
dgemm(
'N',
'N', i, m, m, one, b( 1, j1 ), ldb, ir,
680 $ ldst, zero, work, i )
681 CALL
dlacpy(
'Full', i, m, work, i, b( 1, j1 ), ldb )
subroutine dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dtgsy2(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, IWORK, PQ, INFO)
DTGSY2 solves the generalized Sylvester equation (unblocked algorithm).
subroutine dtgex2(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, J1, N1, N2, WORK, LWORK, INFO)
DTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an orthogonal equ...
subroutine dorg2r(M, N, K, A, LDA, TAU, WORK, INFO)
DORG2R generates all or part of the orthogonal matrix Q from a QR factorization determined by sgeqrf ...
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine dlassq(N, X, INCX, SCALE, SUMSQ)
DLASSQ updates a sum of squares represented in scaled form.
subroutine dgeqr2(M, N, A, LDA, TAU, WORK, INFO)
DGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm...
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dlagv2(A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL, CSR, SNR)
DLAGV2 computes the Generalized Schur factorization of a real 2-by-2 matrix pencil (A...
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine dorm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
DORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
subroutine dormr2(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
DORMR2 multiplies a general matrix by the orthogonal matrix from a RQ factorization determined by sge...
subroutine dgerq2(M, N, A, LDA, TAU, WORK, INFO)
DGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm...
subroutine dorgr2(M, N, K, A, LDA, TAU, WORK, INFO)
DORGR2 generates all or part of the orthogonal matrix Q from an RQ factorization determined by sgerqf...