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
521 IF( .NOT.lquery .OR. ijob.NE.0 )
THEN
523 alpha( k ) = a( k, k )
524 beta( k ) = b( k, k )
535 IF( ijob.EQ.1 .OR. ijob.EQ.2 .OR. ijob.EQ.4 )
THEN
536 lwmin = max( 1, 2*m*( n-m ) )
537 liwmin = max( 1, n+2 )
538 ELSE IF( ijob.EQ.3 .OR. ijob.EQ.5 )
THEN
539 lwmin = max( 1, 4*m*( n-m ) )
540 liwmin = max( 1, 2*m*( n-m ), n+2 )
549 IF( lwork.LT.lwmin .AND. .NOT.lquery )
THEN
551 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery )
THEN
556 CALL xerbla(
'ZTGSEN', -info )
558 ELSE IF( lquery )
THEN
564 IF( m.EQ.n .OR. m.EQ.0 )
THEN
573 CALL zlassq( n, a( 1, i ), 1, dscale, dsum )
574 CALL zlassq( n, b( 1, i ), 1, dscale, dsum )
576 dif( 1 ) = dscale*sqrt( dsum )
584 safmin = dlamch(
'S' )
598 $
CALL ztgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
627 CALL zlacpy(
'Full', n1, n2, a( 1, i ), lda, work, n1 )
628 CALL zlacpy(
'Full', n1, n2, b( 1, i ), ldb, work( n1*n2+1 ),
631 CALL ztgsyl(
'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
632 $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ), n1,
633 $ dscale, dif( 1 ), work( n1*n2*2+1 ),
634 $ lwork-2*n1*n2, iwork, ierr )
641 CALL zlassq( n1*n2, work, 1, rdscal, dsum )
642 pl = rdscal*sqrt( dsum )
643 IF( pl.EQ.zero )
THEN
646 pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) )
650 CALL zlassq( n1*n2, work( n1*n2+1 ), 1, rdscal, dsum )
651 pr = rdscal*sqrt( dsum )
652 IF( pr.EQ.zero )
THEN
655 pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) )
670 CALL ztgsyl(
'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
671 $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ),
672 $ n1, dscale, dif( 1 ), work( n1*n2*2+1 ),
673 $ lwork-2*n1*n2, iwork, ierr )
677 CALL ztgsyl(
'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,
678 $ n2, b( i, i ), ldb, b, ldb, work( n1*n2+1 ),
679 $ n2, dscale, dif( 2 ), work( n1*n2*2+1 ),
680 $ lwork-2*n1*n2, iwork, ierr )
698 CALL zlacn2( mn2, work( mn2+1 ), work, dif( 1 ), kase,
705 CALL ztgsyl(
'N', ijb, n1, n2, a, lda, a( i, i ), lda,
706 $ work, n1, b, ldb, b( i, i ), ldb,
707 $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
708 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
714 CALL ztgsyl(
'C', ijb, n1, n2, a, lda, a( i, i ), lda,
715 $ work, n1, b, ldb, b( i, i ), ldb,
716 $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
717 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
722 dif( 1 ) = dscale / dif( 1 )
727 CALL zlacn2( mn2, work( mn2+1 ), work, dif( 2 ), kase,
734 CALL ztgsyl(
'N', ijb, n2, n1, a( i, i ), lda, a, lda,
735 $ work, n2, b( i, i ), ldb, b, ldb,
736 $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
737 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
743 CALL ztgsyl(
'C', ijb, n2, n1, a( i, i ), lda, a, lda,
744 $ work, n2, b, ldb, b( i, i ), ldb,
745 $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
746 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
751 dif( 2 ) = dscale / dif( 2 )
760 dscale = abs( b( k, k ) )
761 IF( dscale.GT.safmin )
THEN
762 temp1 = dconjg( b( k, k ) / dscale )
763 temp2 = b( k, k ) / dscale
765 CALL zscal( n-k, temp1, b( k, k+1 ), ldb )
766 CALL zscal( n-k+1, temp1, a( k, k ), lda )
768 $
CALL zscal( n, temp2, q( 1, k ), 1 )
770 b( k, k ) = dcmplx( zero, zero )
773 alpha( k ) = a( k, k )
774 beta( k ) = b( k, k )
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine ztgexc(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, INFO)
ZTGEXC
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 xerbla(SRNAME, INFO)
XERBLA
subroutine zlassq(N, X, INCX, SCALE, SUMSQ)
ZLASSQ updates a sum of squares represented in scaled form.
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 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