228 CHARACTER HOWMNY, SIDE
229 INTEGER INFO, LDT, LDVL, LDVR, M, MM, N
233 DOUBLE PRECISION T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
240 DOUBLE PRECISION ZERO, ONE
241 parameter( zero = 0.0d+0, one = 1.0d+0 )
244 LOGICAL ALLV, BOTHV, LEFTV, OVER, PAIR, RIGHTV, SOMEV
245 INTEGER I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI, N2
246 DOUBLE PRECISION BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
247 $ SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR,
253 DOUBLE PRECISION DDOT, DLAMCH
261 INTRINSIC abs, max, sqrt
264 DOUBLE PRECISION X( 2, 2 )
270 bothv =
lsame( side,
'B' )
271 rightv =
lsame( side,
'R' ) .OR. bothv
272 leftv =
lsame( side,
'L' ) .OR. bothv
274 allv =
lsame( howmny,
'A' )
275 over =
lsame( howmny,
'B' )
276 somev =
lsame( howmny,
'S' )
279 IF( .NOT.rightv .AND. .NOT.leftv )
THEN
281 ELSE IF( .NOT.allv .AND. .NOT.over .AND. .NOT.somev )
THEN
283 ELSE IF( n.LT.0 )
THEN
285 ELSE IF( ldt.LT.max( 1, n ) )
THEN
287 ELSE IF( ldvl.LT.1 .OR. ( leftv .AND. ldvl.LT.n ) )
THEN
289 ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) )
THEN
303 SELECT( j ) = .false.
306 IF( t( j+1, j ).EQ.zero )
THEN
311 IF(
SELECT( j ) .OR.
SELECT( j+1 ) )
THEN
331 CALL xerbla(
'DTREVC', -info )
342 unfl =
dlamch(
'Safe minimum' )
345 ulp =
dlamch(
'Precision' )
346 smlnum = unfl*( n / ulp )
347 bignum = ( one-ulp ) / smlnum
356 work( j ) = work( j ) + abs( t( i, j ) )
379 IF( t( ki, ki-1 ).EQ.zero )
386 IF( .NOT.
SELECT( ki ) )
389 IF( .NOT.
SELECT( ki-1 ) )
399 $ wi = sqrt( abs( t( ki, ki-1 ) ) )*
400 $ sqrt( abs( t( ki-1, ki ) ) )
401 smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum )
412 work( k+n ) = -t( k, ki )
419 DO 60 j = ki - 1, 1, -1
426 IF( t( j, j-1 ).NE.zero )
THEN
436 CALL dlaln2( .false., 1, 1, smin, one, t( j, j ),
437 $ ldt, one, one, work( j+n ), n, wr,
438 $ zero, x, 2, scale, xnorm, ierr )
443 IF( xnorm.GT.one )
THEN
444 IF( work( j ).GT.bignum / xnorm )
THEN
445 x( 1, 1 ) = x( 1, 1 ) / xnorm
446 scale = scale / xnorm
453 $
CALL dscal( ki, scale, work( 1+n ), 1 )
454 work( j+n ) = x( 1, 1 )
458 CALL daxpy( j-1, -x( 1, 1 ), t( 1, j ), 1,
465 CALL dlaln2( .false., 2, 1, smin, one,
466 $ t( j-1, j-1 ), ldt, one, one,
467 $ work( j-1+n ), n, wr, zero, x, 2,
468 $ scale, xnorm, ierr )
473 IF( xnorm.GT.one )
THEN
474 beta = max( work( j-1 ), work( j ) )
475 IF( beta.GT.bignum / xnorm )
THEN
476 x( 1, 1 ) = x( 1, 1 ) / xnorm
477 x( 2, 1 ) = x( 2, 1 ) / xnorm
478 scale = scale / xnorm
485 $
CALL dscal( ki, scale, work( 1+n ), 1 )
486 work( j-1+n ) = x( 1, 1 )
487 work( j+n ) = x( 2, 1 )
491 CALL daxpy( j-2, -x( 1, 1 ), t( 1, j-1 ), 1,
493 CALL daxpy( j-2, -x( 2, 1 ), t( 1, j ), 1,
501 CALL dcopy( ki, work( 1+n ), 1, vr( 1, is ), 1 )
503 ii =
idamax( ki, vr( 1, is ), 1 )
504 remax = one / abs( vr( ii, is ) )
505 CALL dscal( ki, remax, vr( 1, is ), 1 )
512 $
CALL dgemv(
'N', n, ki-1, one, vr, ldvr,
513 $ work( 1+n ), 1, work( ki+n ),
516 ii =
idamax( n, vr( 1, ki ), 1 )
517 remax = one / abs( vr( ii, ki ) )
518 CALL dscal( n, remax, vr( 1, ki ), 1 )
529 IF( abs( t( ki-1, ki ) ).GE.abs( t( ki, ki-1 ) ) )
THEN
531 work( ki+n2 ) = wi / t( ki-1, ki )
533 work( ki-1+n ) = -wi / t( ki, ki-1 )
537 work( ki-1+n2 ) = zero
542 work( k+n ) = -work( ki-1+n )*t( k, ki-1 )
543 work( k+n2 ) = -work( ki+n2 )*t( k, ki )
550 DO 90 j = ki - 2, 1, -1
557 IF( t( j, j-1 ).NE.zero )
THEN
567 CALL dlaln2( .false., 1, 2, smin, one, t( j, j ),
568 $ ldt, one, one, work( j+n ), n, wr, wi,
569 $ x, 2, scale, xnorm, ierr )
574 IF( xnorm.GT.one )
THEN
575 IF( work( j ).GT.bignum / xnorm )
THEN
576 x( 1, 1 ) = x( 1, 1 ) / xnorm
577 x( 1, 2 ) = x( 1, 2 ) / xnorm
578 scale = scale / xnorm
584 IF( scale.NE.one )
THEN
585 CALL dscal( ki, scale, work( 1+n ), 1 )
586 CALL dscal( ki, scale, work( 1+n2 ), 1 )
588 work( j+n ) = x( 1, 1 )
589 work( j+n2 ) = x( 1, 2 )
593 CALL daxpy( j-1, -x( 1, 1 ), t( 1, j ), 1,
595 CALL daxpy( j-1, -x( 1, 2 ), t( 1, j ), 1,
602 CALL dlaln2( .false., 2, 2, smin, one,
603 $ t( j-1, j-1 ), ldt, one, one,
604 $ work( j-1+n ), n, wr, wi, x, 2, scale,
610 IF( xnorm.GT.one )
THEN
611 beta = max( work( j-1 ), work( j ) )
612 IF( beta.GT.bignum / xnorm )
THEN
614 x( 1, 1 ) = x( 1, 1 )*rec
615 x( 1, 2 ) = x( 1, 2 )*rec
616 x( 2, 1 ) = x( 2, 1 )*rec
617 x( 2, 2 ) = x( 2, 2 )*rec
624 IF( scale.NE.one )
THEN
625 CALL dscal( ki, scale, work( 1+n ), 1 )
626 CALL dscal( ki, scale, work( 1+n2 ), 1 )
628 work( j-1+n ) = x( 1, 1 )
629 work( j+n ) = x( 2, 1 )
630 work( j-1+n2 ) = x( 1, 2 )
631 work( j+n2 ) = x( 2, 2 )
635 CALL daxpy( j-2, -x( 1, 1 ), t( 1, j-1 ), 1,
637 CALL daxpy( j-2, -x( 2, 1 ), t( 1, j ), 1,
639 CALL daxpy( j-2, -x( 1, 2 ), t( 1, j-1 ), 1,
641 CALL daxpy( j-2, -x( 2, 2 ), t( 1, j ), 1,
649 CALL dcopy( ki, work( 1+n ), 1, vr( 1, is-1 ), 1 )
650 CALL dcopy( ki, work( 1+n2 ), 1, vr( 1, is ), 1 )
654 emax = max( emax, abs( vr( k, is-1 ) )+
655 $ abs( vr( k, is ) ) )
659 CALL dscal( ki, remax, vr( 1, is-1 ), 1 )
660 CALL dscal( ki, remax, vr( 1, is ), 1 )
670 CALL dgemv(
'N', n, ki-2, one, vr, ldvr,
671 $ work( 1+n ), 1, work( ki-1+n ),
673 CALL dgemv(
'N', n, ki-2, one, vr, ldvr,
674 $ work( 1+n2 ), 1, work( ki+n2 ),
677 CALL dscal( n, work( ki-1+n ), vr( 1, ki-1 ), 1 )
678 CALL dscal( n, work( ki+n2 ), vr( 1, ki ), 1 )
683 emax = max( emax, abs( vr( k, ki-1 ) )+
684 $ abs( vr( k, ki ) ) )
687 CALL dscal( n, remax, vr( 1, ki-1 ), 1 )
688 CALL dscal( n, remax, vr( 1, ki ), 1 )
715 IF( t( ki+1, ki ).EQ.zero )
721 IF( .NOT.
SELECT( ki ) )
730 $ wi = sqrt( abs( t( ki, ki+1 ) ) )*
731 $ sqrt( abs( t( ki+1, ki ) ) )
732 smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum )
743 work( k+n ) = -t( ki, k )
760 IF( t( j+1, j ).NE.zero )
THEN
773 IF( work( j ).GT.vcrit )
THEN
775 CALL dscal( n-ki+1, rec, work( ki+n ), 1 )
780 work( j+n ) = work( j+n ) -
781 $
ddot( j-ki-1, t( ki+1, j ), 1,
782 $ work( ki+1+n ), 1 )
786 CALL dlaln2( .false., 1, 1, smin, one, t( j, j ),
787 $ ldt, one, one, work( j+n ), n, wr,
788 $ zero, x, 2, scale, xnorm, ierr )
793 $
CALL dscal( n-ki+1, scale, work( ki+n ), 1 )
794 work( j+n ) = x( 1, 1 )
795 vmax = max( abs( work( j+n ) ), vmax )
796 vcrit = bignum / vmax
805 beta = max( work( j ), work( j+1 ) )
806 IF( beta.GT.vcrit )
THEN
808 CALL dscal( n-ki+1, rec, work( ki+n ), 1 )
813 work( j+n ) = work( j+n ) -
814 $
ddot( j-ki-1, t( ki+1, j ), 1,
815 $ work( ki+1+n ), 1 )
817 work( j+1+n ) = work( j+1+n ) -
818 $
ddot( j-ki-1, t( ki+1, j+1 ), 1,
819 $ work( ki+1+n ), 1 )
825 CALL dlaln2( .true., 2, 1, smin, one, t( j, j ),
826 $ ldt, one, one, work( j+n ), n, wr,
827 $ zero, x, 2, scale, xnorm, ierr )
832 $
CALL dscal( n-ki+1, scale, work( ki+n ), 1 )
833 work( j+n ) = x( 1, 1 )
834 work( j+1+n ) = x( 2, 1 )
836 vmax = max( abs( work( j+n ) ),
837 $ abs( work( j+1+n ) ), vmax )
838 vcrit = bignum / vmax
846 CALL dcopy( n-ki+1, work( ki+n ), 1, vl( ki, is ), 1 )
848 ii =
idamax( n-ki+1, vl( ki, is ), 1 ) + ki - 1
849 remax = one / abs( vl( ii, is ) )
850 CALL dscal( n-ki+1, remax, vl( ki, is ), 1 )
859 $
CALL dgemv(
'N', n, n-ki, one, vl( 1, ki+1 ), ldvl,
860 $ work( ki+1+n ), 1, work( ki+n ),
863 ii =
idamax( n, vl( 1, ki ), 1 )
864 remax = one / abs( vl( ii, ki ) )
865 CALL dscal( n, remax, vl( 1, ki ), 1 )
877 IF( abs( t( ki, ki+1 ) ).GE.abs( t( ki+1, ki ) ) )
THEN
878 work( ki+n ) = wi / t( ki, ki+1 )
879 work( ki+1+n2 ) = one
882 work( ki+1+n2 ) = -wi / t( ki+1, ki )
884 work( ki+1+n ) = zero
890 work( k+n ) = -work( ki+n )*t( ki, k )
891 work( k+n2 ) = -work( ki+1+n2 )*t( ki+1, k )
908 IF( t( j+1, j ).NE.zero )
THEN
921 IF( work( j ).GT.vcrit )
THEN
923 CALL dscal( n-ki+1, rec, work( ki+n ), 1 )
924 CALL dscal( n-ki+1, rec, work( ki+n2 ), 1 )
929 work( j+n ) = work( j+n ) -
930 $
ddot( j-ki-2, t( ki+2, j ), 1,
931 $ work( ki+2+n ), 1 )
932 work( j+n2 ) = work( j+n2 ) -
933 $
ddot( j-ki-2, t( ki+2, j ), 1,
934 $ work( ki+2+n2 ), 1 )
938 CALL dlaln2( .false., 1, 2, smin, one, t( j, j ),
939 $ ldt, one, one, work( j+n ), n, wr,
940 $ -wi, x, 2, scale, xnorm, ierr )
944 IF( scale.NE.one )
THEN
945 CALL dscal( n-ki+1, scale, work( ki+n ), 1 )
946 CALL dscal( n-ki+1, scale, work( ki+n2 ), 1 )
948 work( j+n ) = x( 1, 1 )
949 work( j+n2 ) = x( 1, 2 )
950 vmax = max( abs( work( j+n ) ),
951 $ abs( work( j+n2 ) ), vmax )
952 vcrit = bignum / vmax
961 beta = max( work( j ), work( j+1 ) )
962 IF( beta.GT.vcrit )
THEN
964 CALL dscal( n-ki+1, rec, work( ki+n ), 1 )
965 CALL dscal( n-ki+1, rec, work( ki+n2 ), 1 )
970 work( j+n ) = work( j+n ) -
971 $
ddot( j-ki-2, t( ki+2, j ), 1,
972 $ work( ki+2+n ), 1 )
974 work( j+n2 ) = work( j+n2 ) -
975 $
ddot( j-ki-2, t( ki+2, j ), 1,
976 $ work( ki+2+n2 ), 1 )
978 work( j+1+n ) = work( j+1+n ) -
979 $
ddot( j-ki-2, t( ki+2, j+1 ), 1,
980 $ work( ki+2+n ), 1 )
982 work( j+1+n2 ) = work( j+1+n2 ) -
983 $
ddot( j-ki-2, t( ki+2, j+1 ), 1,
984 $ work( ki+2+n2 ), 1 )
990 CALL dlaln2( .true., 2, 2, smin, one, t( j, j ),
991 $ ldt, one, one, work( j+n ), n, wr,
992 $ -wi, x, 2, scale, xnorm, ierr )
996 IF( scale.NE.one )
THEN
997 CALL dscal( n-ki+1, scale, work( ki+n ), 1 )
998 CALL dscal( n-ki+1, scale, work( ki+n2 ), 1 )
1000 work( j+n ) = x( 1, 1 )
1001 work( j+n2 ) = x( 1, 2 )
1002 work( j+1+n ) = x( 2, 1 )
1003 work( j+1+n2 ) = x( 2, 2 )
1004 vmax = max( abs( x( 1, 1 ) ), abs( x( 1, 2 ) ),
1005 $ abs( x( 2, 1 ) ), abs( x( 2, 2 ) ), vmax )
1006 vcrit = bignum / vmax
1013 IF( .NOT.over )
THEN
1014 CALL dcopy( n-ki+1, work( ki+n ), 1, vl( ki, is ), 1 )
1015 CALL dcopy( n-ki+1, work( ki+n2 ), 1, vl( ki, is+1 ),
1020 emax = max( emax, abs( vl( k, is ) )+
1021 $ abs( vl( k, is+1 ) ) )
1024 CALL dscal( n-ki+1, remax, vl( ki, is ), 1 )
1025 CALL dscal( n-ki+1, remax, vl( ki, is+1 ), 1 )
1027 DO 230 k = 1, ki - 1
1029 vl( k, is+1 ) = zero
1032 IF( ki.LT.n-1 )
THEN
1033 CALL dgemv(
'N', n, n-ki-1, one, vl( 1, ki+2 ),
1034 $ ldvl, work( ki+2+n ), 1, work( ki+n ),
1036 CALL dgemv(
'N', n, n-ki-1, one, vl( 1, ki+2 ),
1037 $ ldvl, work( ki+2+n2 ), 1,
1038 $ work( ki+1+n2 ), vl( 1, ki+1 ), 1 )
1040 CALL dscal( n, work( ki+n ), vl( 1, ki ), 1 )
1041 CALL dscal( n, work( ki+1+n2 ), vl( 1, ki+1 ), 1 )
1046 emax = max( emax, abs( vl( k, ki ) )+
1047 $ abs( vl( k, ki+1 ) ) )
1050 CALL dscal( n, remax, vl( 1, ki ), 1 )
1051 CALL dscal( n, remax, vl( 1, ki+1 ), 1 )
double precision function dlamch(CMACH)
DLAMCH
subroutine dlabad(SMALL, LARGE)
DLABAD
integer function idamax(N, DX, INCX)
IDAMAX
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
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.