432 SUBROUTINE ztgsen( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
433 $ alpha, beta, q, ldq, z, ldz, m, pl, pr, dif,
434 $ work, lwork, iwork, liwork, info )
443 INTEGER IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK,
445 DOUBLE PRECISION PL, PR
450 DOUBLE PRECISION DIF( * )
451 COMPLEX*16 A( lda, * ), ALPHA( * ), B( ldb, * ),
452 $ beta( * ), q( ldq, * ), work( * ), z( ldz, * )
459 parameter( idifjb = 3 )
460 DOUBLE PRECISION ZERO, ONE
461 parameter( zero = 0.0d+0, one = 1.0d+0 )
464 LOGICAL LQUERY, SWAP, WANTD, WANTD1, WANTD2, WANTP
465 INTEGER I, IERR, IJB, K, KASE, KS, LIWMIN, LWMIN, MN2,
467 DOUBLE PRECISION DSCALE, DSUM, RDSCAL, SAFMIN
468 COMPLEX*16 TEMP1, TEMP2
478 INTRINSIC abs, dcmplx, dconjg, max, sqrt
481 DOUBLE PRECISION DLAMCH
489 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
491 IF( ijob.LT.0 .OR. ijob.GT.5 )
THEN
493 ELSE IF( n.LT.0 )
THEN
495 ELSE IF( lda.LT.max( 1, n ) )
THEN
497 ELSE IF( ldb.LT.max( 1, n ) )
THEN
499 ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) )
THEN
501 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
THEN
506 CALL xerbla(
'ZTGSEN', -info )
512 wantp = ijob.EQ.1 .OR. ijob.GE.4
513 wantd1 = ijob.EQ.2 .OR. ijob.EQ.4
514 wantd2 = ijob.EQ.3 .OR. ijob.EQ.5
515 wantd = wantd1 .OR. wantd2
522 alpha( k ) = a( k, k )
523 beta( k ) = b( k, k )
533 IF( ijob.EQ.1 .OR. ijob.EQ.2 .OR. ijob.EQ.4 )
THEN
534 lwmin = max( 1, 2*m*( n-m ) )
535 liwmin = max( 1, n+2 )
536 ELSE IF( ijob.EQ.3 .OR. ijob.EQ.5 )
THEN
537 lwmin = max( 1, 4*m*( n-m ) )
538 liwmin = max( 1, 2*m*( n-m ), n+2 )
547 IF( lwork.LT.lwmin .AND. .NOT.lquery )
THEN
549 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery )
THEN
554 CALL xerbla(
'ZTGSEN', -info )
556 ELSE IF( lquery )
THEN
562 IF( m.EQ.n .OR. m.EQ.0 )
THEN
571 CALL zlassq( n, a( 1, i ), 1, dscale, dsum )
572 CALL zlassq( n, b( 1, i ), 1, dscale, dsum )
574 dif( 1 ) = dscale*sqrt( dsum )
582 safmin = dlamch(
'S' )
596 $
CALL ztgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
625 CALL zlacpy(
'Full', n1, n2, a( 1, i ), lda, work, n1 )
626 CALL zlacpy(
'Full', n1, n2, b( 1, i ), ldb, work( n1*n2+1 ),
629 CALL ztgsyl(
'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
630 $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ), n1,
631 $ dscale, dif( 1 ), work( n1*n2*2+1 ),
632 $ lwork-2*n1*n2, iwork, ierr )
639 CALL zlassq( n1*n2, work, 1, rdscal, dsum )
640 pl = rdscal*sqrt( dsum )
641 IF( pl.EQ.zero )
THEN
644 pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) )
648 CALL zlassq( n1*n2, work( n1*n2+1 ), 1, rdscal, dsum )
649 pr = rdscal*sqrt( dsum )
650 IF( pr.EQ.zero )
THEN
653 pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) )
668 CALL ztgsyl(
'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
669 $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ),
670 $ n1, dscale, dif( 1 ), work( n1*n2*2+1 ),
671 $ lwork-2*n1*n2, iwork, ierr )
675 CALL ztgsyl(
'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,
676 $ n2, b( i, i ), ldb, b, ldb, work( n1*n2+1 ),
677 $ n2, dscale, dif( 2 ), work( n1*n2*2+1 ),
678 $ lwork-2*n1*n2, iwork, ierr )
696 CALL zlacn2( mn2, work( mn2+1 ), work, dif( 1 ), kase,
703 CALL ztgsyl(
'N', ijb, n1, n2, a, lda, a( i, i ), lda,
704 $ work, n1, b, ldb, b( i, i ), ldb,
705 $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
706 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
712 CALL ztgsyl(
'C', ijb, n1, n2, a, lda, a( i, i ), lda,
713 $ work, n1, b, ldb, b( i, i ), ldb,
714 $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
715 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
720 dif( 1 ) = dscale / dif( 1 )
725 CALL zlacn2( mn2, work( mn2+1 ), work, dif( 2 ), kase,
732 CALL ztgsyl(
'N', ijb, n2, n1, a( i, i ), lda, a, lda,
733 $ work, n2, b( i, i ), ldb, b, ldb,
734 $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
735 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
741 CALL ztgsyl(
'C', ijb, n2, n1, a( i, i ), lda, a, lda,
742 $ work, n2, b, ldb, b( i, i ), ldb,
743 $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
744 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
749 dif( 2 ) = dscale / dif( 2 )
758 dscale = abs( b( k, k ) )
759 IF( dscale.GT.safmin )
THEN
760 temp1 = dconjg( b( k, k ) / dscale )
761 temp2 = b( k, k ) / dscale
763 CALL zscal( n-k, temp1, b( k, k+1 ), ldb )
764 CALL zscal( n-k+1, temp1, a( k, k ), lda )
766 $
CALL zscal( n, temp2, q( 1, k ), 1 )
768 b( k, k ) = dcmplx( zero, zero )
771 alpha( k ) = a( k, k )
772 beta( k ) = b( k, k )
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 zlassq(N, X, INCX, SCALE, SUMSQ)
ZLASSQ updates a sum of squares represented in scaled form.
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine xerbla(SRNAME, INFO)
XERBLA
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 ztgexc(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, INFO)
ZTGEXC
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
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL