363 SUBROUTINE sggesx( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA,
364 $ b, ldb, sdim, alphar, alphai, beta, vsl, ldvsl,
365 $ vsr, ldvsr, rconde, rcondv, work, lwork, iwork,
366 $ liwork, bwork, info )
374 CHARACTER JOBVSL, JOBVSR, SENSE, SORT
375 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N,
381 REAL A( lda, * ), ALPHAI( * ), ALPHAR( * ),
382 $ b( ldb, * ), beta( * ), rconde( 2 ),
383 $ rcondv( 2 ), vsl( ldvsl, * ), vsr( ldvsr, * ),
395 parameter ( zero = 0.0e+0, one = 1.0e+0 )
398 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
399 $ lquery, lst2sl, wantsb, wantse, wantsn, wantst,
401 INTEGER I, ICOLS, IERR, IHI, IJOB, IJOBVL, IJOBVR,
402 $ ileft, ilo, ip, iright, irows, itau, iwrk,
403 $ liwmin, lwrk, maxwrk, minwrk
404 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PL,
405 $ pr, safmax, safmin, smlnum
419 EXTERNAL lsame, ilaenv, slamch, slange
422 INTRINSIC abs, max, sqrt
428 IF( lsame( jobvsl,
'N' ) )
THEN
431 ELSE IF( lsame( jobvsl,
'V' ) )
THEN
439 IF( lsame( jobvsr,
'N' ) )
THEN
442 ELSE IF( lsame( jobvsr,
'V' ) )
THEN
450 wantst = lsame( sort,
'S' )
451 wantsn = lsame( sense,
'N' )
452 wantse = lsame( sense,
'E' )
453 wantsv = lsame( sense,
'V' )
454 wantsb = lsame( sense,
'B' )
455 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
458 ELSE IF( wantse )
THEN
460 ELSE IF( wantsv )
THEN
462 ELSE IF( wantsb )
THEN
469 IF( ijobvl.LE.0 )
THEN
471 ELSE IF( ijobvr.LE.0 )
THEN
473 ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort,
'N' ) ) )
THEN
475 ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsv .OR. wantsb ) .OR.
476 $ ( .NOT.wantst .AND. .NOT.wantsn ) )
THEN
478 ELSE IF( n.LT.0 )
THEN
480 ELSE IF( lda.LT.max( 1, n ) )
THEN
482 ELSE IF( ldb.LT.max( 1, n ) )
THEN
484 ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) )
THEN
486 ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) )
THEN
499 minwrk = max( 8*n, 6*n + 16 )
500 maxwrk = minwrk - n +
501 $ n*ilaenv( 1,
'SGEQRF',
' ', n, 1, n, 0 )
502 maxwrk = max( maxwrk, minwrk - n +
503 $ n*ilaenv( 1,
'SORMQR',
' ', n, 1, n, -1 ) )
505 maxwrk = max( maxwrk, minwrk - n +
506 $ n*ilaenv( 1,
'SORGQR',
' ', n, 1, n, -1 ) )
510 $ lwrk = max( lwrk, n*n/2 )
517 IF( wantsn .OR. n.EQ.0 )
THEN
524 IF( lwork.LT.minwrk .AND. .NOT.lquery )
THEN
526 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery )
THEN
532 CALL xerbla(
'SGGESX', -info )
534 ELSE IF (lquery)
THEN
548 safmin = slamch(
'S' )
549 safmax = one / safmin
550 CALL slabad( safmin, safmax )
551 smlnum = sqrt( safmin ) / eps
552 bignum = one / smlnum
556 anrm = slange(
'M', n, n, a, lda, work )
558 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
561 ELSE IF( anrm.GT.bignum )
THEN
566 $
CALL slascl(
'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
570 bnrm = slange(
'M', n, n, b, ldb, work )
572 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
575 ELSE IF( bnrm.GT.bignum )
THEN
580 $
CALL slascl(
'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
588 CALL sggbal(
'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
589 $ work( iright ), work( iwrk ), ierr )
594 irows = ihi + 1 - ilo
598 CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
599 $ work( iwrk ), lwork+1-iwrk, ierr )
604 CALL sormqr(
'L',
'T', irows, icols, irows, b( ilo, ilo ), ldb,
605 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
606 $ lwork+1-iwrk, ierr )
612 CALL slaset(
'Full', n, n, zero, one, vsl, ldvsl )
613 IF( irows.GT.1 )
THEN
614 CALL slacpy(
'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
615 $ vsl( ilo+1, ilo ), ldvsl )
617 CALL sorgqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
618 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
624 $
CALL slaset(
'Full', n, n, zero, one, vsr, ldvsr )
629 CALL sgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
630 $ ldvsl, vsr, ldvsr, ierr )
638 CALL shgeqz(
'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
639 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
640 $ work( iwrk ), lwork+1-iwrk, ierr )
642 IF( ierr.GT.0 .AND. ierr.LE.n )
THEN
644 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n )
THEN
662 CALL slascl(
'G', 0, 0, anrmto, anrm, n, 1, alphar, n,
664 CALL slascl(
'G', 0, 0, anrmto, anrm, n, 1, alphai, n,
668 $
CALL slascl(
'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
673 bwork( i ) = selctg( alphar( i ), alphai( i ), beta( i ) )
679 CALL stgsen( ijob, ilvsl, ilvsr, bwork, n, a, lda, b, ldb,
680 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
681 $ sdim, pl, pr, dif, work( iwrk ), lwork-iwrk+1,
682 $ iwork, liwork, ierr )
685 $ maxwrk = max( maxwrk, 2*sdim*( n-sdim ) )
686 IF( ierr.EQ.-22 )
THEN
692 IF( ijob.EQ.1 .OR. ijob.EQ.4 )
THEN
696 IF( ijob.EQ.2 .OR. ijob.EQ.4 )
THEN
697 rcondv( 1 ) = dif( 1 )
698 rcondv( 2 ) = dif( 2 )
710 $
CALL sggbak(
'P',
'L', n, ilo, ihi, work( ileft ),
711 $ work( iright ), n, vsl, ldvsl, ierr )
714 $
CALL sggbak(
'P',
'R', n, ilo, ihi, work( ileft ),
715 $ work( iright ), n, vsr, ldvsr, ierr )
723 IF( alphai( i ).NE.zero )
THEN
724 IF( ( alphar( i ) / safmax ).GT.( anrmto / anrm ) .OR.
725 $ ( safmin / alphar( i ) ).GT.( anrm / anrmto ) )
727 work( 1 ) = abs( a( i, i ) / alphar( i ) )
728 beta( i ) = beta( i )*work( 1 )
729 alphar( i ) = alphar( i )*work( 1 )
730 alphai( i ) = alphai( i )*work( 1 )
731 ELSE IF( ( alphai( i ) / safmax ).GT.( anrmto / anrm )
732 $ .OR. ( safmin / alphai( i ) ).GT.( anrm / anrmto ) )
734 work( 1 ) = abs( a( i, i+1 ) / alphai( i ) )
735 beta( i ) = beta( i )*work( 1 )
736 alphar( i ) = alphar( i )*work( 1 )
737 alphai( i ) = alphai( i )*work( 1 )
745 IF( alphai( i ).NE.zero )
THEN
746 IF( ( beta( i ) / safmax ).GT.( bnrmto / bnrm ) .OR.
747 $ ( safmin / beta( i ) ).GT.( bnrm / bnrmto ) )
THEN
748 work( 1 ) = abs( b( i, i ) / beta( i ) )
749 beta( i ) = beta( i )*work( 1 )
750 alphar( i ) = alphar( i )*work( 1 )
751 alphai( i ) = alphai( i )*work( 1 )
760 CALL slascl(
'H', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
761 CALL slascl(
'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
762 CALL slascl(
'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
766 CALL slascl(
'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
767 CALL slascl(
'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
779 cursl = selctg( alphar( i ), alphai( i ), beta( i ) )
780 IF( alphai( i ).EQ.zero )
THEN
784 IF( cursl .AND. .NOT.lastsl )
791 cursl = cursl .OR. lastsl
796 IF( cursl .AND. .NOT.lst2sl )
subroutine sggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
SGGBAL
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 sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
subroutine shgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
SHGEQZ
subroutine xerbla(SRNAME, INFO)
XERBLA
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 sggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
SGGBAK
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 slabad(SMALL, LARGE)
SLABAD
subroutine sgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
SGGHRD
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.