359 SUBROUTINE dggesx( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A,
361 $ B, LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL,
362 $ VSR, LDVSR, RCONDE, RCONDV, WORK, LWORK, IWORK,
363 $ LIWORK, BWORK, INFO )
370 CHARACTER JOBVSL, JOBVSR, SENSE, SORT
371 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N,
377 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
378 $ B( LDB, * ), BETA( * ), RCONDE( 2 ),
379 $ rcondv( 2 ), vsl( ldvsl, * ), vsr( ldvsr, * ),
390 DOUBLE PRECISION ZERO, ONE
391 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0 )
394 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
395 $ lquery, lst2sl, wantsb, wantse, wantsn, wantst,
397 INTEGER I, ICOLS, IERR, IHI, IJOB, IJOBVL, IJOBVR,
398 $ ileft, ilo, ip, iright, irows, itau, iwrk,
399 $ liwmin, lwrk, maxwrk, minwrk
400 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PL,
401 $ pr, safmax, safmin, smlnum
404 DOUBLE PRECISION DIF( 2 )
407 EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ,
414 DOUBLE PRECISION DLAMCH, DLANGE
415 EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
418 INTRINSIC abs, max, sqrt
424 IF( lsame( jobvsl,
'N' ) )
THEN
427 ELSE IF( lsame( jobvsl,
'V' ) )
THEN
435 IF( lsame( jobvsr,
'N' ) )
THEN
438 ELSE IF( lsame( jobvsr,
'V' ) )
THEN
446 wantst = lsame( sort,
'S' )
447 wantsn = lsame( sense,
'N' )
448 wantse = lsame( sense,
'E' )
449 wantsv = lsame( sense,
'V' )
450 wantsb = lsame( sense,
'B' )
451 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
454 ELSE IF( wantse )
THEN
456 ELSE IF( wantsv )
THEN
458 ELSE IF( wantsb )
THEN
465 IF( ijobvl.LE.0 )
THEN
467 ELSE IF( ijobvr.LE.0 )
THEN
469 ELSE IF( ( .NOT.wantst ) .AND.
470 $ ( .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,
'DGEQRF',
' ', n, 1, n, 0 )
499 maxwrk = max( maxwrk, minwrk - n +
500 $ n*ilaenv( 1,
'DORMQR',
' ', n, 1, n, -1 ) )
502 maxwrk = max( maxwrk, minwrk - n +
503 $ n*ilaenv( 1,
'DORGQR',
' ', 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(
'DGGESX', -info )
531 ELSE IF (lquery)
THEN
545 safmin = dlamch(
'S' )
546 safmax = one / safmin
547 smlnum = sqrt( safmin ) / eps
548 bignum = one / smlnum
552 anrm = dlange(
'M', n, n, a, lda, work )
554 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
557 ELSE IF( anrm.GT.bignum )
THEN
562 $
CALL dlascl(
'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
566 bnrm = dlange(
'M', n, n, b, ldb, work )
568 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum )
THEN
571 ELSE IF( bnrm.GT.bignum )
THEN
576 $
CALL dlascl(
'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
584 CALL dggbal(
'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
585 $ work( iright ), work( iwrk ), ierr )
590 irows = ihi + 1 - ilo
594 CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
595 $ work( iwrk ), lwork+1-iwrk, ierr )
600 CALL dormqr(
'L',
'T', irows, icols, irows, b( ilo, ilo ), ldb,
601 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
602 $ lwork+1-iwrk, ierr )
608 CALL dlaset(
'Full', n, n, zero, one, vsl, ldvsl )
609 IF( irows.GT.1 )
THEN
610 CALL dlacpy(
'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
611 $ vsl( ilo+1, ilo ), ldvsl )
613 CALL dorgqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
614 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
620 $
CALL dlaset(
'Full', n, n, zero, one, vsr, ldvsr )
625 CALL dgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
626 $ ldvsl, vsr, ldvsr, ierr )
634 CALL dhgeqz(
'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
635 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
636 $ work( iwrk ), lwork+1-iwrk, ierr )
638 IF( ierr.GT.0 .AND. ierr.LE.n )
THEN
640 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n )
THEN
658 CALL dlascl(
'G', 0, 0, anrmto, anrm, n, 1, alphar, n,
660 CALL dlascl(
'G', 0, 0, anrmto, anrm, n, 1, alphai, n,
664 $
CALL dlascl(
'G', 0, 0, bnrmto, bnrm, n, 1, beta, n,
670 bwork( i ) = selctg( alphar( i ), alphai( i ),
677 CALL dtgsen( ijob, ilvsl, ilvsr, bwork, n, a, lda, b, ldb,
678 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
679 $ sdim, pl, pr, dif, work( iwrk ), lwork-iwrk+1,
680 $ iwork, liwork, ierr )
683 $ maxwrk = max( maxwrk, 2*sdim*( n-sdim ) )
684 IF( ierr.EQ.-22 )
THEN
690 IF( ijob.EQ.1 .OR. ijob.EQ.4 )
THEN
694 IF( ijob.EQ.2 .OR. ijob.EQ.4 )
THEN
695 rcondv( 1 ) = dif( 1 )
696 rcondv( 2 ) = dif( 2 )
708 $
CALL dggbak(
'P',
'L', n, ilo, ihi, work( ileft ),
709 $ work( iright ), n, vsl, ldvsl, ierr )
712 $
CALL dggbak(
'P',
'R', n, ilo, ihi, work( ileft ),
713 $ work( iright ), n, vsr, ldvsr, ierr )
721 IF( alphai( i ).NE.zero )
THEN
722 IF( ( alphar( i ) / safmax ).GT.( anrmto / anrm ) .OR.
723 $ ( safmin / alphar( i ) ).GT.( anrm / anrmto ) )
THEN
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.
729 $ ( anrmto / anrm ) .OR.
730 $ ( safmin / alphai( i ) ).GT.( anrm / anrmto ) )
732 work( 1 ) = abs( a( i, i+1 ) / alphai( i ) )
733 beta( i ) = beta( i )*work( 1 )
734 alphar( i ) = alphar( i )*work( 1 )
735 alphai( i ) = alphai( i )*work( 1 )
743 IF( alphai( i ).NE.zero )
THEN
744 IF( ( beta( i ) / safmax ).GT.( bnrmto / bnrm ) .OR.
745 $ ( safmin / beta( i ) ).GT.( bnrm / bnrmto ) )
THEN
746 work( 1 ) = abs( b( i, i ) / beta( i ) )
747 beta( i ) = beta( i )*work( 1 )
748 alphar( i ) = alphar( i )*work( 1 )
749 alphai( i ) = alphai( i )*work( 1 )
758 CALL dlascl(
'H', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
759 CALL dlascl(
'G', 0, 0, anrmto, anrm, n, 1, alphar, n,
761 CALL dlascl(
'G', 0, 0, anrmto, anrm, n, 1, alphai, n,
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 )