293 SUBROUTINE dtgevc( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
294 $ LDVL, VR, LDVR, MM, M, WORK, INFO )
301 CHARACTER HOWMNY, SIDE
302 INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
306 DOUBLE PRECISION P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
307 $ vr( ldvr, * ), work( * )
314 DOUBLE PRECISION ZERO, ONE, SAFETY
315 parameter( zero = 0.0d+0, one = 1.0d+0,
319 LOGICAL COMPL, COMPR, IL2BY2, ILABAD, ILALL, ILBACK,
320 $ ilbbad, ilcomp, ilcplx, lsa, lsb
321 INTEGER I, IBEG, IEIG, IEND, IHWMNY, IINFO, IM, ISIDE,
322 $ j, ja, jc, je, jr, jw, na, nw
323 DOUBLE PRECISION ACOEF, ACOEFA, ANORM, ASCALE, BCOEFA, BCOEFI,
324 $ bcoefr, big, bignum, bnorm, bscale, cim2a,
325 $ cim2b, cimaga, cimagb, cre2a, cre2b, creala,
326 $ crealb, dmin, safmin, salfar, sbeta, scale,
327 $ small, temp, temp2, temp2i, temp2r, ulp, xmax,
331 DOUBLE PRECISION BDIAG( 2 ), SUM( 2, 2 ), SUMS( 2, 2 ),
336 DOUBLE PRECISION DLAMCH
337 EXTERNAL lsame, dlamch
343 INTRINSIC abs, max, min
349 IF( lsame( howmny,
'A' ) )
THEN
353 ELSE IF( lsame( howmny,
'S' ) )
THEN
357 ELSE IF( lsame( howmny,
'B' ) )
THEN
366 IF( lsame( side,
'R' ) )
THEN
370 ELSE IF( lsame( side,
'L' ) )
THEN
374 ELSE IF( lsame( side,
'B' ) )
THEN
383 IF( iside.LT.0 )
THEN
385 ELSE IF( ihwmny.LT.0 )
THEN
387 ELSE IF( n.LT.0 )
THEN
389 ELSE IF( lds.LT.max( 1, n ) )
THEN
391 ELSE IF( ldp.LT.max( 1, n ) )
THEN
395 CALL xerbla(
'DTGEVC', -info )
401 IF( .NOT.ilall )
THEN
410 IF( s( j+1, j ).NE.zero )
414 IF(
SELECT( j ) .OR.
SELECT( j+1 ) )
430 IF( s( j+1, j ).NE.zero )
THEN
431 IF( p( j, j ).EQ.zero .OR. p( j+1, j+1 ).EQ.zero .OR.
432 $ p( j, j+1 ).NE.zero )ilbbad = .true.
434 IF( s( j+2, j+1 ).NE.zero )
442 ELSE IF( ilbbad )
THEN
444 ELSE IF( compl .AND. ldvl.LT.n .OR. ldvl.LT.1 )
THEN
446 ELSE IF( compr .AND. ldvr.LT.n .OR. ldvr.LT.1 )
THEN
448 ELSE IF( mm.LT.im )
THEN
452 CALL xerbla(
'DTGEVC', -info )
464 safmin = dlamch(
'Safe minimum' )
466 CALL dlabad( safmin, big )
467 ulp = dlamch(
'Epsilon' )*dlamch(
'Base' )
468 small = safmin*n / ulp
470 bignum = one / ( safmin*n )
477 anorm = abs( s( 1, 1 ) )
479 $ anorm = anorm + abs( s( 2, 1 ) )
480 bnorm = abs( p( 1, 1 ) )
487 IF( s( j, j-1 ).EQ.zero )
THEN
493 temp = temp + abs( s( i, j ) )
494 temp2 = temp2 + abs( p( i, j ) )
498 DO 40 i = iend + 1, min( j+1, n )
499 temp = temp + abs( s( i, j ) )
500 temp2 = temp2 + abs( p( i, j ) )
502 anorm = max( anorm, temp )
503 bnorm = max( bnorm, temp2 )
506 ascale = one / max( anorm, safmin )
507 bscale = one / max( bnorm, safmin )
530 IF( s( je+1, je ).NE.zero )
THEN
537 ELSE IF( ilcplx )
THEN
538 ilcomp =
SELECT( je ) .OR.
SELECT( je+1 )
540 ilcomp =
SELECT( je )
548 IF( .NOT.ilcplx )
THEN
549 IF( abs( s( je, je ) ).LE.safmin .AND.
550 $ abs( p( je, je ) ).LE.safmin )
THEN
556 vl( jr, ieig ) = zero
558 vl( ieig, ieig ) = one
566 work( 2*n+jr ) = zero
573 IF( .NOT.ilcplx )
THEN
577 temp = one / max( abs( s( je, je ) )*ascale,
578 $ abs( p( je, je ) )*bscale, safmin )
579 salfar = ( temp*s( je, je ) )*ascale
580 sbeta = ( temp*p( je, je ) )*bscale
582 bcoefr = salfar*bscale
588 lsa = abs( sbeta ).GE.safmin .AND. abs( acoef ).LT.small
589 lsb = abs( salfar ).GE.safmin .AND. abs( bcoefr ).LT.
592 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
594 $ scale = max( scale, ( small / abs( salfar ) )*
595 $ min( bnorm, big ) )
596 IF( lsa .OR. lsb )
THEN
597 scale = min( scale, one /
598 $ ( safmin*max( one, abs( acoef ),
599 $ abs( bcoefr ) ) ) )
601 acoef = ascale*( scale*sbeta )
606 bcoefr = bscale*( scale*salfar )
608 bcoefr = scale*bcoefr
611 acoefa = abs( acoef )
612 bcoefa = abs( bcoefr )
622 CALL dlag2( s( je, je ), lds, p( je, je ), ldp,
623 $ safmin*safety, acoef, temp, bcoefr, temp2,
626 IF( bcoefi.EQ.zero )
THEN
633 acoefa = abs( acoef )
634 bcoefa = abs( bcoefr ) + abs( bcoefi )
636 IF( acoefa*ulp.LT.safmin .AND. acoefa.GE.safmin )
637 $ scale = ( safmin / ulp ) / acoefa
638 IF( bcoefa*ulp.LT.safmin .AND. bcoefa.GE.safmin )
639 $ scale = max( scale, ( safmin / ulp ) / bcoefa )
640 IF( safmin*acoefa.GT.ascale )
641 $ scale = ascale / ( safmin*acoefa )
642 IF( safmin*bcoefa.GT.bscale )
643 $ scale = min( scale, bscale / ( safmin*bcoefa ) )
644 IF( scale.NE.one )
THEN
646 acoefa = abs( acoef )
647 bcoefr = scale*bcoefr
648 bcoefi = scale*bcoefi
649 bcoefa = abs( bcoefr ) + abs( bcoefi )
654 temp = acoef*s( je+1, je )
655 temp2r = acoef*s( je, je ) - bcoefr*p( je, je )
656 temp2i = -bcoefi*p( je, je )
657 IF( abs( temp ).GT.abs( temp2r )+abs( temp2i ) )
THEN
659 work( 3*n+je ) = zero
660 work( 2*n+je+1 ) = -temp2r / temp
661 work( 3*n+je+1 ) = -temp2i / temp
663 work( 2*n+je+1 ) = one
664 work( 3*n+je+1 ) = zero
665 temp = acoef*s( je, je+1 )
666 work( 2*n+je ) = ( bcoefr*p( je+1, je+1 )-acoef*
667 $ s( je+1, je+1 ) ) / temp
668 work( 3*n+je ) = bcoefi*p( je+1, je+1 ) / temp
670 xmax = max( abs( work( 2*n+je ) )+abs( work( 3*n+je ) ),
671 $ abs( work( 2*n+je+1 ) )+abs( work( 3*n+je+1 ) ) )
674 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
684 DO 160 j = je + nw, n
691 bdiag( 1 ) = p( j, j )
693 IF( s( j+1, j ).NE.zero )
THEN
695 bdiag( 2 ) = p( j+1, j+1 )
702 xscale = one / max( one, xmax )
703 temp = max( work( j ), work( n+j ),
704 $ acoefa*work( j )+bcoefa*work( n+j ) )
706 $ temp = max( temp, work( j+1 ), work( n+j+1 ),
707 $ acoefa*work( j+1 )+bcoefa*work( n+j+1 ) )
708 IF( temp.GT.bignum*xscale )
THEN
711 work( ( jw+2 )*n+jr ) = xscale*
712 $ work( ( jw+2 )*n+jr )
736 sums( ja, jw ) = zero
737 sump( ja, jw ) = zero
739 DO 100 jr = je, j - 1
740 sums( ja, jw ) = sums( ja, jw ) +
742 $ work( ( jw+1 )*n+jr )
743 sump( ja, jw ) = sump( ja, jw ) +
745 $ work( ( jw+1 )*n+jr )
752 sum( ja, 1 ) = -acoef*sums( ja, 1 ) +
753 $ bcoefr*sump( ja, 1 ) -
754 $ bcoefi*sump( ja, 2 )
755 sum( ja, 2 ) = -acoef*sums( ja, 2 ) +
756 $ bcoefr*sump( ja, 2 ) +
757 $ bcoefi*sump( ja, 1 )
759 sum( ja, 1 ) = -acoef*sums( ja, 1 ) +
760 $ bcoefr*sump( ja, 1 )
768 CALL dlaln2( .true., na, nw, dmin, acoef, s( j, j ), lds,
769 $ bdiag( 1 ), bdiag( 2 ), sum, 2, bcoefr,
770 $ bcoefi, work( 2*n+j ), n, scale, temp,
772 IF( scale.LT.one )
THEN
773 DO 150 jw = 0, nw - 1
774 DO 140 jr = je, j - 1
775 work( ( jw+2 )*n+jr ) = scale*
776 $ work( ( jw+2 )*n+jr )
781 xmax = max( xmax, temp )
789 DO 170 jw = 0, nw - 1
790 CALL dgemv(
'N', n, n+1-je, one, vl( 1, je ), ldvl,
791 $ work( ( jw+2 )*n+je ), 1, zero,
792 $ work( ( jw+4 )*n+1 ), 1 )
794 CALL dlacpy(
' ', n, nw, work( 4*n+1 ), n, vl( 1, je ),
798 CALL dlacpy(
' ', n, nw, work( 2*n+1 ), n, vl( 1, ieig ),
808 xmax = max( xmax, abs( vl( j, ieig ) )+
809 $ abs( vl( j, ieig+1 ) ) )
813 xmax = max( xmax, abs( vl( j, ieig ) ) )
817 IF( xmax.GT.safmin )
THEN
820 DO 210 jw = 0, nw - 1
822 vl( jr, ieig+jw ) = xscale*vl( jr, ieig+jw )
855 IF( s( je, je-1 ).NE.zero )
THEN
862 ELSE IF( ilcplx )
THEN
863 ilcomp =
SELECT( je ) .OR.
SELECT( je-1 )
865 ilcomp =
SELECT( je )
873 IF( .NOT.ilcplx )
THEN
874 IF( abs( s( je, je ) ).LE.safmin .AND.
875 $ abs( p( je, je ) ).LE.safmin )
THEN
881 vr( jr, ieig ) = zero
883 vr( ieig, ieig ) = one
890 DO 250 jw = 0, nw - 1
892 work( ( jw+2 )*n+jr ) = zero
900 IF( .NOT.ilcplx )
THEN
904 temp = one / max( abs( s( je, je ) )*ascale,
905 $ abs( p( je, je ) )*bscale, safmin )
906 salfar = ( temp*s( je, je ) )*ascale
907 sbeta = ( temp*p( je, je ) )*bscale
909 bcoefr = salfar*bscale
915 lsa = abs( sbeta ).GE.safmin .AND. abs( acoef ).LT.small
916 lsb = abs( salfar ).GE.safmin .AND. abs( bcoefr ).LT.
919 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
921 $ scale = max( scale, ( small / abs( salfar ) )*
922 $ min( bnorm, big ) )
923 IF( lsa .OR. lsb )
THEN
924 scale = min( scale, one /
925 $ ( safmin*max( one, abs( acoef ),
926 $ abs( bcoefr ) ) ) )
928 acoef = ascale*( scale*sbeta )
933 bcoefr = bscale*( scale*salfar )
935 bcoefr = scale*bcoefr
938 acoefa = abs( acoef )
939 bcoefa = abs( bcoefr )
949 DO 260 jr = 1, je - 1
950 work( 2*n+jr ) = bcoefr*p( jr, je ) -
957 CALL dlag2( s( je-1, je-1 ), lds, p( je-1, je-1 ), ldp,
958 $ safmin*safety, acoef, temp, bcoefr, temp2,
960 IF( bcoefi.EQ.zero )
THEN
967 acoefa = abs( acoef )
968 bcoefa = abs( bcoefr ) + abs( bcoefi )
970 IF( acoefa*ulp.LT.safmin .AND. acoefa.GE.safmin )
971 $ scale = ( safmin / ulp ) / acoefa
972 IF( bcoefa*ulp.LT.safmin .AND. bcoefa.GE.safmin )
973 $ scale = max( scale, ( safmin / ulp ) / bcoefa )
974 IF( safmin*acoefa.GT.ascale )
975 $ scale = ascale / ( safmin*acoefa )
976 IF( safmin*bcoefa.GT.bscale )
977 $ scale = min( scale, bscale / ( safmin*bcoefa ) )
978 IF( scale.NE.one )
THEN
980 acoefa = abs( acoef )
981 bcoefr = scale*bcoefr
982 bcoefi = scale*bcoefi
983 bcoefa = abs( bcoefr ) + abs( bcoefi )
989 temp = acoef*s( je, je-1 )
990 temp2r = acoef*s( je, je ) - bcoefr*p( je, je )
991 temp2i = -bcoefi*p( je, je )
992 IF( abs( temp ).GE.abs( temp2r )+abs( temp2i ) )
THEN
994 work( 3*n+je ) = zero
995 work( 2*n+je-1 ) = -temp2r / temp
996 work( 3*n+je-1 ) = -temp2i / temp
998 work( 2*n+je-1 ) = one
999 work( 3*n+je-1 ) = zero
1000 temp = acoef*s( je-1, je )
1001 work( 2*n+je ) = ( bcoefr*p( je-1, je-1 )-acoef*
1002 $ s( je-1, je-1 ) ) / temp
1003 work( 3*n+je ) = bcoefi*p( je-1, je-1 ) / temp
1006 xmax = max( abs( work( 2*n+je ) )+abs( work( 3*n+je ) ),
1007 $ abs( work( 2*n+je-1 ) )+abs( work( 3*n+je-1 ) ) )
1012 creala = acoef*work( 2*n+je-1 )
1013 cimaga = acoef*work( 3*n+je-1 )
1014 crealb = bcoefr*work( 2*n+je-1 ) -
1015 $ bcoefi*work( 3*n+je-1 )
1016 cimagb = bcoefi*work( 2*n+je-1 ) +
1017 $ bcoefr*work( 3*n+je-1 )
1018 cre2a = acoef*work( 2*n+je )
1019 cim2a = acoef*work( 3*n+je )
1020 cre2b = bcoefr*work( 2*n+je ) - bcoefi*work( 3*n+je )
1021 cim2b = bcoefi*work( 2*n+je ) + bcoefr*work( 3*n+je )
1022 DO 270 jr = 1, je - 2
1023 work( 2*n+jr ) = -creala*s( jr, je-1 ) +
1024 $ crealb*p( jr, je-1 ) -
1025 $ cre2a*s( jr, je ) + cre2b*p( jr, je )
1026 work( 3*n+jr ) = -cimaga*s( jr, je-1 ) +
1027 $ cimagb*p( jr, je-1 ) -
1028 $ cim2a*s( jr, je ) + cim2b*p( jr, je )
1032 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
1037 DO 370 j = je - nw, 1, -1
1042 IF( .NOT.il2by2 .AND. j.GT.1 )
THEN
1043 IF( s( j, j-1 ).NE.zero )
THEN
1048 bdiag( 1 ) = p( j, j )
1051 bdiag( 2 ) = p( j+1, j+1 )
1058 CALL dlaln2( .false., na, nw, dmin, acoef, s( j, j ),
1059 $ lds, bdiag( 1 ), bdiag( 2 ), work( 2*n+j ),
1060 $ n, bcoefr, bcoefi, sum, 2, scale, temp,
1062 IF( scale.LT.one )
THEN
1064 DO 290 jw = 0, nw - 1
1066 work( ( jw+2 )*n+jr ) = scale*
1067 $ work( ( jw+2 )*n+jr )
1071 xmax = max( scale*xmax, temp )
1075 work( ( jw+1 )*n+j+ja-1 ) = sum( ja, jw )
1085 xscale = one / max( one, xmax )
1086 temp = acoefa*work( j ) + bcoefa*work( n+j )
1088 $ temp = max( temp, acoefa*work( j+1 )+bcoefa*
1090 temp = max( temp, acoefa, bcoefa )
1091 IF( temp.GT.bignum*xscale )
THEN
1093 DO 330 jw = 0, nw - 1
1095 work( ( jw+2 )*n+jr ) = xscale*
1096 $ work( ( jw+2 )*n+jr )
1109 creala = acoef*work( 2*n+j+ja-1 )
1110 cimaga = acoef*work( 3*n+j+ja-1 )
1111 crealb = bcoefr*work( 2*n+j+ja-1 ) -
1112 $ bcoefi*work( 3*n+j+ja-1 )
1113 cimagb = bcoefi*work( 2*n+j+ja-1 ) +
1114 $ bcoefr*work( 3*n+j+ja-1 )
1115 DO 340 jr = 1, j - 1
1116 work( 2*n+jr ) = work( 2*n+jr ) -
1117 $ creala*s( jr, j+ja-1 ) +
1118 $ crealb*p( jr, j+ja-1 )
1119 work( 3*n+jr ) = work( 3*n+jr ) -
1120 $ cimaga*s( jr, j+ja-1 ) +
1121 $ cimagb*p( jr, j+ja-1 )
1124 creala = acoef*work( 2*n+j+ja-1 )
1125 crealb = bcoefr*work( 2*n+j+ja-1 )
1126 DO 350 jr = 1, j - 1
1127 work( 2*n+jr ) = work( 2*n+jr ) -
1128 $ creala*s( jr, j+ja-1 ) +
1129 $ crealb*p( jr, j+ja-1 )
1144 DO 410 jw = 0, nw - 1
1146 work( ( jw+4 )*n+jr ) = work( ( jw+2 )*n+1 )*
1156 work( ( jw+4 )*n+jr ) = work( ( jw+4 )*n+jr ) +
1157 $ work( ( jw+2 )*n+jc )*vr( jr, jc )
1162 DO 430 jw = 0, nw - 1
1164 vr( jr, ieig+jw ) = work( ( jw+4 )*n+jr )
1170 DO 450 jw = 0, nw - 1
1172 vr( jr, ieig+jw ) = work( ( jw+2 )*n+jr )
1184 xmax = max( xmax, abs( vr( j, ieig ) )+
1185 $ abs( vr( j, ieig+1 ) ) )
1189 xmax = max( xmax, abs( vr( j, ieig ) ) )
1193 IF( xmax.GT.safmin )
THEN
1195 DO 490 jw = 0, nw - 1
1197 vr( jr, ieig+jw ) = xscale*vr( jr, ieig+jw )
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
subroutine dtgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
DTGEVC
subroutine dlag2(A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2, WI)
DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
subroutine dlaln2(LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, LDB, WR, WI, X, LDX, SCALE, XNORM, INFO)
DLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.