380 SUBROUTINE dtgsna( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL,
381 $ ldvl, vr, ldvr, s, dif, mm, m, work, lwork,
390 CHARACTER HOWMNY, JOB
391 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
396 DOUBLE PRECISION A( lda, * ), B( ldb, * ), DIF( * ), S( * ),
397 $ vl( ldvl, * ), vr( ldvr, * ), work( * )
404 parameter ( difdri = 3 )
405 DOUBLE PRECISION ZERO, ONE, TWO, FOUR
406 parameter ( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0,
410 LOGICAL LQUERY, PAIR, SOMCON, WANTBH, WANTDF, WANTS
411 INTEGER I, IERR, IFST, ILST, IZ, K, KS, LWMIN, N1, N2
412 DOUBLE PRECISION ALPHAI, ALPHAR, ALPRQT, BETA, C1, C2, COND,
413 $ eps, lnrm, rnrm, root1, root2, scale, smlnum,
414 $ tmpii, tmpir, tmpri, tmprr, uhav, uhavi, uhbv,
418 DOUBLE PRECISION DUMMY( 1 ), DUMMY1( 1 )
422 DOUBLE PRECISION DDOT, DLAMCH, DLAPY2, DNRM2
423 EXTERNAL lsame, ddot, dlamch, dlapy2, dnrm2
429 INTRINSIC max, min, sqrt
435 wantbh = lsame( job,
'B' )
436 wants = lsame( job,
'E' ) .OR. wantbh
437 wantdf = lsame( job,
'V' ) .OR. wantbh
439 somcon = lsame( howmny,
'S' )
442 lquery = ( lwork.EQ.-1 )
444 IF( .NOT.wants .AND. .NOT.wantdf )
THEN
446 ELSE IF( .NOT.lsame( howmny,
'A' ) .AND. .NOT.somcon )
THEN
448 ELSE IF( n.LT.0 )
THEN
450 ELSE IF( lda.LT.max( 1, n ) )
THEN
452 ELSE IF( ldb.LT.max( 1, n ) )
THEN
454 ELSE IF( wants .AND. ldvl.LT.n )
THEN
456 ELSE IF( wants .AND. ldvr.LT.n )
THEN
471 IF( a( k+1, k ).EQ.zero )
THEN
476 IF(
SELECT( k ) .OR.
SELECT( k+1 ) )
491 ELSE IF( lsame( job,
'V' ) .OR. lsame( job,
'B' ) )
THEN
492 lwmin = 2*n*( n + 2 ) + 16
500 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery )
THEN
506 CALL xerbla(
'DTGSNA', -info )
508 ELSE IF( lquery )
THEN
520 smlnum = dlamch(
'S' ) / eps
533 $ pair = a( k+1, k ).NE.zero
541 IF( .NOT.
SELECT( k ) .AND. .NOT.
SELECT( k+1 ) )
544 IF( .NOT.
SELECT( k ) )
560 rnrm = dlapy2( dnrm2( n, vr( 1, ks ), 1 ),
561 $ dnrm2( n, vr( 1, ks+1 ), 1 ) )
562 lnrm = dlapy2( dnrm2( n, vl( 1, ks ), 1 ),
563 $ dnrm2( n, vl( 1, ks+1 ), 1 ) )
564 CALL dgemv(
'N', n, n, one, a, lda, vr( 1, ks ), 1, zero,
566 tmprr = ddot( n, work, 1, vl( 1, ks ), 1 )
567 tmpri = ddot( n, work, 1, vl( 1, ks+1 ), 1 )
568 CALL dgemv(
'N', n, n, one, a, lda, vr( 1, ks+1 ), 1,
570 tmpii = ddot( n, work, 1, vl( 1, ks+1 ), 1 )
571 tmpir = ddot( n, work, 1, vl( 1, ks ), 1 )
573 uhavi = tmpir - tmpri
574 CALL dgemv(
'N', n, n, one, b, ldb, vr( 1, ks ), 1, zero,
576 tmprr = ddot( n, work, 1, vl( 1, ks ), 1 )
577 tmpri = ddot( n, work, 1, vl( 1, ks+1 ), 1 )
578 CALL dgemv(
'N', n, n, one, b, ldb, vr( 1, ks+1 ), 1,
580 tmpii = ddot( n, work, 1, vl( 1, ks+1 ), 1 )
581 tmpir = ddot( n, work, 1, vl( 1, ks ), 1 )
583 uhbvi = tmpir - tmpri
584 uhav = dlapy2( uhav, uhavi )
585 uhbv = dlapy2( uhbv, uhbvi )
586 cond = dlapy2( uhav, uhbv )
587 s( ks ) = cond / ( rnrm*lnrm )
594 rnrm = dnrm2( n, vr( 1, ks ), 1 )
595 lnrm = dnrm2( n, vl( 1, ks ), 1 )
596 CALL dgemv(
'N', n, n, one, a, lda, vr( 1, ks ), 1, zero,
598 uhav = ddot( n, work, 1, vl( 1, ks ), 1 )
599 CALL dgemv(
'N', n, n, one, b, ldb, vr( 1, ks ), 1, zero,
601 uhbv = ddot( n, work, 1, vl( 1, ks ), 1 )
602 cond = dlapy2( uhav, uhbv )
603 IF( cond.EQ.zero )
THEN
606 s( ks ) = cond / ( rnrm*lnrm )
613 dif( ks ) = dlapy2( a( 1, 1 ), b( 1, 1 ) )
624 work( 1 ) = a( k, k )
625 work( 2 ) = a( k+1, k )
626 work( 3 ) = a( k, k+1 )
627 work( 4 ) = a( k+1, k+1 )
628 work( 5 ) = b( k, k )
629 work( 6 ) = b( k+1, k )
630 work( 7 ) = b( k, k+1 )
631 work( 8 ) = b( k+1, k+1 )
632 CALL dlag2( work, 2, work( 5 ), 2, smlnum*eps, beta,
633 $ dummy1( 1 ), alphar, dummy( 1 ), alphai )
635 c1 = two*( alphar*alphar+alphai*alphai+beta*beta )
636 c2 = four*beta*beta*alphai*alphai
637 root1 = c1 + sqrt( c1*c1-4.0d0*c2 )
640 cond = min( sqrt( root1 ), sqrt( root2 ) )
646 CALL dlacpy(
'Full', n, n, a, lda, work, n )
647 CALL dlacpy(
'Full', n, n, b, ldb, work( n*n+1 ), n )
651 CALL dtgexc( .false., .false., n, work, n, work( n*n+1 ), n,
652 $ dummy, 1, dummy1, 1, ifst, ilst,
653 $ work( n*n*2+1 ), lwork-2*n*n, ierr )
669 IF( work( 2 ).NE.zero )
677 CALL dtgsyl(
'N', difdri, n2, n1, work( n*n1+n1+1 ),
678 $ n, work, n, work( n1+1 ), n,
679 $ work( n*n1+n1+i ), n, work( i ), n,
680 $ work( n1+i ), n, scale, dif( ks ),
681 $ work( iz+1 ), lwork-2*n*n, iwork, ierr )
684 $ dif( ks ) = min( max( one, alprqt )*dif( ks ),
689 $ dif( ks+1 ) = dif( ks )
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
subroutine dtgexc(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, WORK, LWORK, INFO)
DTGEXC
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlag2(A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2, WI)
DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
subroutine dtgsyl(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, IWORK, INFO)
DTGSYL
subroutine dtgsna(JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK, IWORK, INFO)
DTGSNA