231 CHARACTER howmny, side
232 INTEGER info, ldt, ldvl, ldvr, m, mm, n
236 DOUBLE PRECISION t( ldt, * ), vl( ldvl, * ), vr( ldvr, * ),
243 DOUBLE PRECISION zero, one
244 parameter ( zero = 0.0d+0, one = 1.0d+0 )
247 LOGICAL allv, bothv, leftv, over, pair, rightv, somev
248 INTEGER i, ierr, ii, ip, is, j, j1, j2, jnxt, k, ki, n2
249 DOUBLE PRECISION beta, bignum, emax, ovfl, rec, remax, scale,
250 $ smin, smlnum, ulp, unfl, vcrit, vmax, wi, wr,
263 INTRINSIC abs, max, sqrt
266 DOUBLE PRECISION x( 2, 2 )
272 bothv =
lsame( side,
'B' )
273 rightv =
lsame( side,
'R' ) .OR. bothv
274 leftv =
lsame( side,
'L' ) .OR. bothv
276 allv =
lsame( howmny,
'A' )
277 over =
lsame( howmny,
'B' )
278 somev =
lsame( howmny,
'S' )
281 IF( .NOT.rightv .AND. .NOT.leftv )
THEN
283 ELSE IF( .NOT.allv .AND. .NOT.over .AND. .NOT.somev )
THEN
285 ELSE IF( n.LT.0 )
THEN
287 ELSE IF( ldt.LT.max( 1, n ) )
THEN
289 ELSE IF( ldvl.LT.1 .OR. ( leftv .AND. ldvl.LT.n ) )
THEN
291 ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) )
THEN
305 SELECT( j ) = .false.
308 IF( t( j+1, j ).EQ.zero )
THEN
313 IF(
SELECT( j ) .OR.
SELECT( j+1 ) )
THEN
333 CALL xerbla(
'DTREVC', -info )
344 unfl =
dlamch(
'Safe minimum' )
347 ulp =
dlamch(
'Precision' )
348 smlnum = unfl*( n / ulp )
349 bignum = ( one-ulp ) / smlnum
358 work( j ) = work( j ) + abs( t( i, j ) )
381 IF( t( ki, ki-1 ).EQ.zero )
388 IF( .NOT.
SELECT( ki ) )
391 IF( .NOT.
SELECT( ki-1 ) )
401 $ wi = sqrt( abs( t( ki, ki-1 ) ) )*
402 $ sqrt( abs( t( ki-1, ki ) ) )
403 smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum )
414 work( k+n ) = -t( k, ki )
421 DO 60 j = ki - 1, 1, -1
428 IF( t( j, j-1 ).NE.zero )
THEN
438 CALL dlaln2( .false., 1, 1, smin, one, t( j, j ),
439 $ ldt, one, one, work( j+n ), n, wr,
440 $ zero, x, 2, scale, xnorm, ierr )
445 IF( xnorm.GT.one )
THEN
446 IF( work( j ).GT.bignum / xnorm )
THEN
447 x( 1, 1 ) = x( 1, 1 ) / xnorm
448 scale = scale / xnorm
455 $
CALL dscal( ki, scale, work( 1+n ), 1 )
456 work( j+n ) = x( 1, 1 )
460 CALL daxpy( j-1, -x( 1, 1 ), t( 1, j ), 1,
467 CALL dlaln2( .false., 2, 1, smin, one,
468 $ t( j-1, j-1 ), ldt, one, one,
469 $ work( j-1+n ), n, wr, zero, x, 2,
470 $ scale, xnorm, ierr )
475 IF( xnorm.GT.one )
THEN
476 beta = max( work( j-1 ), work( j ) )
477 IF( beta.GT.bignum / xnorm )
THEN
478 x( 1, 1 ) = x( 1, 1 ) / xnorm
479 x( 2, 1 ) = x( 2, 1 ) / xnorm
480 scale = scale / xnorm
487 $
CALL dscal( ki, scale, work( 1+n ), 1 )
488 work( j-1+n ) = x( 1, 1 )
489 work( j+n ) = x( 2, 1 )
493 CALL daxpy( j-2, -x( 1, 1 ), t( 1, j-1 ), 1,
495 CALL daxpy( j-2, -x( 2, 1 ), t( 1, j ), 1,
503 CALL dcopy( ki, work( 1+n ), 1, vr( 1, is ), 1 )
505 ii =
idamax( ki, vr( 1, is ), 1 )
506 remax = one / abs( vr( ii, is ) )
507 CALL dscal( ki, remax, vr( 1, is ), 1 )
514 $
CALL dgemv(
'N', n, ki-1, one, vr, ldvr,
515 $ work( 1+n ), 1, work( ki+n ),
518 ii =
idamax( n, vr( 1, ki ), 1 )
519 remax = one / abs( vr( ii, ki ) )
520 CALL dscal( n, remax, vr( 1, ki ), 1 )
531 IF( abs( t( ki-1, ki ) ).GE.abs( t( ki, ki-1 ) ) )
THEN
533 work( ki+n2 ) = wi / t( ki-1, ki )
535 work( ki-1+n ) = -wi / t( ki, ki-1 )
539 work( ki-1+n2 ) = zero
544 work( k+n ) = -work( ki-1+n )*t( k, ki-1 )
545 work( k+n2 ) = -work( ki+n2 )*t( k, ki )
552 DO 90 j = ki - 2, 1, -1
559 IF( t( j, j-1 ).NE.zero )
THEN
569 CALL dlaln2( .false., 1, 2, smin, one, t( j, j ),
570 $ ldt, one, one, work( j+n ), n, wr, wi,
571 $ x, 2, scale, xnorm, ierr )
576 IF( xnorm.GT.one )
THEN
577 IF( work( j ).GT.bignum / xnorm )
THEN
578 x( 1, 1 ) = x( 1, 1 ) / xnorm
579 x( 1, 2 ) = x( 1, 2 ) / xnorm
580 scale = scale / xnorm
586 IF( scale.NE.one )
THEN
587 CALL dscal( ki, scale, work( 1+n ), 1 )
588 CALL dscal( ki, scale, work( 1+n2 ), 1 )
590 work( j+n ) = x( 1, 1 )
591 work( j+n2 ) = x( 1, 2 )
595 CALL daxpy( j-1, -x( 1, 1 ), t( 1, j ), 1,
597 CALL daxpy( j-1, -x( 1, 2 ), t( 1, j ), 1,
604 CALL dlaln2( .false., 2, 2, smin, one,
605 $ t( j-1, j-1 ), ldt, one, one,
606 $ work( j-1+n ), n, wr, wi, x, 2, scale,
612 IF( xnorm.GT.one )
THEN
613 beta = max( work( j-1 ), work( j ) )
614 IF( beta.GT.bignum / xnorm )
THEN
616 x( 1, 1 ) = x( 1, 1 )*rec
617 x( 1, 2 ) = x( 1, 2 )*rec
618 x( 2, 1 ) = x( 2, 1 )*rec
619 x( 2, 2 ) = x( 2, 2 )*rec
626 IF( scale.NE.one )
THEN
627 CALL dscal( ki, scale, work( 1+n ), 1 )
628 CALL dscal( ki, scale, work( 1+n2 ), 1 )
630 work( j-1+n ) = x( 1, 1 )
631 work( j+n ) = x( 2, 1 )
632 work( j-1+n2 ) = x( 1, 2 )
633 work( j+n2 ) = x( 2, 2 )
637 CALL daxpy( j-2, -x( 1, 1 ), t( 1, j-1 ), 1,
639 CALL daxpy( j-2, -x( 2, 1 ), t( 1, j ), 1,
641 CALL daxpy( j-2, -x( 1, 2 ), t( 1, j-1 ), 1,
643 CALL daxpy( j-2, -x( 2, 2 ), t( 1, j ), 1,
651 CALL dcopy( ki, work( 1+n ), 1, vr( 1, is-1 ), 1 )
652 CALL dcopy( ki, work( 1+n2 ), 1, vr( 1, is ), 1 )
656 emax = max( emax, abs( vr( k, is-1 ) )+
657 $ abs( vr( k, is ) ) )
661 CALL dscal( ki, remax, vr( 1, is-1 ), 1 )
662 CALL dscal( ki, remax, vr( 1, is ), 1 )
672 CALL dgemv(
'N', n, ki-2, one, vr, ldvr,
673 $ work( 1+n ), 1, work( ki-1+n ),
675 CALL dgemv(
'N', n, ki-2, one, vr, ldvr,
676 $ work( 1+n2 ), 1, work( ki+n2 ),
679 CALL dscal( n, work( ki-1+n ), vr( 1, ki-1 ), 1 )
680 CALL dscal( n, work( ki+n2 ), vr( 1, ki ), 1 )
685 emax = max( emax, abs( vr( k, ki-1 ) )+
686 $ abs( vr( k, ki ) ) )
689 CALL dscal( n, remax, vr( 1, ki-1 ), 1 )
690 CALL dscal( n, remax, vr( 1, ki ), 1 )
717 IF( t( ki+1, ki ).EQ.zero )
723 IF( .NOT.
SELECT( ki ) )
732 $ wi = sqrt( abs( t( ki, ki+1 ) ) )*
733 $ sqrt( abs( t( ki+1, ki ) ) )
734 smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum )
745 work( k+n ) = -t( ki, k )
762 IF( t( j+1, j ).NE.zero )
THEN
775 IF( work( j ).GT.vcrit )
THEN
777 CALL dscal( n-ki+1, rec, work( ki+n ), 1 )
782 work( j+n ) = work( j+n ) -
783 $
ddot( j-ki-1, t( ki+1, j ), 1,
784 $ work( ki+1+n ), 1 )
788 CALL dlaln2( .false., 1, 1, smin, one, t( j, j ),
789 $ ldt, one, one, work( j+n ), n, wr,
790 $ zero, x, 2, scale, xnorm, ierr )
795 $
CALL dscal( n-ki+1, scale, work( ki+n ), 1 )
796 work( j+n ) = x( 1, 1 )
797 vmax = max( abs( work( j+n ) ), vmax )
798 vcrit = bignum / vmax
807 beta = max( work( j ), work( j+1 ) )
808 IF( beta.GT.vcrit )
THEN
810 CALL dscal( n-ki+1, rec, work( ki+n ), 1 )
815 work( j+n ) = work( j+n ) -
816 $
ddot( j-ki-1, t( ki+1, j ), 1,
817 $ work( ki+1+n ), 1 )
819 work( j+1+n ) = work( j+1+n ) -
820 $
ddot( j-ki-1, t( ki+1, j+1 ), 1,
821 $ work( ki+1+n ), 1 )
827 CALL dlaln2( .true., 2, 1, smin, one, t( j, j ),
828 $ ldt, one, one, work( j+n ), n, wr,
829 $ zero, x, 2, scale, xnorm, ierr )
834 $
CALL dscal( n-ki+1, scale, work( ki+n ), 1 )
835 work( j+n ) = x( 1, 1 )
836 work( j+1+n ) = x( 2, 1 )
838 vmax = max( abs( work( j+n ) ),
839 $ abs( work( j+1+n ) ), vmax )
840 vcrit = bignum / vmax
848 CALL dcopy( n-ki+1, work( ki+n ), 1, vl( ki, is ), 1 )
850 ii =
idamax( n-ki+1, vl( ki, is ), 1 ) + ki - 1
851 remax = one / abs( vl( ii, is ) )
852 CALL dscal( n-ki+1, remax, vl( ki, is ), 1 )
861 $
CALL dgemv(
'N', n, n-ki, one, vl( 1, ki+1 ), ldvl,
862 $ work( ki+1+n ), 1, work( ki+n ),
865 ii =
idamax( n, vl( 1, ki ), 1 )
866 remax = one / abs( vl( ii, ki ) )
867 CALL dscal( n, remax, vl( 1, ki ), 1 )
879 IF( abs( t( ki, ki+1 ) ).GE.abs( t( ki+1, ki ) ) )
THEN
880 work( ki+n ) = wi / t( ki, ki+1 )
881 work( ki+1+n2 ) = one
884 work( ki+1+n2 ) = -wi / t( ki+1, ki )
886 work( ki+1+n ) = zero
892 work( k+n ) = -work( ki+n )*t( ki, k )
893 work( k+n2 ) = -work( ki+1+n2 )*t( ki+1, k )
910 IF( t( j+1, j ).NE.zero )
THEN
923 IF( work( j ).GT.vcrit )
THEN
925 CALL dscal( n-ki+1, rec, work( ki+n ), 1 )
926 CALL dscal( n-ki+1, rec, work( ki+n2 ), 1 )
931 work( j+n ) = work( j+n ) -
932 $
ddot( j-ki-2, t( ki+2, j ), 1,
933 $ work( ki+2+n ), 1 )
934 work( j+n2 ) = work( j+n2 ) -
935 $
ddot( j-ki-2, t( ki+2, j ), 1,
936 $ work( ki+2+n2 ), 1 )
940 CALL dlaln2( .false., 1, 2, smin, one, t( j, j ),
941 $ ldt, one, one, work( j+n ), n, wr,
942 $ -wi, x, 2, scale, xnorm, ierr )
946 IF( scale.NE.one )
THEN
947 CALL dscal( n-ki+1, scale, work( ki+n ), 1 )
948 CALL dscal( n-ki+1, scale, work( ki+n2 ), 1 )
950 work( j+n ) = x( 1, 1 )
951 work( j+n2 ) = x( 1, 2 )
952 vmax = max( abs( work( j+n ) ),
953 $ abs( work( j+n2 ) ), vmax )
954 vcrit = bignum / vmax
963 beta = max( work( j ), work( j+1 ) )
964 IF( beta.GT.vcrit )
THEN
966 CALL dscal( n-ki+1, rec, work( ki+n ), 1 )
967 CALL dscal( n-ki+1, rec, work( ki+n2 ), 1 )
972 work( j+n ) = work( j+n ) -
973 $
ddot( j-ki-2, t( ki+2, j ), 1,
974 $ work( ki+2+n ), 1 )
976 work( j+n2 ) = work( j+n2 ) -
977 $
ddot( j-ki-2, t( ki+2, j ), 1,
978 $ work( ki+2+n2 ), 1 )
980 work( j+1+n ) = work( j+1+n ) -
981 $
ddot( j-ki-2, t( ki+2, j+1 ), 1,
982 $ work( ki+2+n ), 1 )
984 work( j+1+n2 ) = work( j+1+n2 ) -
985 $
ddot( j-ki-2, t( ki+2, j+1 ), 1,
986 $ work( ki+2+n2 ), 1 )
992 CALL dlaln2( .true., 2, 2, smin, one, t( j, j ),
993 $ ldt, one, one, work( j+n ), n, wr,
994 $ -wi, x, 2, scale, xnorm, ierr )
998 IF( scale.NE.one )
THEN
999 CALL dscal( n-ki+1, scale, work( ki+n ), 1 )
1000 CALL dscal( n-ki+1, scale, work( ki+n2 ), 1 )
1002 work( j+n ) = x( 1, 1 )
1003 work( j+n2 ) = x( 1, 2 )
1004 work( j+1+n ) = x( 2, 1 )
1005 work( j+1+n2 ) = x( 2, 2 )
1006 vmax = max( abs( x( 1, 1 ) ), abs( x( 1, 2 ) ),
1007 $ abs( x( 2, 1 ) ), abs( x( 2, 2 ) ), vmax )
1008 vcrit = bignum / vmax
1015 IF( .NOT.over )
THEN
1016 CALL dcopy( n-ki+1, work( ki+n ), 1, vl( ki, is ), 1 )
1017 CALL dcopy( n-ki+1, work( ki+n2 ), 1, vl( ki, is+1 ),
1022 emax = max( emax, abs( vl( k, is ) )+
1023 $ abs( vl( k, is+1 ) ) )
1026 CALL dscal( n-ki+1, remax, vl( ki, is ), 1 )
1027 CALL dscal( n-ki+1, remax, vl( ki, is+1 ), 1 )
1029 DO 230 k = 1, ki - 1
1031 vl( k, is+1 ) = zero
1034 IF( ki.LT.n-1 )
THEN
1035 CALL dgemv(
'N', n, n-ki-1, one, vl( 1, ki+2 ),
1036 $ ldvl, work( ki+2+n ), 1, work( ki+n ),
1038 CALL dgemv(
'N', n, n-ki-1, one, vl( 1, ki+2 ),
1039 $ ldvl, work( ki+2+n2 ), 1,
1040 $ work( ki+1+n2 ), vl( 1, ki+1 ), 1 )
1042 CALL dscal( n, work( ki+n ), vl( 1, ki ), 1 )
1043 CALL dscal( n, work( ki+1+n2 ), vl( 1, ki+1 ), 1 )
1048 emax = max( emax, abs( vl( k, ki ) )+
1049 $ abs( vl( k, ki+1 ) ) )
1052 CALL dscal( n, remax, vl( 1, ki ), 1 )
1053 CALL dscal( n, remax, vl( 1, ki+1 ), 1 )
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
integer function idamax(N, DX, INCX)
IDAMAX
double precision function dlamch(CMACH)
DLAMCH
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlabad(SMALL, LARGE)
DLABAD
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.
subroutine dscal(N, DA, DX, INCX)
DSCAL
logical function lsame(CA, CB)
LSAME