430 SUBROUTINE ctgsen( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
431 $ ALPHA, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF,
432 $ WORK, LWORK, IWORK, LIWORK, INFO )
440 INTEGER IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK,
448 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ),
449 $ beta( * ), q( ldq, * ), work( * ), z( ldz, * )
456 PARAMETER ( IDIFJB = 3 )
458 parameter( zero = 0.0e+0, one = 1.0e+0 )
461 LOGICAL LQUERY, SWAP, WANTD, WANTD1, WANTD2, WANTP
462 INTEGER I, IERR, IJB, K, KASE, KS, LIWMIN, LWMIN, MN2,
464 REAL DSCALE, DSUM, RDSCAL, SAFMIN
472 EXTERNAL CLACN2, CLACPY, CLASSQ, CSCAL, CTGEXC, CTGSYL,
476 INTRINSIC abs, cmplx, conjg, max, sqrt
483 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
485 IF( ijob.LT.0 .OR. ijob.GT.5 )
THEN
487 ELSE IF( n.LT.0 )
THEN
489 ELSE IF( lda.LT.max( 1, n ) )
THEN
491 ELSE IF( ldb.LT.max( 1, n ) )
THEN
493 ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) )
THEN
495 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
THEN
500 CALL xerbla(
'CTGSEN', -info )
506 wantp = ijob.EQ.1 .OR. ijob.GE.4
507 wantd1 = ijob.EQ.2 .OR. ijob.EQ.4
508 wantd2 = ijob.EQ.3 .OR. ijob.EQ.5
509 wantd = wantd1 .OR. wantd2
515 IF( .NOT.lquery .OR. ijob.NE.0 )
THEN
517 alpha( k ) = a( k, k )
518 beta( k ) = b( k, k )
529 IF( ijob.EQ.1 .OR. ijob.EQ.2 .OR. ijob.EQ.4 )
THEN
530 lwmin = max( 1, 2*m*(n-m) )
531 liwmin = max( 1, n+2 )
532 ELSE IF( ijob.EQ.3 .OR. ijob.EQ.5 )
THEN
533 lwmin = max( 1, 4*m*(n-m) )
534 liwmin = max( 1, 2*m*(n-m), n+2 )
543 IF( lwork.LT.lwmin .AND. .NOT.lquery )
THEN
545 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery )
THEN
550 CALL xerbla(
'CTGSEN', -info )
552 ELSE IF( lquery )
THEN
558 IF( m.EQ.n .OR. m.EQ.0 )
THEN
567 CALL classq( n, a( 1, i ), 1, dscale, dsum )
568 CALL classq( n, b( 1, i ), 1, dscale, dsum )
570 dif( 1 ) = dscale*sqrt( dsum )
578 safmin = slamch(
'S' )
592 $
CALL ctgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
621 CALL clacpy(
'Full', n1, n2, a( 1, i ), lda, work, n1 )
622 CALL clacpy(
'Full', n1, n2, b( 1, i ), ldb, work( n1*n2+1 ),
625 CALL ctgsyl(
'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
626 $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ), n1,
627 $ dscale, dif( 1 ), work( n1*n2*2+1 ),
628 $ lwork-2*n1*n2, iwork, ierr )
635 CALL classq( n1*n2, work, 1, rdscal, dsum )
636 pl = rdscal*sqrt( dsum )
637 IF( pl.EQ.zero )
THEN
640 pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) )
644 CALL classq( n1*n2, work( n1*n2+1 ), 1, rdscal, dsum )
645 pr = rdscal*sqrt( dsum )
646 IF( pr.EQ.zero )
THEN
649 pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) )
664 CALL ctgsyl(
'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
665 $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ),
666 $ n1, dscale, dif( 1 ), work( n1*n2*2+1 ),
667 $ lwork-2*n1*n2, iwork, ierr )
671 CALL ctgsyl(
'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,
672 $ n2, b( i, i ), ldb, b, ldb, work( n1*n2+1 ),
673 $ n2, dscale, dif( 2 ), work( n1*n2*2+1 ),
674 $ lwork-2*n1*n2, iwork, ierr )
692 CALL clacn2( mn2, work( mn2+1 ), work, dif( 1 ), kase,
699 CALL ctgsyl(
'N', ijb, n1, n2, a, lda, a( i, i ), lda,
700 $ work, n1, b, ldb, b( i, i ), ldb,
701 $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
702 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
708 CALL ctgsyl(
'C', ijb, n1, n2, a, lda, a( i, i ), lda,
709 $ work, n1, b, ldb, b( i, i ), ldb,
710 $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
711 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
716 dif( 1 ) = dscale / dif( 1 )
721 CALL clacn2( mn2, work( mn2+1 ), work, dif( 2 ), kase,
728 CALL ctgsyl(
'N', ijb, n2, n1, a( i, i ), lda, a, lda,
729 $ work, n2, b( i, i ), ldb, b, ldb,
730 $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
731 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
737 CALL ctgsyl(
'C', ijb, n2, n1, a( i, i ), lda, a, lda,
738 $ work, n2, b, ldb, b( i, i ), ldb,
739 $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
740 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
745 dif( 2 ) = dscale / dif( 2 )
754 dscale = abs( b( k, k ) )
755 IF( dscale.GT.safmin )
THEN
756 temp1 = conjg( b( k, k ) / dscale )
757 temp2 = b( k, k ) / dscale
759 CALL cscal( n-k, temp1, b( k, k+1 ), ldb )
760 CALL cscal( n-k+1, temp1, a( k, k ), lda )
762 $
CALL cscal( n, temp2, q( 1, k ), 1 )
764 b( k, k ) = cmplx( zero, zero )
767 alpha( k ) = a( k, k )
768 beta( k ) = b( k, k )
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine ctgsen(IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB, ALPHA, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF, WORK, LWORK, IWORK, LIWORK, INFO)
CTGSEN