363 SUBROUTINE dggesx( 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 DOUBLE PRECISION a( lda, * ), alphai( * ), alphar( * ),
382 $
b( ldb, * ), beta( * ), rconde( 2 ),
383 $ rcondv( 2 ), vsl( ldvsl, * ), vsr( ldvsr, * ),
394 DOUBLE PRECISION zero, one
395 parameter( zero = 0.0d+0, one = 1.0d+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 DOUBLE PRECISION anrm, anrmto, bignum, bnrm, bnrmto, eps, pl,
405 $ pr, safmax, safmin, smlnum
408 DOUBLE PRECISION dif( 2 )
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,
'DGEQRF',
' ', n, 1, n, 0 )
502 maxwrk = max( maxwrk, minwrk - n +
503 $ n*
ilaenv( 1,
'DORMQR',
' ', n, 1, n, -1 ) )
505 maxwrk = max( maxwrk, minwrk - n +
506 $ n*
ilaenv( 1,
'DORGQR',
' ', 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(
'DGGESX', -info )
534 ELSE IF (lquery)
THEN
549 safmax = one / safmin
550 CALL
dlabad( safmin, safmax )
551 smlnum = sqrt( safmin ) / eps
552 bignum = one / smlnum
556 anrm =
dlange(
'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
dlascl(
'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
570 bnrm =
dlange(
'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
dlascl(
'G', 0, 0, bnrm, bnrmto, n, n,
b, ldb, ierr )
588 CALL
dggbal(
'P', n, a, lda,
b, ldb, ilo, ihi, work( ileft ),
589 $ work( iright ), work( iwrk ), ierr )
594 irows = ihi + 1 - ilo
598 CALL
dgeqrf( irows, icols,
b( ilo, ilo ), ldb, work( itau ),
599 $ work( iwrk ), lwork+1-iwrk, ierr )
604 CALL
dormqr(
'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
dlaset(
'Full', n, n, zero, one, vsl, ldvsl )
613 IF( irows.GT.1 )
THEN
614 CALL
dlacpy(
'L', irows-1, irows-1,
b( ilo+1, ilo ), ldb,
615 $ vsl( ilo+1, ilo ), ldvsl )
617 CALL
dorgqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
618 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
624 $ CALL
dlaset(
'Full', n, n, zero, one, vsr, ldvsr )
629 CALL
dgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda,
b, ldb, vsl,
630 $ ldvsl, vsr, ldvsr, ierr )
638 CALL
dhgeqz(
'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
dlascl(
'G', 0, 0, anrmto, anrm, n, 1, alphar, n,
664 CALL
dlascl(
'G', 0, 0, anrmto, anrm, n, 1, alphai, n,
668 $ CALL
dlascl(
'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
673 bwork( i ) = selctg( alphar( i ), alphai( i ), beta( i ) )
679 CALL
dtgsen( 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
dggbak(
'P',
'L', n, ilo, ihi, work( ileft ),
711 $ work( iright ), n, vsl, ldvsl, ierr )
714 $ CALL
dggbak(
'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 ) )
THEN
726 work( 1 ) = abs( a( i, i ) / alphar( i ) )
727 beta( i ) = beta( i )*work( 1 )
728 alphar( i ) = alphar( i )*work( 1 )
729 alphai( i ) = alphai( i )*work( 1 )
730 ELSE IF( ( alphai( i ) / safmax ).GT.
731 $ ( anrmto / anrm ) .OR.
732 $ ( 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
dlascl(
'H', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
761 CALL
dlascl(
'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
762 CALL
dlascl(
'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
766 CALL
dlascl(
'U', 0, 0, bnrmto, bnrm, n, n,
b, ldb, ierr )
767 CALL
dlascl(
'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 dggesx(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)
DGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE...
subroutine dgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
DGGHRD
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine dlabad(SMALL, LARGE)
DLABAD
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
logical function lsame(CA, CB)
LSAME
subroutine dorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGQR
subroutine dhgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
DHGEQZ
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine dggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
DGGBAK
double precision function dlamch(CMACH)
DLAMCH
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine dggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
DGGBAL
subroutine dtgsen(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)
DTGSEN