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
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
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 slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
LOGICAL function lsame(CA, CB)
LSAME
subroutine sggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
SGGBAL
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
REAL function slamch(CMACH)
SLAMCH
subroutine xerbla(SRNAME, INFO)
XERBLA
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
subroutine shgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
SHGEQZ
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 slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
INTEGER function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
REAL function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
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 slabad(SMALL, LARGE)
SLABAD
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
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 sggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
SGGBAK