450 SUBROUTINE stgsen( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
451 $ alphar, alphai, beta, q, ldq, z, ldz, m, pl,
452 $ pr, dif, work, lwork, iwork, liwork, info )
461 INTEGER IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK,
468 REAL A( lda, * ), ALPHAI( * ), ALPHAR( * ),
469 $ b( ldb, * ), beta( * ), dif( * ), q( ldq, * ),
470 $ work( * ), z( ldz, * )
477 parameter ( idifjb = 3 )
479 parameter ( zero = 0.0e+0, one = 1.0e+0 )
482 LOGICAL LQUERY, PAIR, SWAP, WANTD, WANTD1, WANTD2,
484 INTEGER I, IERR, IJB, K, KASE, KK, KS, LIWMIN, LWMIN,
486 REAL DSCALE, DSUM, EPS, RDSCAL, SMLNUM
500 INTRINSIC max, sign, sqrt
507 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
509 IF( ijob.LT.0 .OR. ijob.GT.5 )
THEN
511 ELSE IF( n.LT.0 )
THEN
513 ELSE IF( lda.LT.max( 1, n ) )
THEN
515 ELSE IF( ldb.LT.max( 1, n ) )
THEN
517 ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) )
THEN
519 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
THEN
524 CALL xerbla(
'STGSEN', -info )
531 smlnum = slamch(
'S' ) / eps
534 wantp = ijob.EQ.1 .OR. ijob.GE.4
535 wantd1 = ijob.EQ.2 .OR. ijob.EQ.4
536 wantd2 = ijob.EQ.3 .OR. ijob.EQ.5
537 wantd = wantd1 .OR. wantd2
544 IF( .NOT.lquery .OR. ijob.NE.0 )
THEN
550 IF( a( k+1, k ).EQ.zero )
THEN
555 IF(
SELECT( k ) .OR.
SELECT( k+1 ) )
566 IF( ijob.EQ.1 .OR. ijob.EQ.2 .OR. ijob.EQ.4 )
THEN
567 lwmin = max( 1, 4*n+16, 2*m*(n-m) )
568 liwmin = max( 1, n+6 )
569 ELSE IF( ijob.EQ.3 .OR. ijob.EQ.5 )
THEN
570 lwmin = max( 1, 4*n+16, 4*m*(n-m) )
571 liwmin = max( 1, 2*m*(n-m), n+6 )
573 lwmin = max( 1, 4*n+16 )
580 IF( lwork.LT.lwmin .AND. .NOT.lquery )
THEN
582 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery )
THEN
587 CALL xerbla(
'STGSEN', -info )
589 ELSE IF( lquery )
THEN
595 IF( m.EQ.n .OR. m.EQ.0 )
THEN
604 CALL slassq( n, a( 1, i ), 1, dscale, dsum )
605 CALL slassq( n, b( 1, i ), 1, dscale, dsum )
607 dif( 1 ) = dscale*sqrt( dsum )
624 IF( a( k+1, k ).NE.zero )
THEN
626 swap = swap .OR.
SELECT( k+1 )
640 $
CALL stgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq,
641 $ z, ldz, kk, ks, work, lwork, ierr )
673 CALL slacpy(
'Full', n1, n2, a( 1, i ), lda, work, n1 )
674 CALL slacpy(
'Full', n1, n2, b( 1, i ), ldb, work( n1*n2+1 ),
676 CALL stgsyl(
'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
677 $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ), n1,
678 $ dscale, dif( 1 ), work( n1*n2*2+1 ),
679 $ lwork-2*n1*n2, iwork, ierr )
686 CALL slassq( n1*n2, work, 1, rdscal, dsum )
687 pl = rdscal*sqrt( dsum )
688 IF( pl.EQ.zero )
THEN
691 pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) )
695 CALL slassq( n1*n2, work( n1*n2+1 ), 1, rdscal, dsum )
696 pr = rdscal*sqrt( dsum )
697 IF( pr.EQ.zero )
THEN
700 pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) )
716 CALL stgsyl(
'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
717 $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ),
718 $ n1, dscale, dif( 1 ), work( 2*n1*n2+1 ),
719 $ lwork-2*n1*n2, iwork, ierr )
723 CALL stgsyl(
'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,
724 $ n2, b( i, i ), ldb, b, ldb, work( n1*n2+1 ),
725 $ n2, dscale, dif( 2 ), work( 2*n1*n2+1 ),
726 $ lwork-2*n1*n2, iwork, ierr )
745 CALL slacn2( mn2, work( mn2+1 ), work, iwork, dif( 1 ),
752 CALL stgsyl(
'N', ijb, n1, n2, a, lda, a( i, i ), lda,
753 $ work, n1, b, ldb, b( i, i ), ldb,
754 $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
755 $ work( 2*n1*n2+1 ), lwork-2*n1*n2, iwork,
761 CALL stgsyl(
'T', ijb, n1, n2, a, lda, a( i, i ), lda,
762 $ work, n1, b, ldb, b( i, i ), ldb,
763 $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
764 $ work( 2*n1*n2+1 ), lwork-2*n1*n2, iwork,
769 dif( 1 ) = dscale / dif( 1 )
774 CALL slacn2( mn2, work( mn2+1 ), work, iwork, dif( 2 ),
781 CALL stgsyl(
'N', ijb, n2, n1, a( i, i ), lda, a, lda,
782 $ work, n2, b( i, i ), ldb, b, ldb,
783 $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
784 $ work( 2*n1*n2+1 ), lwork-2*n1*n2, iwork,
790 CALL stgsyl(
'T', ijb, n2, n1, a( i, i ), lda, a, lda,
791 $ work, n2, b( i, i ), ldb, b, ldb,
792 $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
793 $ work( 2*n1*n2+1 ), lwork-2*n1*n2, iwork,
798 dif( 2 ) = dscale / dif( 2 )
815 IF( a( k+1, k ).NE.zero )
THEN
824 work( 1 ) = a( k, k )
825 work( 2 ) = a( k+1, k )
826 work( 3 ) = a( k, k+1 )
827 work( 4 ) = a( k+1, k+1 )
828 work( 5 ) = b( k, k )
829 work( 6 ) = b( k+1, k )
830 work( 7 ) = b( k, k+1 )
831 work( 8 ) = b( k+1, k+1 )
832 CALL slag2( work, 2, work( 5 ), 2, smlnum*eps, beta( k ),
833 $ beta( k+1 ), alphar( k ), alphar( k+1 ),
835 alphai( k+1 ) = -alphai( k )
839 IF( sign( one, b( k, k ) ).LT.zero )
THEN
844 a( k, i ) = -a( k, i )
845 b( k, i ) = -b( k, i )
846 IF( wantq ) q( i, k ) = -q( i, k )
850 alphar( k ) = a( k, k )
852 beta( k ) = b( k, k )
subroutine stgsen(IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF, WORK, LWORK, IWORK, LIWORK, INFO)
STGSEN
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 slacn2(N, V, X, ISGN, EST, KASE, ISAVE)
SLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
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 slassq(N, X, INCX, SCALE, SUMSQ)
SLASSQ updates a sum of squares represented in scaled form.