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 )
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 ) )
561 $
dnrm2( n, vr( 1, ks+1 ), 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 )