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
549 IF( a( k+1, k ).EQ.zero )
THEN
554 IF(
SELECT( k ) .OR.
SELECT( k+1 ) )
564 IF( ijob.EQ.1 .OR. ijob.EQ.2 .OR. ijob.EQ.4 )
THEN
565 lwmin = max( 1, 4*n+16, 2*m*(n-m) )
566 liwmin = max( 1, n+6 )
567 ELSE IF( ijob.EQ.3 .OR. ijob.EQ.5 )
THEN
568 lwmin = max( 1, 4*n+16, 4*m*(n-m) )
569 liwmin = max( 1, 2*m*(n-m), n+6 )
571 lwmin = max( 1, 4*n+16 )
578 IF( lwork.LT.lwmin .AND. .NOT.lquery )
THEN
580 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery )
THEN
585 CALL xerbla(
'STGSEN', -info )
587 ELSE IF( lquery )
THEN
593 IF( m.EQ.n .OR. m.EQ.0 )
THEN
602 CALL slassq( n, a( 1, i ), 1, dscale, dsum )
603 CALL slassq( n, b( 1, i ), 1, dscale, dsum )
605 dif( 1 ) = dscale*sqrt( dsum )
622 IF( a( k+1, k ).NE.zero )
THEN
624 swap = swap .OR.
SELECT( k+1 )
638 $
CALL stgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq,
639 $ z, ldz, kk, ks, work, lwork, ierr )
671 CALL slacpy(
'Full', n1, n2, a( 1, i ), lda, work, n1 )
672 CALL slacpy(
'Full', n1, n2, b( 1, i ), ldb, work( n1*n2+1 ),
674 CALL stgsyl(
'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
675 $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ), n1,
676 $ dscale, dif( 1 ), work( n1*n2*2+1 ),
677 $ lwork-2*n1*n2, iwork, ierr )
684 CALL slassq( n1*n2, work, 1, rdscal, dsum )
685 pl = rdscal*sqrt( dsum )
686 IF( pl.EQ.zero )
THEN
689 pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) )
693 CALL slassq( n1*n2, work( n1*n2+1 ), 1, rdscal, dsum )
694 pr = rdscal*sqrt( dsum )
695 IF( pr.EQ.zero )
THEN
698 pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) )
714 CALL stgsyl(
'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
715 $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ),
716 $ n1, dscale, dif( 1 ), work( 2*n1*n2+1 ),
717 $ lwork-2*n1*n2, iwork, ierr )
721 CALL stgsyl(
'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,
722 $ n2, b( i, i ), ldb, b, ldb, work( n1*n2+1 ),
723 $ n2, dscale, dif( 2 ), work( 2*n1*n2+1 ),
724 $ lwork-2*n1*n2, iwork, ierr )
743 CALL slacn2( mn2, work( mn2+1 ), work, iwork, dif( 1 ),
750 CALL stgsyl(
'N', ijb, n1, n2, a, lda, a( i, i ), lda,
751 $ work, n1, b, ldb, b( i, i ), ldb,
752 $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
753 $ work( 2*n1*n2+1 ), lwork-2*n1*n2, iwork,
759 CALL stgsyl(
'T', ijb, n1, n2, a, lda, a( i, i ), lda,
760 $ work, n1, b, ldb, b( i, i ), ldb,
761 $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
762 $ work( 2*n1*n2+1 ), lwork-2*n1*n2, iwork,
767 dif( 1 ) = dscale / dif( 1 )
772 CALL slacn2( mn2, work( mn2+1 ), work, iwork, dif( 2 ),
779 CALL stgsyl(
'N', ijb, n2, n1, a( i, i ), lda, a, lda,
780 $ work, n2, b( i, i ), ldb, b, ldb,
781 $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
782 $ work( 2*n1*n2+1 ), lwork-2*n1*n2, iwork,
788 CALL stgsyl(
'T', ijb, n2, n1, a( i, i ), lda, a, lda,
789 $ work, n2, b( i, i ), ldb, b, ldb,
790 $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
791 $ work( 2*n1*n2+1 ), lwork-2*n1*n2, iwork,
796 dif( 2 ) = dscale / dif( 2 )
813 IF( a( k+1, k ).NE.zero )
THEN
822 work( 1 ) = a( k, k )
823 work( 2 ) = a( k+1, k )
824 work( 3 ) = a( k, k+1 )
825 work( 4 ) = a( k+1, k+1 )
826 work( 5 ) = b( k, k )
827 work( 6 ) = b( k+1, k )
828 work( 7 ) = b( k, k+1 )
829 work( 8 ) = b( k+1, k+1 )
830 CALL slag2( work, 2, work( 5 ), 2, smlnum*eps, beta( k ),
831 $ beta( k+1 ), alphar( k ), alphar( k+1 ),
833 alphai( k+1 ) = -alphai( k )
837 IF( sign( one, b( k, k ) ).LT.zero )
THEN
842 a( k, i ) = -a( k, i )
843 b( k, i ) = -b( k, i )
844 IF( wantq ) q( i, k ) = -q( i, k )
848 alphar( k ) = a( k, k )
850 beta( k ) = b( k, k )
subroutine slassq(N, X, INCX, SCALE, SUMSQ)
SLASSQ updates a sum of squares represented in scaled form.
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 xerbla(SRNAME, INFO)
XERBLA
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
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 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