432 SUBROUTINE ctgsen( 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,
451 COMPLEX A( lda, * ), ALPHA( * ), B( ldb, * ),
452 $ beta( * ), q( ldq, * ), work( * ), z( ldz, * )
459 parameter( idifjb = 3 )
461 parameter( zero = 0.0e+0, one = 1.0e+0 )
464 LOGICAL LQUERY, SWAP, WANTD, WANTD1, WANTD2, WANTP
465 INTEGER I, IERR, IJB, K, KASE, KS, LIWMIN, LWMIN, MN2,
467 REAL DSCALE, DSUM, RDSCAL, SAFMIN
479 INTRINSIC abs, cmplx, conjg, max, sqrt
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(
'CTGSEN', -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
519 alpha( k ) = a( k, k )
520 beta( k ) = b( k, k )
530 IF( ijob.EQ.1 .OR. ijob.EQ.2 .OR. ijob.EQ.4 )
THEN
531 lwmin = max( 1, 2*m*(n-m) )
532 liwmin = max( 1, n+2 )
533 ELSE IF( ijob.EQ.3 .OR. ijob.EQ.5 )
THEN
534 lwmin = max( 1, 4*m*(n-m) )
535 liwmin = max( 1, 2*m*(n-m), n+2 )
544 IF( lwork.LT.lwmin .AND. .NOT.lquery )
THEN
546 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery )
THEN
551 CALL xerbla(
'CTGSEN', -info )
553 ELSE IF( lquery )
THEN
559 IF( m.EQ.n .OR. m.EQ.0 )
THEN
568 CALL classq( n, a( 1, i ), 1, dscale, dsum )
569 CALL classq( n, b( 1, i ), 1, dscale, dsum )
571 dif( 1 ) = dscale*sqrt( dsum )
579 safmin = slamch(
'S' )
593 $
CALL ctgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
622 CALL clacpy(
'Full', n1, n2, a( 1, i ), lda, work, n1 )
623 CALL clacpy(
'Full', n1, n2, b( 1, i ), ldb, work( n1*n2+1 ),
626 CALL ctgsyl(
'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
627 $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ), n1,
628 $ dscale, dif( 1 ), work( n1*n2*2+1 ),
629 $ lwork-2*n1*n2, iwork, ierr )
636 CALL classq( n1*n2, work, 1, rdscal, dsum )
637 pl = rdscal*sqrt( dsum )
638 IF( pl.EQ.zero )
THEN
641 pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) )
645 CALL classq( n1*n2, work( n1*n2+1 ), 1, rdscal, dsum )
646 pr = rdscal*sqrt( dsum )
647 IF( pr.EQ.zero )
THEN
650 pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) )
665 CALL ctgsyl(
'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
666 $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ),
667 $ n1, dscale, dif( 1 ), work( n1*n2*2+1 ),
668 $ lwork-2*n1*n2, iwork, ierr )
672 CALL ctgsyl(
'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,
673 $ n2, b( i, i ), ldb, b, ldb, work( n1*n2+1 ),
674 $ n2, dscale, dif( 2 ), work( n1*n2*2+1 ),
675 $ lwork-2*n1*n2, iwork, ierr )
693 CALL clacn2( mn2, work( mn2+1 ), work, dif( 1 ), kase,
700 CALL ctgsyl(
'N', ijb, n1, n2, a, lda, a( i, i ), lda,
701 $ work, n1, b, ldb, b( i, i ), ldb,
702 $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
703 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
709 CALL ctgsyl(
'C', ijb, n1, n2, a, lda, a( i, i ), lda,
710 $ work, n1, b, ldb, b( i, i ), ldb,
711 $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
712 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
717 dif( 1 ) = dscale / dif( 1 )
722 CALL clacn2( mn2, work( mn2+1 ), work, dif( 2 ), kase,
729 CALL ctgsyl(
'N', ijb, n2, n1, a( i, i ), lda, a, lda,
730 $ work, n2, b( i, i ), ldb, b, ldb,
731 $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
732 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
738 CALL ctgsyl(
'C', ijb, n2, n1, a( i, i ), lda, a, lda,
739 $ work, n2, b, ldb, b( i, i ), ldb,
740 $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
741 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
746 dif( 2 ) = dscale / dif( 2 )
755 dscale = abs( b( k, k ) )
756 IF( dscale.GT.safmin )
THEN
757 temp1 = conjg( b( k, k ) / dscale )
758 temp2 = b( k, k ) / dscale
760 CALL cscal( n-k, temp1, b( k, k+1 ), ldb )
761 CALL cscal( n-k+1, temp1, a( k, k ), lda )
763 $
CALL cscal( n, temp2, q( 1, k ), 1 )
765 b( k, k ) = cmplx( zero, zero )
768 alpha( k ) = a( k, k )
769 beta( k ) = b( k, k )
subroutine cscal(N, CA, CX, INCX)
CSCAL
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine ctgexc(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, INFO)
CTGEXC
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine ctgsyl(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, IWORK, INFO)
CTGSYL
subroutine classq(N, X, INCX, SCALE, SUMSQ)
CLASSQ updates a sum of squares represented in scaled form.
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
subroutine clacn2(N, V, X, EST, KASE, ISAVE)
CLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...