361 SUBROUTINE sggesx( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA,
362 $ B, LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL,
363 $ VSR, LDVSR, RCONDE, RCONDV, WORK, LWORK, IWORK,
364 $ LIWORK, BWORK, INFO )
371 CHARACTER JOBVSL, JOBVSR, SENSE, SORT
372 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N,
378 REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
379 $ b( ldb, * ), beta( * ), rconde( 2 ),
380 $ rcondv( 2 ), vsl( ldvsl, * ), vsr( ldvsr, * ),
392 PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0 )
395 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
396 $ LQUERY, LST2SL, WANTSB, WANTSE, WANTSN, WANTST,
398 INTEGER I, ICOLS, IERR, IHI, IJOB, IJOBVL, IJOBVR,
399 $ ILEFT, ILO, IP, IRIGHT, IROWS, ITAU, IWRK,
400 $ LIWMIN, LWRK, MAXWRK, MINWRK
401 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PL,
402 $ PR, SAFMAX, SAFMIN, SMLNUM
416 EXTERNAL lsame, ilaenv, slamch, slange
419 INTRINSIC abs, max, sqrt
425 IF( lsame( jobvsl,
'N' ) )
THEN
428 ELSE IF( lsame( jobvsl,
'V' ) )
THEN
436 IF( lsame( jobvsr,
'N' ) )
THEN
439 ELSE IF( lsame( jobvsr,
'V' ) )
THEN
447 wantst = lsame( sort,
'S' )
448 wantsn = lsame( sense,
'N' )
449 wantse = lsame( sense,
'E' )
450 wantsv = lsame( sense,
'V' )
451 wantsb = lsame( sense,
'B' )
452 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
455 ELSE IF( wantse )
THEN
457 ELSE IF( wantsv )
THEN
459 ELSE IF( wantsb )
THEN
466 IF( ijobvl.LE.0 )
THEN
468 ELSE IF( ijobvr.LE.0 )
THEN
470 ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort,
'N' ) ) )
THEN
472 ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsv .OR. wantsb ) .OR.
473 $ ( .NOT.wantst .AND. .NOT.wantsn ) )
THEN
475 ELSE IF( n.LT.0 )
THEN
477 ELSE IF( lda.LT.max( 1, n ) )
THEN
479 ELSE IF( ldb.LT.max( 1, n ) )
THEN
481 ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) )
THEN
483 ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) )
THEN
496 minwrk = max( 8*n, 6*n + 16 )
497 maxwrk = minwrk - n +
498 $ n*ilaenv( 1,
'SGEQRF',
' ', n, 1, n, 0 )
499 maxwrk = max( maxwrk, minwrk - n +
500 $ n*ilaenv( 1,
'SORMQR',
' ', n, 1, n, -1 ) )
502 maxwrk = max( maxwrk, minwrk - n +
503 $ n*ilaenv( 1,
'SORGQR',
' ', n, 1, n, -1 ) )
507 $ lwrk = max( lwrk, n*n/2 )
514 IF( wantsn .OR. n.EQ.0 )
THEN
521 IF( lwork.LT.minwrk .AND. .NOT.lquery )
THEN
523 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery )
THEN
529 CALL xerbla(
'SGGESX', -info )
531 ELSE IF (lquery)
THEN
545 safmin = slamch(
'S' )
546 safmax = one / safmin
547 CALL slabad( safmin, safmax )
548 smlnum = sqrt( safmin ) / eps
549 bignum = one / smlnum
553 anrm = slange(
'M', n, n, a, lda, work )
555 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
558 ELSE IF( anrm.GT.bignum )
THEN
563 $
CALL slascl(
'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
567 bnrm = slange(
'M', n, n, b, ldb, work )
569 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
572 ELSE IF( bnrm.GT.bignum )
THEN
577 $
CALL slascl(
'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
585 CALL sggbal(
'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
586 $ work( iright ), work( iwrk ), ierr )
591 irows = ihi + 1 - ilo
595 CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
596 $ work( iwrk ), lwork+1-iwrk, ierr )
601 CALL sormqr(
'L',
'T', irows, icols, irows, b( ilo, ilo ), ldb,
602 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
603 $ lwork+1-iwrk, ierr )
609 CALL slaset(
'Full', n, n, zero, one, vsl, ldvsl )
610 IF( irows.GT.1 )
THEN
611 CALL slacpy(
'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
612 $ vsl( ilo+1, ilo ), ldvsl )
614 CALL sorgqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
615 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
621 $
CALL slaset(
'Full', n, n, zero, one, vsr, ldvsr )
626 CALL sgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
627 $ ldvsl, vsr, ldvsr, ierr )
635 CALL shgeqz(
'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
636 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
637 $ work( iwrk ), lwork+1-iwrk, ierr )
639 IF( ierr.GT.0 .AND. ierr.LE.n )
THEN
641 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n )
THEN
659 CALL slascl(
'G', 0, 0, anrmto, anrm, n, 1, alphar, n,
661 CALL slascl(
'G', 0, 0, anrmto, anrm, n, 1, alphai, n,
665 $
CALL slascl(
'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
670 bwork( i ) = selctg( alphar( i ), alphai( i ), beta( i ) )
676 CALL stgsen( ijob, ilvsl, ilvsr, bwork, n, a, lda, b, ldb,
677 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
678 $ sdim, pl, pr, dif, work( iwrk ), lwork-iwrk+1,
679 $ iwork, liwork, ierr )
682 $ maxwrk = max( maxwrk, 2*sdim*( n-sdim ) )
683 IF( ierr.EQ.-22 )
THEN
689 IF( ijob.EQ.1 .OR. ijob.EQ.4 )
THEN
693 IF( ijob.EQ.2 .OR. ijob.EQ.4 )
THEN
694 rcondv( 1 ) = dif( 1 )
695 rcondv( 2 ) = dif( 2 )
707 $
CALL sggbak(
'P',
'L', n, ilo, ihi, work( ileft ),
708 $ work( iright ), n, vsl, ldvsl, ierr )
711 $
CALL sggbak(
'P',
'R', n, ilo, ihi, work( ileft ),
712 $ work( iright ), n, vsr, ldvsr, ierr )
720 IF( alphai( i ).NE.zero )
THEN
721 IF( ( alphar( i ) / safmax ).GT.( anrmto / anrm ) .OR.
722 $ ( safmin / alphar( i ) ).GT.( anrm / anrmto ) )
724 work( 1 ) = abs( a( i, i ) / alphar( i ) )
725 beta( i ) = beta( i )*work( 1 )
726 alphar( i ) = alphar( i )*work( 1 )
727 alphai( i ) = alphai( i )*work( 1 )
728 ELSE IF( ( alphai( i ) / safmax ).GT.( anrmto / anrm )
729 $ .OR. ( safmin / alphai( i ) ).GT.( anrm / anrmto ) )
731 work( 1 ) = abs( a( i, i+1 ) / alphai( i ) )
732 beta( i ) = beta( i )*work( 1 )
733 alphar( i ) = alphar( i )*work( 1 )
734 alphai( i ) = alphai( i )*work( 1 )
742 IF( alphai( i ).NE.zero )
THEN
743 IF( ( beta( i ) / safmax ).GT.( bnrmto / bnrm ) .OR.
744 $ ( safmin / beta( i ) ).GT.( bnrm / bnrmto ) )
THEN
745 work( 1 ) = abs( b( i, i ) / beta( i ) )
746 beta( i ) = beta( i )*work( 1 )
747 alphar( i ) = alphar( i )*work( 1 )
748 alphai( i ) = alphai( i )*work( 1 )
757 CALL slascl(
'H', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
758 CALL slascl(
'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
759 CALL slascl(
'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
763 CALL slascl(
'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
764 CALL slascl(
'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
776 cursl = selctg( alphar( i ), alphai( i ), beta( i ) )
777 IF( alphai( i ).EQ.zero )
THEN
781 IF( cursl .AND. .NOT.lastsl )
788 cursl = cursl .OR. lastsl
793 IF( cursl .AND. .NOT.lst2sl )
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine sggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
SGGBAK
subroutine sggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
SGGBAL
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
subroutine shgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
SHGEQZ
subroutine sggesx(JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA, B, LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, RCONDE, RCONDV, WORK, LWORK, IWORK, LIWORK, BWORK, INFO)
SGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE ...
subroutine stgsen(IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF, WORK, LWORK, IWORK, LIWORK, INFO)
STGSEN
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
subroutine sgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
SGGHRD