302 SUBROUTINE dgegv( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR,
303 $ ALPHAI, BETA, VL, LDVL, VR, LDVR, WORK, LWORK,
311 CHARACTER JOBVL, JOBVR
312 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
315 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
316 $ B( LDB, * ), BETA( * ), VL( LDVL, * ),
317 $ vr( ldvr, * ), work( * )
323 DOUBLE PRECISION ZERO, ONE
324 PARAMETER ( ZERO = 0.0d0, one = 1.0d0 )
327 LOGICAL ILIMIT, ILV, ILVL, ILVR, LQUERY
329 INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
330 $ in, iright, irows, itau, iwork, jc, jr, lopt,
331 $ lwkmin, lwkopt, nb, nb1, nb2, nb3
332 DOUBLE PRECISION ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM,
333 $ BNRM1, BNRM2, EPS, ONEPLS, SAFMAX, SAFMIN,
334 $ salfai, salfar, sbeta, scale, temp
346 DOUBLE PRECISION DLAMCH, DLANGE
347 EXTERNAL lsame, ilaenv, dlamch, dlange
350 INTRINSIC abs, int, max
356 IF( lsame( jobvl,
'N' ) )
THEN
359 ELSE IF( lsame( jobvl,
'V' ) )
THEN
367 IF( lsame( jobvr,
'N' ) )
THEN
370 ELSE IF( lsame( jobvr,
'V' ) )
THEN
381 lwkmin = max( 8*n, 1 )
384 lquery = ( lwork.EQ.-1 )
386 IF( ijobvl.LE.0 )
THEN
388 ELSE IF( ijobvr.LE.0 )
THEN
390 ELSE IF( n.LT.0 )
THEN
392 ELSE IF( lda.LT.max( 1, n ) )
THEN
394 ELSE IF( ldb.LT.max( 1, n ) )
THEN
396 ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) )
THEN
398 ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) )
THEN
400 ELSE IF( lwork.LT.lwkmin .AND. .NOT.lquery )
THEN
405 nb1 = ilaenv( 1,
'DGEQRF',
' ', n, n, -1, -1 )
406 nb2 = ilaenv( 1,
'DORMQR',
' ', n, n, n, -1 )
407 nb3 = ilaenv( 1,
'DORGQR',
' ', n, n, n, -1 )
408 nb = max( nb1, nb2, nb3 )
409 lopt = 2*n + max( 6*n, n*( nb+1 ) )
414 CALL xerbla(
'DGEGV ', -info )
416 ELSE IF( lquery )
THEN
427 eps = dlamch(
'E' )*dlamch(
'B' )
428 safmin = dlamch(
'S' )
429 safmin = safmin + safmin
430 safmax = one / safmin
431 onepls = one + ( 4*eps )
435 anrm = dlange(
'M', n, n, a, lda, work )
438 IF( anrm.LT.one )
THEN
439 IF( safmax*anrm.LT.one )
THEN
445 IF( anrm.GT.zero )
THEN
446 CALL dlascl(
'G', -1, -1, anrm, one, n, n, a, lda, iinfo )
447 IF( iinfo.NE.0 )
THEN
455 bnrm = dlange(
'M', n, n, b, ldb, work )
458 IF( bnrm.LT.one )
THEN
459 IF( safmax*bnrm.LT.one )
THEN
465 IF( bnrm.GT.zero )
THEN
466 CALL dlascl(
'G', -1, -1, bnrm, one, n, n, b, ldb, iinfo )
467 IF( iinfo.NE.0 )
THEN
480 CALL dggbal(
'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
481 $ work( iright ), work( iwork ), iinfo )
482 IF( iinfo.NE.0 )
THEN
491 irows = ihi + 1 - ilo
499 CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
500 $ work( iwork ), lwork+1-iwork, iinfo )
502 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
503 IF( iinfo.NE.0 )
THEN
508 CALL dormqr(
'L',
'T', irows, icols, irows, b( ilo, ilo ), ldb,
509 $ work( itau ), a( ilo, ilo ), lda, work( iwork ),
510 $ lwork+1-iwork, iinfo )
512 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
513 IF( iinfo.NE.0 )
THEN
519 CALL dlaset(
'Full', n, n, zero, one, vl, ldvl )
520 CALL dlacpy(
'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
521 $ vl( ilo+1, ilo ), ldvl )
522 CALL dorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
523 $ work( itau ), work( iwork ), lwork+1-iwork,
526 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
527 IF( iinfo.NE.0 )
THEN
534 $
CALL dlaset(
'Full', n, n, zero, one, vr, ldvr )
542 CALL dgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
543 $ ldvl, vr, ldvr, iinfo )
545 CALL dgghrd(
'N',
'N', irows, 1, irows, a( ilo, ilo ), lda,
546 $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, iinfo )
548 IF( iinfo.NE.0 )
THEN
563 CALL dhgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
564 $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
565 $ work( iwork ), lwork+1-iwork, iinfo )
567 $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
568 IF( iinfo.NE.0 )
THEN
569 IF( iinfo.GT.0 .AND. iinfo.LE.n )
THEN
571 ELSE IF( iinfo.GT.n .AND. iinfo.LE.2*n )
THEN
593 CALL dtgevc( chtemp,
'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
594 $ vr, ldvr, n, in, work( iwork ), iinfo )
595 IF( iinfo.NE.0 )
THEN
603 CALL dggbak(
'P',
'L', n, ilo, ihi, work( ileft ),
604 $ work( iright ), n, vl, ldvl, iinfo )
605 IF( iinfo.NE.0 )
THEN
610 IF( alphai( jc ).LT.zero )
613 IF( alphai( jc ).EQ.zero )
THEN
615 temp = max( temp, abs( vl( jr, jc ) ) )
619 temp = max( temp, abs( vl( jr, jc ) )+
620 $ abs( vl( jr, jc+1 ) ) )
626 IF( alphai( jc ).EQ.zero )
THEN
628 vl( jr, jc ) = vl( jr, jc )*temp
632 vl( jr, jc ) = vl( jr, jc )*temp
633 vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
639 CALL dggbak(
'P',
'R', n, ilo, ihi, work( ileft ),
640 $ work( iright ), n, vr, ldvr, iinfo )
641 IF( iinfo.NE.0 )
THEN
646 IF( alphai( jc ).LT.zero )
649 IF( alphai( jc ).EQ.zero )
THEN
651 temp = max( temp, abs( vr( jr, jc ) ) )
655 temp = max( temp, abs( vr( jr, jc ) )+
656 $ abs( vr( jr, jc+1 ) ) )
662 IF( alphai( jc ).EQ.zero )
THEN
664 vr( jr, jc ) = vr( jr, jc )*temp
668 vr( jr, jc ) = vr( jr, jc )*temp
669 vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
688 absar = abs( alphar( jc ) )
689 absai = abs( alphai( jc ) )
690 absb = abs( beta( jc ) )
691 salfar = anrm*alphar( jc )
692 salfai = anrm*alphai( jc )
693 sbeta = bnrm*beta( jc )
699 IF( abs( salfai ).LT.safmin .AND. absai.GE.
700 $ max( safmin, eps*absar, eps*absb ) )
THEN
702 scale = ( onepls*safmin / anrm1 ) /
703 $ max( onepls*safmin, anrm2*absai )
705 ELSE IF( salfai.EQ.zero )
THEN
710 IF( alphai( jc ).LT.zero .AND. jc.GT.1 )
THEN
711 alphai( jc-1 ) = zero
712 ELSE IF( alphai( jc ).GT.zero .AND. jc.LT.n )
THEN
713 alphai( jc+1 ) = zero
719 IF( abs( salfar ).LT.safmin .AND. absar.GE.
720 $ max( safmin, eps*absai, eps*absb ) )
THEN
722 scale = max( scale, ( onepls*safmin / anrm1 ) /
723 $ max( onepls*safmin, anrm2*absar ) )
728 IF( abs( sbeta ).LT.safmin .AND. absb.GE.
729 $ max( safmin, eps*absar, eps*absai ) )
THEN
731 scale = max( scale, ( onepls*safmin / bnrm1 ) /
732 $ max( onepls*safmin, bnrm2*absb ) )
738 temp = ( scale*safmin )*max( abs( salfar ), abs( salfai ),
741 $ scale = scale / temp
749 salfar = ( scale*alphar( jc ) )*anrm
750 salfai = ( scale*alphai( jc ) )*anrm
751 sbeta = ( scale*beta( jc ) )*bnrm
753 alphar( jc ) = salfar
754 alphai( jc ) = salfai