430 SUBROUTINE ztgsen( 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,
442 DOUBLE PRECISION PL, PR
447 DOUBLE PRECISION DIF( * )
448 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
449 $ beta( * ), q( ldq, * ), work( * ), z( ldz, * )
456 PARAMETER ( IDIFJB = 3 )
457 DOUBLE PRECISION ZERO, ONE
458 parameter( zero = 0.0d+0, one = 1.0d+0 )
461 LOGICAL LQUERY, SWAP, WANTD, WANTD1, WANTD2, WANTP
462 INTEGER I, IERR, IJB, K, KASE, KS, LIWMIN, LWMIN, MN2,
464 DOUBLE PRECISION DSCALE, DSUM, RDSCAL, SAFMIN
465 COMPLEX*16 TEMP1, TEMP2
475 INTRINSIC abs, dcmplx, dconjg, max, sqrt
478 DOUBLE PRECISION DLAMCH
486 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
488 IF( ijob.LT.0 .OR. ijob.GT.5 )
THEN
490 ELSE IF( n.LT.0 )
THEN
492 ELSE IF( lda.LT.max( 1, n ) )
THEN
494 ELSE IF( ldb.LT.max( 1, n ) )
THEN
496 ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) )
THEN
498 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
THEN
503 CALL xerbla(
'ZTGSEN', -info )
509 wantp = ijob.EQ.1 .OR. ijob.GE.4
510 wantd1 = ijob.EQ.2 .OR. ijob.EQ.4
511 wantd2 = ijob.EQ.3 .OR. ijob.EQ.5
512 wantd = wantd1 .OR. wantd2
518 IF( .NOT.lquery .OR. ijob.NE.0 )
THEN
520 alpha( k ) = a( k, k )
521 beta( k ) = b( k, k )
532 IF( ijob.EQ.1 .OR. ijob.EQ.2 .OR. ijob.EQ.4 )
THEN
533 lwmin = max( 1, 2*m*( n-m ) )
534 liwmin = max( 1, n+2 )
535 ELSE IF( ijob.EQ.3 .OR. ijob.EQ.5 )
THEN
536 lwmin = max( 1, 4*m*( n-m ) )
537 liwmin = max( 1, 2*m*( n-m ), n+2 )
546 IF( lwork.LT.lwmin .AND. .NOT.lquery )
THEN
548 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery )
THEN
553 CALL xerbla(
'ZTGSEN', -info )
555 ELSE IF( lquery )
THEN
561 IF( m.EQ.n .OR. m.EQ.0 )
THEN
570 CALL zlassq( n, a( 1, i ), 1, dscale, dsum )
571 CALL zlassq( n, b( 1, i ), 1, dscale, dsum )
573 dif( 1 ) = dscale*sqrt( dsum )
581 safmin = dlamch(
'S' )
595 $
CALL ztgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
624 CALL zlacpy(
'Full', n1, n2, a( 1, i ), lda, work, n1 )
625 CALL zlacpy(
'Full', n1, n2, b( 1, i ), ldb, work( n1*n2+1 ),
628 CALL ztgsyl(
'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
629 $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ), n1,
630 $ dscale, dif( 1 ), work( n1*n2*2+1 ),
631 $ lwork-2*n1*n2, iwork, ierr )
638 CALL zlassq( n1*n2, work, 1, rdscal, dsum )
639 pl = rdscal*sqrt( dsum )
640 IF( pl.EQ.zero )
THEN
643 pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) )
647 CALL zlassq( n1*n2, work( n1*n2+1 ), 1, rdscal, dsum )
648 pr = rdscal*sqrt( dsum )
649 IF( pr.EQ.zero )
THEN
652 pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) )
667 CALL ztgsyl(
'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
668 $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ),
669 $ n1, dscale, dif( 1 ), work( n1*n2*2+1 ),
670 $ lwork-2*n1*n2, iwork, ierr )
674 CALL ztgsyl(
'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,
675 $ n2, b( i, i ), ldb, b, ldb, work( n1*n2+1 ),
676 $ n2, dscale, dif( 2 ), work( n1*n2*2+1 ),
677 $ lwork-2*n1*n2, iwork, ierr )
695 CALL zlacn2( mn2, work( mn2+1 ), work, dif( 1 ), kase,
702 CALL ztgsyl(
'N', ijb, n1, n2, a, lda, a( i, i ), lda,
703 $ work, n1, b, ldb, b( i, i ), ldb,
704 $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
705 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
711 CALL ztgsyl(
'C', ijb, n1, n2, a, lda, a( i, i ), lda,
712 $ work, n1, b, ldb, b( i, i ), ldb,
713 $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
714 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
719 dif( 1 ) = dscale / dif( 1 )
724 CALL zlacn2( mn2, work( mn2+1 ), work, dif( 2 ), kase,
731 CALL ztgsyl(
'N', ijb, n2, n1, a( i, i ), lda, a, lda,
732 $ work, n2, b( i, i ), ldb, b, ldb,
733 $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
734 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
740 CALL ztgsyl(
'C', ijb, n2, n1, a( i, i ), lda, a, lda,
741 $ work, n2, b, ldb, b( i, i ), ldb,
742 $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
743 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
748 dif( 2 ) = dscale / dif( 2 )
757 dscale = abs( b( k, k ) )
758 IF( dscale.GT.safmin )
THEN
759 temp1 = dconjg( b( k, k ) / dscale )
760 temp2 = b( k, k ) / dscale
762 CALL zscal( n-k, temp1, b( k, k+1 ), ldb )
763 CALL zscal( n-k+1, temp1, a( k, k ), lda )
765 $
CALL zscal( n, temp2, q( 1, k ), 1 )
767 b( k, k ) = dcmplx( zero, zero )
770 alpha( k ) = a( k, k )
771 beta( k ) = b( k, k )
subroutine zlassq(n, x, incx, scl, sumsq)
ZLASSQ updates a sum of squares represented in scaled form.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
subroutine ztgexc(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, INFO)
ZTGEXC
subroutine zlacn2(N, V, X, EST, KASE, ISAVE)
ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine ztgsen(IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB, ALPHA, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF, WORK, LWORK, IWORK, LIWORK, INFO)
ZTGSEN
subroutine ztgsyl(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, IWORK, INFO)
ZTGSYL