378 SUBROUTINE stgsna( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL,
379 $ LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK,
387 CHARACTER HOWMNY, JOB
388 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
393 REAL A( LDA, * ), B( LDB, * ), DIF( * ), S( * ),
394 $ vl( ldvl, * ), vr( ldvr, * ), work( * )
401 PARAMETER ( DIFDRI = 3 )
402 REAL ZERO, ONE, TWO, FOUR
403 parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0,
407 LOGICAL LQUERY, PAIR, SOMCON, WANTBH, WANTDF, WANTS
408 INTEGER I, IERR, IFST, ILST, IZ, K, KS, LWMIN, N1, N2
409 REAL ALPHAI, ALPHAR, ALPRQT, BETA, C1, C2, COND,
410 $ eps, lnrm, rnrm, root1, root2, scale, smlnum,
411 $ tmpii, tmpir, tmpri, tmprr, uhav, uhavi, uhbv,
415 REAL DUMMY( 1 ), DUMMY1( 1 )
419 REAL SDOT, SLAMCH, SLAPY2, SNRM2, SROUNDUP_LWORK
420 EXTERNAL lsame, sdot, slamch, slapy2, snrm2,
427 INTRINSIC max, min, sqrt
433 wantbh = lsame( job,
'B' )
434 wants = lsame( job,
'E' ) .OR. wantbh
435 wantdf = lsame( job,
'V' ) .OR. wantbh
437 somcon = lsame( howmny,
'S' )
440 lquery = ( lwork.EQ.-1 )
442 IF( .NOT.wants .AND. .NOT.wantdf )
THEN
444 ELSE IF( .NOT.lsame( howmny,
'A' ) .AND. .NOT.somcon )
THEN
446 ELSE IF( n.LT.0 )
THEN
448 ELSE IF( lda.LT.max( 1, n ) )
THEN
450 ELSE IF( ldb.LT.max( 1, n ) )
THEN
452 ELSE IF( wants .AND. ldvl.LT.n )
THEN
454 ELSE IF( wants .AND. ldvr.LT.n )
THEN
469 IF( a( k+1, k ).EQ.zero )
THEN
474 IF(
SELECT( k ) .OR.
SELECT( k+1 ) )
489 ELSE IF( lsame( job,
'V' ) .OR. lsame( job,
'B' ) )
THEN
490 lwmin = 2*n*( n + 2 ) + 16
494 work( 1 ) = sroundup_lwork(lwmin)
498 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery )
THEN
504 CALL xerbla(
'STGSNA', -info )
506 ELSE IF( lquery )
THEN
518 smlnum = slamch(
'S' ) / eps
531 $ pair = a( k+1, k ).NE.zero
539 IF( .NOT.
SELECT( k ) .AND. .NOT.
SELECT( k+1 ) )
542 IF( .NOT.
SELECT( k ) )
558 rnrm = slapy2( snrm2( n, vr( 1, ks ), 1 ),
559 $ snrm2( n, vr( 1, ks+1 ), 1 ) )
560 lnrm = slapy2( snrm2( n, vl( 1, ks ), 1 ),
561 $ snrm2( n, vl( 1, ks+1 ), 1 ) )
562 CALL sgemv(
'N', n, n, one, a, lda, vr( 1, ks ), 1, zero,
564 tmprr = sdot( n, work, 1, vl( 1, ks ), 1 )
565 tmpri = sdot( n, work, 1, vl( 1, ks+1 ), 1 )
566 CALL sgemv(
'N', n, n, one, a, lda, vr( 1, ks+1 ), 1,
568 tmpii = sdot( n, work, 1, vl( 1, ks+1 ), 1 )
569 tmpir = sdot( n, work, 1, vl( 1, ks ), 1 )
571 uhavi = tmpir - tmpri
572 CALL sgemv(
'N', n, n, one, b, ldb, vr( 1, ks ), 1, zero,
574 tmprr = sdot( n, work, 1, vl( 1, ks ), 1 )
575 tmpri = sdot( n, work, 1, vl( 1, ks+1 ), 1 )
576 CALL sgemv(
'N', n, n, one, b, ldb, vr( 1, ks+1 ), 1,
578 tmpii = sdot( n, work, 1, vl( 1, ks+1 ), 1 )
579 tmpir = sdot( n, work, 1, vl( 1, ks ), 1 )
581 uhbvi = tmpir - tmpri
582 uhav = slapy2( uhav, uhavi )
583 uhbv = slapy2( uhbv, uhbvi )
584 cond = slapy2( uhav, uhbv )
585 s( ks ) = cond / ( rnrm*lnrm )
592 rnrm = snrm2( n, vr( 1, ks ), 1 )
593 lnrm = snrm2( n, vl( 1, ks ), 1 )
594 CALL sgemv(
'N', n, n, one, a, lda, vr( 1, ks ), 1, zero,
596 uhav = sdot( n, work, 1, vl( 1, ks ), 1 )
597 CALL sgemv(
'N', n, n, one, b, ldb, vr( 1, ks ), 1, zero,
599 uhbv = sdot( n, work, 1, vl( 1, ks ), 1 )
600 cond = slapy2( uhav, uhbv )
601 IF( cond.EQ.zero )
THEN
604 s( ks ) = cond / ( rnrm*lnrm )
611 dif( ks ) = slapy2( a( 1, 1 ), b( 1, 1 ) )
622 work( 1 ) = a( k, k )
623 work( 2 ) = a( k+1, k )
624 work( 3 ) = a( k, k+1 )
625 work( 4 ) = a( k+1, k+1 )
626 work( 5 ) = b( k, k )
627 work( 6 ) = b( k+1, k )
628 work( 7 ) = b( k, k+1 )
629 work( 8 ) = b( k+1, k+1 )
630 CALL slag2( work, 2, work( 5 ), 2, smlnum*eps, beta,
631 $ dummy1( 1 ), alphar, dummy( 1 ), alphai )
633 c1 = two*( alphar*alphar+alphai*alphai+beta*beta )
634 c2 = four*beta*beta*alphai*alphai
635 root1 = c1 + sqrt( c1*c1-4.0*c2 )
638 cond = min( sqrt( root1 ), sqrt( root2 ) )
644 CALL slacpy(
'Full', n, n, a, lda, work, n )
645 CALL slacpy(
'Full', n, n, b, ldb, work( n*n+1 ), n )
649 CALL stgexc( .false., .false., n, work, n, work( n*n+1 ), n,
650 $ dummy, 1, dummy1, 1, ifst, ilst,
651 $ work( n*n*2+1 ), lwork-2*n*n, ierr )
667 IF( work( 2 ).NE.zero )
675 CALL stgsyl(
'N', difdri, n2, n1, work( n*n1+n1+1 ),
676 $ n, work, n, work( n1+1 ), n,
677 $ work( n*n1+n1+i ), n, work( i ), n,
678 $ work( n1+i ), n, scale, dif( ks ),
679 $ work( iz+1 ), lwork-2*n*n, iwork, ierr )
682 $ dif( ks ) = min( max( one, alprqt )*dif( ks ),
687 $ dif( ks+1 ) = dif( ks )
693 work( 1 ) = sroundup_lwork(lwmin)
subroutine xerbla(srname, info)
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
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 stgexc(wantq, wantz, n, a, lda, b, ldb, q, ldq, z, ldz, ifst, ilst, work, lwork, info)
STGEXC
subroutine stgsna(job, howmny, select, n, a, lda, b, ldb, vl, ldvl, vr, ldvr, s, dif, mm, m, work, lwork, iwork, info)
STGSNA
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