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 )
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 ) )
561 $
snrm2( n, vr( 1, ks+1 ), 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 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 ...
real function sdot(N, SX, INCX, SY, INCY)
SDOT
subroutine xerbla(SRNAME, INFO)
XERBLA
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
logical function lsame(CA, CB)
LSAME
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.
real function slamch(CMACH)
SLAMCH
subroutine stgexc(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, WORK, LWORK, INFO)
STGEXC
real function snrm2(N, X, INCX)
SNRM2
subroutine stgsna(JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK, IWORK, INFO)
STGSNA
real function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).
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