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
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 )
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...
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
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
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL