231 CHARACTER howmny, side
232 INTEGER info, ldt, ldvl, ldvr, m, mm, n
236 REAL t( ldt, * ), vl( ldvl, * ), vr( ldvr, * ),
244 parameter ( zero = 0.0e+0, one = 1.0e+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 REAL beta, bignum, emax, ovfl, rec, remax, scale,
250 $ smin, smlnum, ulp, unfl, vcrit, vmax, wi, wr,
264 INTRINSIC abs, max, sqrt
273 bothv =
lsame( side,
'B' )
274 rightv =
lsame( side,
'R' ) .OR. bothv
275 leftv =
lsame( side,
'L' ) .OR. bothv
277 allv =
lsame( howmny,
'A' )
278 over =
lsame( howmny,
'B' )
279 somev =
lsame( howmny,
'S' )
282 IF( .NOT.rightv .AND. .NOT.leftv )
THEN
284 ELSE IF( .NOT.allv .AND. .NOT.over .AND. .NOT.somev )
THEN
286 ELSE IF( n.LT.0 )
THEN
288 ELSE IF( ldt.LT.max( 1, n ) )
THEN
290 ELSE IF( ldvl.LT.1 .OR. ( leftv .AND. ldvl.LT.n ) )
THEN
292 ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) )
THEN
306 SELECT( j ) = .false.
309 IF( t( j+1, j ).EQ.zero )
THEN
314 IF(
SELECT( j ) .OR.
SELECT( j+1 ) )
THEN
334 CALL xerbla(
'STREVC', -info )
345 unfl =
slamch(
'Safe minimum' )
348 ulp =
slamch(
'Precision' )
349 smlnum = unfl*( n / ulp )
350 bignum = ( one-ulp ) / smlnum
359 work( j ) = work( j ) + abs( t( i, j ) )
382 IF( t( ki, ki-1 ).EQ.zero )
389 IF( .NOT.
SELECT( ki ) )
392 IF( .NOT.
SELECT( ki-1 ) )
402 $ wi = sqrt( abs( t( ki, ki-1 ) ) )*
403 $ sqrt( abs( t( ki-1, ki ) ) )
404 smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum )
415 work( k+n ) = -t( k, ki )
422 DO 60 j = ki - 1, 1, -1
429 IF( t( j, j-1 ).NE.zero )
THEN
439 CALL slaln2( .false., 1, 1, smin, one, t( j, j ),
440 $ ldt, one, one, work( j+n ), n, wr,
441 $ zero, x, 2, scale, xnorm, ierr )
446 IF( xnorm.GT.one )
THEN
447 IF( work( j ).GT.bignum / xnorm )
THEN
448 x( 1, 1 ) = x( 1, 1 ) / xnorm
449 scale = scale / xnorm
456 $
CALL sscal( ki, scale, work( 1+n ), 1 )
457 work( j+n ) = x( 1, 1 )
461 CALL saxpy( j-1, -x( 1, 1 ), t( 1, j ), 1,
468 CALL slaln2( .false., 2, 1, smin, one,
469 $ t( j-1, j-1 ), ldt, one, one,
470 $ work( j-1+n ), n, wr, zero, x, 2,
471 $ scale, xnorm, ierr )
476 IF( xnorm.GT.one )
THEN
477 beta = max( work( j-1 ), work( j ) )
478 IF( beta.GT.bignum / xnorm )
THEN
479 x( 1, 1 ) = x( 1, 1 ) / xnorm
480 x( 2, 1 ) = x( 2, 1 ) / xnorm
481 scale = scale / xnorm
488 $
CALL sscal( ki, scale, work( 1+n ), 1 )
489 work( j-1+n ) = x( 1, 1 )
490 work( j+n ) = x( 2, 1 )
494 CALL saxpy( j-2, -x( 1, 1 ), t( 1, j-1 ), 1,
496 CALL saxpy( j-2, -x( 2, 1 ), t( 1, j ), 1,
504 CALL scopy( ki, work( 1+n ), 1, vr( 1, is ), 1 )
506 ii =
isamax( ki, vr( 1, is ), 1 )
507 remax = one / abs( vr( ii, is ) )
508 CALL sscal( ki, remax, vr( 1, is ), 1 )
515 $
CALL sgemv(
'N', n, ki-1, one, vr, ldvr,
516 $ work( 1+n ), 1, work( ki+n ),
519 ii =
isamax( n, vr( 1, ki ), 1 )
520 remax = one / abs( vr( ii, ki ) )
521 CALL sscal( n, remax, vr( 1, ki ), 1 )
532 IF( abs( t( ki-1, ki ) ).GE.abs( t( ki, ki-1 ) ) )
THEN
534 work( ki+n2 ) = wi / t( ki-1, ki )
536 work( ki-1+n ) = -wi / t( ki, ki-1 )
540 work( ki-1+n2 ) = zero
545 work( k+n ) = -work( ki-1+n )*t( k, ki-1 )
546 work( k+n2 ) = -work( ki+n2 )*t( k, ki )
553 DO 90 j = ki - 2, 1, -1
560 IF( t( j, j-1 ).NE.zero )
THEN
570 CALL slaln2( .false., 1, 2, smin, one, t( j, j ),
571 $ ldt, one, one, work( j+n ), n, wr, wi,
572 $ x, 2, scale, xnorm, ierr )
577 IF( xnorm.GT.one )
THEN
578 IF( work( j ).GT.bignum / xnorm )
THEN
579 x( 1, 1 ) = x( 1, 1 ) / xnorm
580 x( 1, 2 ) = x( 1, 2 ) / xnorm
581 scale = scale / xnorm
587 IF( scale.NE.one )
THEN
588 CALL sscal( ki, scale, work( 1+n ), 1 )
589 CALL sscal( ki, scale, work( 1+n2 ), 1 )
591 work( j+n ) = x( 1, 1 )
592 work( j+n2 ) = x( 1, 2 )
596 CALL saxpy( j-1, -x( 1, 1 ), t( 1, j ), 1,
598 CALL saxpy( j-1, -x( 1, 2 ), t( 1, j ), 1,
605 CALL slaln2( .false., 2, 2, smin, one,
606 $ t( j-1, j-1 ), ldt, one, one,
607 $ work( j-1+n ), n, wr, wi, x, 2, scale,
613 IF( xnorm.GT.one )
THEN
614 beta = max( work( j-1 ), work( j ) )
615 IF( beta.GT.bignum / xnorm )
THEN
617 x( 1, 1 ) = x( 1, 1 )*rec
618 x( 1, 2 ) = x( 1, 2 )*rec
619 x( 2, 1 ) = x( 2, 1 )*rec
620 x( 2, 2 ) = x( 2, 2 )*rec
627 IF( scale.NE.one )
THEN
628 CALL sscal( ki, scale, work( 1+n ), 1 )
629 CALL sscal( ki, scale, work( 1+n2 ), 1 )
631 work( j-1+n ) = x( 1, 1 )
632 work( j+n ) = x( 2, 1 )
633 work( j-1+n2 ) = x( 1, 2 )
634 work( j+n2 ) = x( 2, 2 )
638 CALL saxpy( j-2, -x( 1, 1 ), t( 1, j-1 ), 1,
640 CALL saxpy( j-2, -x( 2, 1 ), t( 1, j ), 1,
642 CALL saxpy( j-2, -x( 1, 2 ), t( 1, j-1 ), 1,
644 CALL saxpy( j-2, -x( 2, 2 ), t( 1, j ), 1,
652 CALL scopy( ki, work( 1+n ), 1, vr( 1, is-1 ), 1 )
653 CALL scopy( ki, work( 1+n2 ), 1, vr( 1, is ), 1 )
657 emax = max( emax, abs( vr( k, is-1 ) )+
658 $ abs( vr( k, is ) ) )
662 CALL sscal( ki, remax, vr( 1, is-1 ), 1 )
663 CALL sscal( ki, remax, vr( 1, is ), 1 )
673 CALL sgemv(
'N', n, ki-2, one, vr, ldvr,
674 $ work( 1+n ), 1, work( ki-1+n ),
676 CALL sgemv(
'N', n, ki-2, one, vr, ldvr,
677 $ work( 1+n2 ), 1, work( ki+n2 ),
680 CALL sscal( n, work( ki-1+n ), vr( 1, ki-1 ), 1 )
681 CALL sscal( n, work( ki+n2 ), vr( 1, ki ), 1 )
686 emax = max( emax, abs( vr( k, ki-1 ) )+
687 $ abs( vr( k, ki ) ) )
690 CALL sscal( n, remax, vr( 1, ki-1 ), 1 )
691 CALL sscal( n, remax, vr( 1, ki ), 1 )
718 IF( t( ki+1, ki ).EQ.zero )
724 IF( .NOT.
SELECT( ki ) )
733 $ wi = sqrt( abs( t( ki, ki+1 ) ) )*
734 $ sqrt( abs( t( ki+1, ki ) ) )
735 smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum )
746 work( k+n ) = -t( ki, k )
763 IF( t( j+1, j ).NE.zero )
THEN
776 IF( work( j ).GT.vcrit )
THEN
778 CALL sscal( n-ki+1, rec, work( ki+n ), 1 )
783 work( j+n ) = work( j+n ) -
784 $
sdot( j-ki-1, t( ki+1, j ), 1,
785 $ work( ki+1+n ), 1 )
789 CALL slaln2( .false., 1, 1, smin, one, t( j, j ),
790 $ ldt, one, one, work( j+n ), n, wr,
791 $ zero, x, 2, scale, xnorm, ierr )
796 $
CALL sscal( n-ki+1, scale, work( ki+n ), 1 )
797 work( j+n ) = x( 1, 1 )
798 vmax = max( abs( work( j+n ) ), vmax )
799 vcrit = bignum / vmax
808 beta = max( work( j ), work( j+1 ) )
809 IF( beta.GT.vcrit )
THEN
811 CALL sscal( n-ki+1, rec, work( ki+n ), 1 )
816 work( j+n ) = work( j+n ) -
817 $
sdot( j-ki-1, t( ki+1, j ), 1,
818 $ work( ki+1+n ), 1 )
820 work( j+1+n ) = work( j+1+n ) -
821 $
sdot( j-ki-1, t( ki+1, j+1 ), 1,
822 $ work( ki+1+n ), 1 )
828 CALL slaln2( .true., 2, 1, smin, one, t( j, j ),
829 $ ldt, one, one, work( j+n ), n, wr,
830 $ zero, x, 2, scale, xnorm, ierr )
835 $
CALL sscal( n-ki+1, scale, work( ki+n ), 1 )
836 work( j+n ) = x( 1, 1 )
837 work( j+1+n ) = x( 2, 1 )
839 vmax = max( abs( work( j+n ) ),
840 $ abs( work( j+1+n ) ), vmax )
841 vcrit = bignum / vmax
849 CALL scopy( n-ki+1, work( ki+n ), 1, vl( ki, is ), 1 )
851 ii =
isamax( n-ki+1, vl( ki, is ), 1 ) + ki - 1
852 remax = one / abs( vl( ii, is ) )
853 CALL sscal( n-ki+1, remax, vl( ki, is ), 1 )
862 $
CALL sgemv(
'N', n, n-ki, one, vl( 1, ki+1 ), ldvl,
863 $ work( ki+1+n ), 1, work( ki+n ),
866 ii =
isamax( n, vl( 1, ki ), 1 )
867 remax = one / abs( vl( ii, ki ) )
868 CALL sscal( n, remax, vl( 1, ki ), 1 )
880 IF( abs( t( ki, ki+1 ) ).GE.abs( t( ki+1, ki ) ) )
THEN
881 work( ki+n ) = wi / t( ki, ki+1 )
882 work( ki+1+n2 ) = one
885 work( ki+1+n2 ) = -wi / t( ki+1, ki )
887 work( ki+1+n ) = zero
893 work( k+n ) = -work( ki+n )*t( ki, k )
894 work( k+n2 ) = -work( ki+1+n2 )*t( ki+1, k )
911 IF( t( j+1, j ).NE.zero )
THEN
924 IF( work( j ).GT.vcrit )
THEN
926 CALL sscal( n-ki+1, rec, work( ki+n ), 1 )
927 CALL sscal( n-ki+1, rec, work( ki+n2 ), 1 )
932 work( j+n ) = work( j+n ) -
933 $
sdot( j-ki-2, t( ki+2, j ), 1,
934 $ work( ki+2+n ), 1 )
935 work( j+n2 ) = work( j+n2 ) -
936 $
sdot( j-ki-2, t( ki+2, j ), 1,
937 $ work( ki+2+n2 ), 1 )
941 CALL slaln2( .false., 1, 2, smin, one, t( j, j ),
942 $ ldt, one, one, work( j+n ), n, wr,
943 $ -wi, x, 2, scale, xnorm, ierr )
947 IF( scale.NE.one )
THEN
948 CALL sscal( n-ki+1, scale, work( ki+n ), 1 )
949 CALL sscal( n-ki+1, scale, work( ki+n2 ), 1 )
951 work( j+n ) = x( 1, 1 )
952 work( j+n2 ) = x( 1, 2 )
953 vmax = max( abs( work( j+n ) ),
954 $ abs( work( j+n2 ) ), vmax )
955 vcrit = bignum / vmax
964 beta = max( work( j ), work( j+1 ) )
965 IF( beta.GT.vcrit )
THEN
967 CALL sscal( n-ki+1, rec, work( ki+n ), 1 )
968 CALL sscal( n-ki+1, rec, work( ki+n2 ), 1 )
973 work( j+n ) = work( j+n ) -
974 $
sdot( j-ki-2, t( ki+2, j ), 1,
975 $ work( ki+2+n ), 1 )
977 work( j+n2 ) = work( j+n2 ) -
978 $
sdot( j-ki-2, t( ki+2, j ), 1,
979 $ work( ki+2+n2 ), 1 )
981 work( j+1+n ) = work( j+1+n ) -
982 $
sdot( j-ki-2, t( ki+2, j+1 ), 1,
983 $ work( ki+2+n ), 1 )
985 work( j+1+n2 ) = work( j+1+n2 ) -
986 $
sdot( j-ki-2, t( ki+2, j+1 ), 1,
987 $ work( ki+2+n2 ), 1 )
993 CALL slaln2( .true., 2, 2, smin, one, t( j, j ),
994 $ ldt, one, one, work( j+n ), n, wr,
995 $ -wi, x, 2, scale, xnorm, ierr )
999 IF( scale.NE.one )
THEN
1000 CALL sscal( n-ki+1, scale, work( ki+n ), 1 )
1001 CALL sscal( n-ki+1, scale, work( ki+n2 ), 1 )
1003 work( j+n ) = x( 1, 1 )
1004 work( j+n2 ) = x( 1, 2 )
1005 work( j+1+n ) = x( 2, 1 )
1006 work( j+1+n2 ) = x( 2, 2 )
1007 vmax = max( abs( x( 1, 1 ) ), abs( x( 1, 2 ) ),
1008 $ abs( x( 2, 1 ) ), abs( x( 2, 2 ) ), vmax )
1009 vcrit = bignum / vmax
1016 IF( .NOT.over )
THEN
1017 CALL scopy( n-ki+1, work( ki+n ), 1, vl( ki, is ), 1 )
1018 CALL scopy( n-ki+1, work( ki+n2 ), 1, vl( ki, is+1 ),
1023 emax = max( emax, abs( vl( k, is ) )+
1024 $ abs( vl( k, is+1 ) ) )
1027 CALL sscal( n-ki+1, remax, vl( ki, is ), 1 )
1028 CALL sscal( n-ki+1, remax, vl( ki, is+1 ), 1 )
1030 DO 230 k = 1, ki - 1
1032 vl( k, is+1 ) = zero
1035 IF( ki.LT.n-1 )
THEN
1036 CALL sgemv(
'N', n, n-ki-1, one, vl( 1, ki+2 ),
1037 $ ldvl, work( ki+2+n ), 1, work( ki+n ),
1039 CALL sgemv(
'N', n, n-ki-1, one, vl( 1, ki+2 ),
1040 $ ldvl, work( ki+2+n2 ), 1,
1041 $ work( ki+1+n2 ), vl( 1, ki+1 ), 1 )
1043 CALL sscal( n, work( ki+n ), vl( 1, ki ), 1 )
1044 CALL sscal( n, work( ki+1+n2 ), vl( 1, ki+1 ), 1 )
1049 emax = max( emax, abs( vl( k, ki ) )+
1050 $ abs( vl( k, ki+1 ) ) )
1053 CALL sscal( n, remax, vl( 1, ki ), 1 )
1054 CALL sscal( n, remax, vl( 1, ki+1 ), 1 )
real function sdot(N, SX, INCX, SY, INCY)
SDOT
integer function isamax(N, SX, INCX)
ISAMAX
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
subroutine slaln2(LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, LDB, WR, WI, X, LDX, SCALE, XNORM, INFO)
SLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
subroutine sscal(N, SA, SX, INCX)
SSCAL
real function slamch(CMACH)
SLAMCH
logical function lsame(CA, CB)
LSAME
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY