239 SUBROUTINE strevc3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
240 $ vr, ldvr, mm, m, work, lwork, info )
249 CHARACTER HOWMNY, SIDE
250 INTEGER INFO, LDT, LDVL, LDVR, LWORK, M, MM, N
254 REAL T( ldt, * ), VL( ldvl, * ), VR( ldvr, * ),
262 parameter ( zero = 0.0e+0, one = 1.0e+0 )
264 parameter ( nbmin = 8, nbmax = 128 )
267 LOGICAL ALLV, BOTHV, LEFTV, LQUERY, OVER, PAIR,
269 INTEGER I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI,
270 $ iv, maxwrk, nb, ki2
271 REAL BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
272 $ smin, smlnum, ulp, unfl, vcrit, vmax, wi, wr,
277 INTEGER ISAMAX, ILAENV
279 EXTERNAL lsame, isamax, ilaenv, sdot, slamch
286 INTRINSIC abs, max, sqrt
290 INTEGER ISCOMPLEX( nbmax )
296 bothv = lsame( side,
'B' )
297 rightv = lsame( side,
'R' ) .OR. bothv
298 leftv = lsame( side,
'L' ) .OR. bothv
300 allv = lsame( howmny,
'A' )
301 over = lsame( howmny,
'B' )
302 somev = lsame( howmny,
'S' )
305 nb = ilaenv( 1,
'STREVC', side // howmny, n, -1, -1, -1 )
308 lquery = ( lwork.EQ.-1 )
309 IF( .NOT.rightv .AND. .NOT.leftv )
THEN
311 ELSE IF( .NOT.allv .AND. .NOT.over .AND. .NOT.somev )
THEN
313 ELSE IF( n.LT.0 )
THEN
315 ELSE IF( ldt.LT.max( 1, n ) )
THEN
317 ELSE IF( ldvl.LT.1 .OR. ( leftv .AND. ldvl.LT.n ) )
THEN
319 ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) )
THEN
321 ELSE IF( lwork.LT.max( 1, 3*n ) .AND. .NOT.lquery )
THEN
335 SELECT( j ) = .false.
338 IF( t( j+1, j ).EQ.zero )
THEN
343 IF(
SELECT( j ) .OR.
SELECT( j+1 ) )
THEN
363 CALL xerbla(
'STREVC3', -info )
365 ELSE IF( lquery )
THEN
377 IF( over .AND. lwork .GE. n + 2*n*nbmin )
THEN
378 nb = (lwork - n) / (2*n)
379 nb = min( nb, nbmax )
380 CALL slaset(
'F', n, 1+2*nb, zero, zero, work, n )
387 unfl = slamch(
'Safe minimum' )
390 ulp = slamch(
'Precision' )
391 smlnum = unfl*( n / ulp )
392 bignum = ( one-ulp ) / smlnum
401 work( j ) = work( j ) + abs( t( i, j ) )
434 ELSE IF( ki.EQ.1 )
THEN
437 ELSE IF( t( ki, ki-1 ).EQ.zero )
THEN
447 IF( .NOT.
SELECT( ki ) )
450 IF( .NOT.
SELECT( ki-1 ) )
460 $ wi = sqrt( abs( t( ki, ki-1 ) ) )*
461 $ sqrt( abs( t( ki-1, ki ) ) )
462 smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum )
469 work( ki + iv*n ) = one
474 work( k + iv*n ) = -t( k, ki )
481 DO 60 j = ki - 1, 1, -1
488 IF( t( j, j-1 ).NE.zero )
THEN
498 CALL slaln2( .false., 1, 1, smin, one, t( j, j ),
499 $ ldt, one, one, work( j+iv*n ), n, wr,
500 $ zero, x, 2, scale, xnorm, ierr )
505 IF( xnorm.GT.one )
THEN
506 IF( work( j ).GT.bignum / xnorm )
THEN
507 x( 1, 1 ) = x( 1, 1 ) / xnorm
508 scale = scale / xnorm
515 $
CALL sscal( ki, scale, work( 1+iv*n ), 1 )
516 work( j+iv*n ) = x( 1, 1 )
520 CALL saxpy( j-1, -x( 1, 1 ), t( 1, j ), 1,
521 $ work( 1+iv*n ), 1 )
527 CALL slaln2( .false., 2, 1, smin, one,
528 $ t( j-1, j-1 ), ldt, one, one,
529 $ work( j-1+iv*n ), n, wr, zero, x, 2,
530 $ scale, xnorm, ierr )
535 IF( xnorm.GT.one )
THEN
536 beta = max( work( j-1 ), work( j ) )
537 IF( beta.GT.bignum / xnorm )
THEN
538 x( 1, 1 ) = x( 1, 1 ) / xnorm
539 x( 2, 1 ) = x( 2, 1 ) / xnorm
540 scale = scale / xnorm
547 $
CALL sscal( ki, scale, work( 1+iv*n ), 1 )
548 work( j-1+iv*n ) = x( 1, 1 )
549 work( j +iv*n ) = x( 2, 1 )
553 CALL saxpy( j-2, -x( 1, 1 ), t( 1, j-1 ), 1,
554 $ work( 1+iv*n ), 1 )
555 CALL saxpy( j-2, -x( 2, 1 ), t( 1, j ), 1,
556 $ work( 1+iv*n ), 1 )
565 CALL scopy( ki, work( 1 + iv*n ), 1, vr( 1, is ), 1 )
567 ii = isamax( ki, vr( 1, is ), 1 )
568 remax = one / abs( vr( ii, is ) )
569 CALL sscal( ki, remax, vr( 1, is ), 1 )
575 ELSE IF( nb.EQ.1 )
THEN
579 $
CALL sgemv(
'N', n, ki-1, one, vr, ldvr,
580 $ work( 1 + iv*n ), 1, work( ki + iv*n ),
583 ii = isamax( n, vr( 1, ki ), 1 )
584 remax = one / abs( vr( ii, ki ) )
585 CALL sscal( n, remax, vr( 1, ki ), 1 )
592 work( k + iv*n ) = zero
606 IF( abs( t( ki-1, ki ) ).GE.abs( t( ki, ki-1 ) ) )
THEN
607 work( ki-1 + (iv-1)*n ) = one
608 work( ki + (iv )*n ) = wi / t( ki-1, ki )
610 work( ki-1 + (iv-1)*n ) = -wi / t( ki, ki-1 )
611 work( ki + (iv )*n ) = one
613 work( ki + (iv-1)*n ) = zero
614 work( ki-1 + (iv )*n ) = zero
619 work( k+(iv-1)*n ) = -work( ki-1+(iv-1)*n )*t(k,ki-1)
620 work( k+(iv )*n ) = -work( ki +(iv )*n )*t(k,ki )
627 DO 90 j = ki - 2, 1, -1
634 IF( t( j, j-1 ).NE.zero )
THEN
644 CALL slaln2( .false., 1, 2, smin, one, t( j, j ),
645 $ ldt, one, one, work( j+(iv-1)*n ), n,
646 $ wr, wi, x, 2, scale, xnorm, ierr )
651 IF( xnorm.GT.one )
THEN
652 IF( work( j ).GT.bignum / xnorm )
THEN
653 x( 1, 1 ) = x( 1, 1 ) / xnorm
654 x( 1, 2 ) = x( 1, 2 ) / xnorm
655 scale = scale / xnorm
661 IF( scale.NE.one )
THEN
662 CALL sscal( ki, scale, work( 1+(iv-1)*n ), 1 )
663 CALL sscal( ki, scale, work( 1+(iv )*n ), 1 )
665 work( j+(iv-1)*n ) = x( 1, 1 )
666 work( j+(iv )*n ) = x( 1, 2 )
670 CALL saxpy( j-1, -x( 1, 1 ), t( 1, j ), 1,
671 $ work( 1+(iv-1)*n ), 1 )
672 CALL saxpy( j-1, -x( 1, 2 ), t( 1, j ), 1,
673 $ work( 1+(iv )*n ), 1 )
679 CALL slaln2( .false., 2, 2, smin, one,
680 $ t( j-1, j-1 ), ldt, one, one,
681 $ work( j-1+(iv-1)*n ), n, wr, wi, x, 2,
682 $ scale, xnorm, ierr )
687 IF( xnorm.GT.one )
THEN
688 beta = max( work( j-1 ), work( j ) )
689 IF( beta.GT.bignum / xnorm )
THEN
691 x( 1, 1 ) = x( 1, 1 )*rec
692 x( 1, 2 ) = x( 1, 2 )*rec
693 x( 2, 1 ) = x( 2, 1 )*rec
694 x( 2, 2 ) = x( 2, 2 )*rec
701 IF( scale.NE.one )
THEN
702 CALL sscal( ki, scale, work( 1+(iv-1)*n ), 1 )
703 CALL sscal( ki, scale, work( 1+(iv )*n ), 1 )
705 work( j-1+(iv-1)*n ) = x( 1, 1 )
706 work( j +(iv-1)*n ) = x( 2, 1 )
707 work( j-1+(iv )*n ) = x( 1, 2 )
708 work( j +(iv )*n ) = x( 2, 2 )
712 CALL saxpy( j-2, -x( 1, 1 ), t( 1, j-1 ), 1,
713 $ work( 1+(iv-1)*n ), 1 )
714 CALL saxpy( j-2, -x( 2, 1 ), t( 1, j ), 1,
715 $ work( 1+(iv-1)*n ), 1 )
716 CALL saxpy( j-2, -x( 1, 2 ), t( 1, j-1 ), 1,
717 $ work( 1+(iv )*n ), 1 )
718 CALL saxpy( j-2, -x( 2, 2 ), t( 1, j ), 1,
719 $ work( 1+(iv )*n ), 1 )
728 CALL scopy( ki, work( 1+(iv-1)*n ), 1, vr(1,is-1), 1 )
729 CALL scopy( ki, work( 1+(iv )*n ), 1, vr(1,is ), 1 )
733 emax = max( emax, abs( vr( k, is-1 ) )+
734 $ abs( vr( k, is ) ) )
737 CALL sscal( ki, remax, vr( 1, is-1 ), 1 )
738 CALL sscal( ki, remax, vr( 1, is ), 1 )
745 ELSE IF( nb.EQ.1 )
THEN
749 CALL sgemv(
'N', n, ki-2, one, vr, ldvr,
750 $ work( 1 + (iv-1)*n ), 1,
751 $ work( ki-1 + (iv-1)*n ), vr(1,ki-1), 1)
752 CALL sgemv(
'N', n, ki-2, one, vr, ldvr,
753 $ work( 1 + (iv)*n ), 1,
754 $ work( ki + (iv)*n ), vr( 1, ki ), 1 )
756 CALL sscal( n, work(ki-1+(iv-1)*n), vr(1,ki-1), 1)
757 CALL sscal( n, work(ki +(iv )*n), vr(1,ki ), 1)
762 emax = max( emax, abs( vr( k, ki-1 ) )+
763 $ abs( vr( k, ki ) ) )
766 CALL sscal( n, remax, vr( 1, ki-1 ), 1 )
767 CALL sscal( n, remax, vr( 1, ki ), 1 )
774 work( k + (iv-1)*n ) = zero
775 work( k + (iv )*n ) = zero
777 iscomplex( iv-1 ) = -ip
797 IF( (iv.LE.2) .OR. (ki2.EQ.1) )
THEN
798 CALL sgemm(
'N',
'N', n, nb-iv+1, ki2+nb-iv, one,
800 $ work( 1 + (iv)*n ), n,
802 $ work( 1 + (nb+iv)*n ), n )
805 IF( iscomplex(k).EQ.0 )
THEN
807 ii = isamax( n, work( 1 + (nb+k)*n ), 1 )
808 remax = one / abs( work( ii + (nb+k)*n ) )
809 ELSE IF( iscomplex(k).EQ.1 )
THEN
814 $ abs( work( ii + (nb+k )*n ) )+
815 $ abs( work( ii + (nb+k+1)*n ) ) )
822 CALL sscal( n, remax, work( 1 + (nb+k)*n ), 1 )
824 CALL slacpy(
'F', n, nb-iv+1,
825 $ work( 1 + (nb+iv)*n ), n,
826 $ vr( 1, ki2 ), ldvr )
858 ELSE IF( ki.EQ.n )
THEN
861 ELSE IF( t( ki+1, ki ).EQ.zero )
THEN
870 IF( .NOT.
SELECT( ki ) )
879 $ wi = sqrt( abs( t( ki, ki+1 ) ) )*
880 $ sqrt( abs( t( ki+1, ki ) ) )
881 smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum )
888 work( ki + iv*n ) = one
893 work( k + iv*n ) = -t( ki, k )
910 IF( t( j+1, j ).NE.zero )
THEN
923 IF( work( j ).GT.vcrit )
THEN
925 CALL sscal( n-ki+1, rec, work( ki+iv*n ), 1 )
930 work( j+iv*n ) = work( j+iv*n ) -
931 $ sdot( j-ki-1, t( ki+1, j ), 1,
932 $ work( ki+1+iv*n ), 1 )
936 CALL slaln2( .false., 1, 1, smin, one, t( j, j ),
937 $ ldt, one, one, work( j+iv*n ), n, wr,
938 $ zero, x, 2, scale, xnorm, ierr )
943 $
CALL sscal( n-ki+1, scale, work( ki+iv*n ), 1 )
944 work( j+iv*n ) = x( 1, 1 )
945 vmax = max( abs( work( j+iv*n ) ), vmax )
946 vcrit = bignum / vmax
955 beta = max( work( j ), work( j+1 ) )
956 IF( beta.GT.vcrit )
THEN
958 CALL sscal( n-ki+1, rec, work( ki+iv*n ), 1 )
963 work( j+iv*n ) = work( j+iv*n ) -
964 $ sdot( j-ki-1, t( ki+1, j ), 1,
965 $ work( ki+1+iv*n ), 1 )
967 work( j+1+iv*n ) = work( j+1+iv*n ) -
968 $ sdot( j-ki-1, t( ki+1, j+1 ), 1,
969 $ work( ki+1+iv*n ), 1 )
975 CALL slaln2( .true., 2, 1, smin, one, t( j, j ),
976 $ ldt, one, one, work( j+iv*n ), n, wr,
977 $ zero, x, 2, scale, xnorm, ierr )
982 $
CALL sscal( n-ki+1, scale, work( ki+iv*n ), 1 )
983 work( j +iv*n ) = x( 1, 1 )
984 work( j+1+iv*n ) = x( 2, 1 )
986 vmax = max( abs( work( j +iv*n ) ),
987 $ abs( work( j+1+iv*n ) ), vmax )
988 vcrit = bignum / vmax
998 CALL scopy( n-ki+1, work( ki + iv*n ), 1,
1001 ii = isamax( n-ki+1, vl( ki, is ), 1 ) + ki - 1
1002 remax = one / abs( vl( ii, is ) )
1003 CALL sscal( n-ki+1, remax, vl( ki, is ), 1 )
1005 DO 180 k = 1, ki - 1
1009 ELSE IF( nb.EQ.1 )
THEN
1013 $
CALL sgemv(
'N', n, n-ki, one,
1014 $ vl( 1, ki+1 ), ldvl,
1015 $ work( ki+1 + iv*n ), 1,
1016 $ work( ki + iv*n ), vl( 1, ki ), 1 )
1018 ii = isamax( n, vl( 1, ki ), 1 )
1019 remax = one / abs( vl( ii, ki ) )
1020 CALL sscal( n, remax, vl( 1, ki ), 1 )
1028 work( k + iv*n ) = zero
1030 iscomplex( iv ) = ip
1042 IF( abs( t( ki, ki+1 ) ).GE.abs( t( ki+1, ki ) ) )
THEN
1043 work( ki + (iv )*n ) = wi / t( ki, ki+1 )
1044 work( ki+1 + (iv+1)*n ) = one
1046 work( ki + (iv )*n ) = one
1047 work( ki+1 + (iv+1)*n ) = -wi / t( ki+1, ki )
1049 work( ki+1 + (iv )*n ) = zero
1050 work( ki + (iv+1)*n ) = zero
1054 DO 190 k = ki + 2, n
1055 work( k+(iv )*n ) = -work( ki +(iv )*n )*t(ki, k)
1056 work( k+(iv+1)*n ) = -work( ki+1+(iv+1)*n )*t(ki+1,k)
1066 DO 200 j = ki + 2, n
1073 IF( t( j+1, j ).NE.zero )
THEN
1086 IF( work( j ).GT.vcrit )
THEN
1088 CALL sscal( n-ki+1, rec, work(ki+(iv )*n), 1 )
1089 CALL sscal( n-ki+1, rec, work(ki+(iv+1)*n), 1 )
1094 work( j+(iv )*n ) = work( j+(iv)*n ) -
1095 $ sdot( j-ki-2, t( ki+2, j ), 1,
1096 $ work( ki+2+(iv)*n ), 1 )
1097 work( j+(iv+1)*n ) = work( j+(iv+1)*n ) -
1098 $ sdot( j-ki-2, t( ki+2, j ), 1,
1099 $ work( ki+2+(iv+1)*n ), 1 )
1103 CALL slaln2( .false., 1, 2, smin, one, t( j, j ),
1104 $ ldt, one, one, work( j+iv*n ), n, wr,
1105 $ -wi, x, 2, scale, xnorm, ierr )
1109 IF( scale.NE.one )
THEN
1110 CALL sscal( n-ki+1, scale, work(ki+(iv )*n), 1)
1111 CALL sscal( n-ki+1, scale, work(ki+(iv+1)*n), 1)
1113 work( j+(iv )*n ) = x( 1, 1 )
1114 work( j+(iv+1)*n ) = x( 1, 2 )
1115 vmax = max( abs( work( j+(iv )*n ) ),
1116 $ abs( work( j+(iv+1)*n ) ), vmax )
1117 vcrit = bignum / vmax
1126 beta = max( work( j ), work( j+1 ) )
1127 IF( beta.GT.vcrit )
THEN
1129 CALL sscal( n-ki+1, rec, work(ki+(iv )*n), 1 )
1130 CALL sscal( n-ki+1, rec, work(ki+(iv+1)*n), 1 )
1135 work( j +(iv )*n ) = work( j+(iv)*n ) -
1136 $ sdot( j-ki-2, t( ki+2, j ), 1,
1137 $ work( ki+2+(iv)*n ), 1 )
1139 work( j +(iv+1)*n ) = work( j+(iv+1)*n ) -
1140 $ sdot( j-ki-2, t( ki+2, j ), 1,
1141 $ work( ki+2+(iv+1)*n ), 1 )
1143 work( j+1+(iv )*n ) = work( j+1+(iv)*n ) -
1144 $ sdot( j-ki-2, t( ki+2, j+1 ), 1,
1145 $ work( ki+2+(iv)*n ), 1 )
1147 work( j+1+(iv+1)*n ) = work( j+1+(iv+1)*n ) -
1148 $ sdot( j-ki-2, t( ki+2, j+1 ), 1,
1149 $ work( ki+2+(iv+1)*n ), 1 )
1155 CALL slaln2( .true., 2, 2, smin, one, t( j, j ),
1156 $ ldt, one, one, work( j+iv*n ), n, wr,
1157 $ -wi, x, 2, scale, xnorm, ierr )
1161 IF( scale.NE.one )
THEN
1162 CALL sscal( n-ki+1, scale, work(ki+(iv )*n), 1)
1163 CALL sscal( n-ki+1, scale, work(ki+(iv+1)*n), 1)
1165 work( j +(iv )*n ) = x( 1, 1 )
1166 work( j +(iv+1)*n ) = x( 1, 2 )
1167 work( j+1+(iv )*n ) = x( 2, 1 )
1168 work( j+1+(iv+1)*n ) = x( 2, 2 )
1169 vmax = max( abs( x( 1, 1 ) ), abs( x( 1, 2 ) ),
1170 $ abs( x( 2, 1 ) ), abs( x( 2, 2 ) ),
1172 vcrit = bignum / vmax
1179 IF( .NOT.over )
THEN
1182 CALL scopy( n-ki+1, work( ki + (iv )*n ), 1,
1184 CALL scopy( n-ki+1, work( ki + (iv+1)*n ), 1,
1185 $ vl( ki, is+1 ), 1 )
1189 emax = max( emax, abs( vl( k, is ) )+
1190 $ abs( vl( k, is+1 ) ) )
1193 CALL sscal( n-ki+1, remax, vl( ki, is ), 1 )
1194 CALL sscal( n-ki+1, remax, vl( ki, is+1 ), 1 )
1196 DO 230 k = 1, ki - 1
1198 vl( k, is+1 ) = zero
1201 ELSE IF( nb.EQ.1 )
THEN
1204 IF( ki.LT.n-1 )
THEN
1205 CALL sgemv(
'N', n, n-ki-1, one,
1206 $ vl( 1, ki+2 ), ldvl,
1207 $ work( ki+2 + (iv)*n ), 1,
1208 $ work( ki + (iv)*n ),
1210 CALL sgemv(
'N', n, n-ki-1, one,
1211 $ vl( 1, ki+2 ), ldvl,
1212 $ work( ki+2 + (iv+1)*n ), 1,
1213 $ work( ki+1 + (iv+1)*n ),
1214 $ vl( 1, ki+1 ), 1 )
1216 CALL sscal( n, work(ki+ (iv )*n), vl(1, ki ), 1)
1217 CALL sscal( n, work(ki+1+(iv+1)*n), vl(1, ki+1), 1)
1222 emax = max( emax, abs( vl( k, ki ) )+
1223 $ abs( vl( k, ki+1 ) ) )
1226 CALL sscal( n, remax, vl( 1, ki ), 1 )
1227 CALL sscal( n, remax, vl( 1, ki+1 ), 1 )
1235 work( k + (iv )*n ) = zero
1236 work( k + (iv+1)*n ) = zero
1238 iscomplex( iv ) = ip
1239 iscomplex( iv+1 ) = -ip
1258 IF( (iv.GE.nb-1) .OR. (ki2.EQ.n) )
THEN
1259 CALL sgemm(
'N',
'N', n, iv, n-ki2+iv, one,
1260 $ vl( 1, ki2-iv+1 ), ldvl,
1261 $ work( ki2-iv+1 + (1)*n ), n,
1263 $ work( 1 + (nb+1)*n ), n )
1266 IF( iscomplex(k).EQ.0)
THEN
1268 ii = isamax( n, work( 1 + (nb+k)*n ), 1 )
1269 remax = one / abs( work( ii + (nb+k)*n ) )
1270 ELSE IF( iscomplex(k).EQ.1)
THEN
1275 $ abs( work( ii + (nb+k )*n ) )+
1276 $ abs( work( ii + (nb+k+1)*n ) ) )
1283 CALL sscal( n, remax, work( 1 + (nb+k)*n ), 1 )
1286 $ work( 1 + (nb+1)*n ), n,
1287 $ vl( 1, ki2-iv+1 ), ldvl )
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
subroutine strevc3(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, LWORK, INFO)
STREVC3
subroutine slabad(SMALL, LARGE)
SLABAD
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 sscal(N, SA, SX, INCX)
SSCAL
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY