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 )
418 DOUBLE PRECISION DLAMCH, DLANGE
419 EXTERNAL lsame, ilaenv, dlamch, dlange
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
548 safmin = dlamch(
'S' )
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 dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dhgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
DHGEQZ
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 dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
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
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 dorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGQR
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
DGGBAL
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
subroutine dggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
DGGBAK
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