380 SUBROUTINE stgsna( 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 REAL A( lda, * ), B( ldb, * ), DIF( * ), S( * ),
397 $ vl( ldvl, * ), vr( ldvr, * ), work( * )
404 parameter ( difdri = 3 )
405 REAL ZERO, ONE, TWO, FOUR
406 parameter ( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0,
410 LOGICAL LQUERY, PAIR, SOMCON, WANTBH, WANTDF, WANTS
411 INTEGER I, IERR, IFST, ILST, IZ, K, KS, LWMIN, N1, N2
412 REAL ALPHAI, ALPHAR, ALPRQT, BETA, C1, C2, COND,
413 $ eps, lnrm, rnrm, root1, root2, scale, smlnum,
414 $ tmpii, tmpir, tmpri, tmprr, uhav, uhavi, uhbv,
418 REAL DUMMY( 1 ), DUMMY1( 1 )
422 REAL SDOT, SLAMCH, SLAPY2, SNRM2
423 EXTERNAL lsame, sdot, slamch, slapy2, snrm2
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(
'STGSNA', -info )
508 ELSE IF( lquery )
THEN
520 smlnum = slamch(
'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 = slapy2( snrm2( n, vr( 1, ks ), 1 ),
561 $ snrm2( n, vr( 1, ks+1 ), 1 ) )
562 lnrm = slapy2( snrm2( n, vl( 1, ks ), 1 ),
563 $ snrm2( n, vl( 1, ks+1 ), 1 ) )
564 CALL sgemv(
'N', n, n, one, a, lda, vr( 1, ks ), 1, zero,
566 tmprr = sdot( n, work, 1, vl( 1, ks ), 1 )
567 tmpri = sdot( n, work, 1, vl( 1, ks+1 ), 1 )
568 CALL sgemv(
'N', n, n, one, a, lda, vr( 1, ks+1 ), 1,
570 tmpii = sdot( n, work, 1, vl( 1, ks+1 ), 1 )
571 tmpir = sdot( n, work, 1, vl( 1, ks ), 1 )
573 uhavi = tmpir - tmpri
574 CALL sgemv(
'N', n, n, one, b, ldb, vr( 1, ks ), 1, zero,
576 tmprr = sdot( n, work, 1, vl( 1, ks ), 1 )
577 tmpri = sdot( n, work, 1, vl( 1, ks+1 ), 1 )
578 CALL sgemv(
'N', n, n, one, b, ldb, vr( 1, ks+1 ), 1,
580 tmpii = sdot( n, work, 1, vl( 1, ks+1 ), 1 )
581 tmpir = sdot( n, work, 1, vl( 1, ks ), 1 )
583 uhbvi = tmpir - tmpri
584 uhav = slapy2( uhav, uhavi )
585 uhbv = slapy2( uhbv, uhbvi )
586 cond = slapy2( uhav, uhbv )
587 s( ks ) = cond / ( rnrm*lnrm )
594 rnrm = snrm2( n, vr( 1, ks ), 1 )
595 lnrm = snrm2( n, vl( 1, ks ), 1 )
596 CALL sgemv(
'N', n, n, one, a, lda, vr( 1, ks ), 1, zero,
598 uhav = sdot( n, work, 1, vl( 1, ks ), 1 )
599 CALL sgemv(
'N', n, n, one, b, ldb, vr( 1, ks ), 1, zero,
601 uhbv = sdot( n, work, 1, vl( 1, ks ), 1 )
602 cond = slapy2( uhav, uhbv )
603 IF( cond.EQ.zero )
THEN
606 s( ks ) = cond / ( rnrm*lnrm )
613 dif( ks ) = slapy2( 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 slag2( 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.0*c2 )
640 cond = min( sqrt( root1 ), sqrt( root2 ) )
646 CALL slacpy(
'Full', n, n, a, lda, work, n )
647 CALL slacpy(
'Full', n, n, b, ldb, work( n*n+1 ), n )
651 CALL stgexc( .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 stgsyl(
'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 sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
subroutine stgsyl(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, IWORK, INFO)
STGSYL
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine stgexc(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, WORK, LWORK, INFO)
STGEXC
subroutine slag2(A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2, WI)
SLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine stgsna(JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK, IWORK, INFO)
STGSNA