235 SUBROUTINE strevc3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
236 $ VR, LDVR, MM, M, WORK, LWORK, INFO )
244 CHARACTER HOWMNY, SIDE
245 INTEGER INFO, LDT, LDVL, LDVR, LWORK, M, MM, N
249 REAL T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
257 parameter( zero = 0.0e+0, one = 1.0e+0 )
259 parameter( nbmin = 8, nbmax = 128 )
262 LOGICAL ALLV, BOTHV, LEFTV, LQUERY, OVER, PAIR,
264 INTEGER I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI,
265 $ iv, maxwrk, nb, ki2
266 REAL BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
267 $ smin, smlnum, ulp, unfl, vcrit, vmax, wi, wr,
272 INTEGER ISAMAX, ILAENV
274 EXTERNAL lsame, isamax, ilaenv, sdot, slamch
281 INTRINSIC abs, max, sqrt
285 INTEGER ISCOMPLEX( NBMAX )
291 bothv = lsame( side,
'B' )
292 rightv = lsame( side,
'R' ) .OR. bothv
293 leftv = lsame( side,
'L' ) .OR. bothv
295 allv = lsame( howmny,
'A' )
296 over = lsame( howmny,
'B' )
297 somev = lsame( howmny,
'S' )
300 nb = ilaenv( 1,
'STREVC', side // howmny, n, -1, -1, -1 )
303 lquery = ( lwork.EQ.-1 )
304 IF( .NOT.rightv .AND. .NOT.leftv )
THEN
306 ELSE IF( .NOT.allv .AND. .NOT.over .AND. .NOT.somev )
THEN
308 ELSE IF( n.LT.0 )
THEN
310 ELSE IF( ldt.LT.max( 1, n ) )
THEN
312 ELSE IF( ldvl.LT.1 .OR. ( leftv .AND. ldvl.LT.n ) )
THEN
314 ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) )
THEN
316 ELSE IF( lwork.LT.max( 1, 3*n ) .AND. .NOT.lquery )
THEN
330 SELECT( j ) = .false.
333 IF( t( j+1, j ).EQ.zero )
THEN
338 IF(
SELECT( j ) .OR.
SELECT( j+1 ) )
THEN
358 CALL xerbla(
'STREVC3', -info )
360 ELSE IF( lquery )
THEN
372 IF( over .AND. lwork .GE. n + 2*n*nbmin )
THEN
373 nb = (lwork - n) / (2*n)
374 nb = min( nb, nbmax )
375 CALL slaset(
'F', n, 1+2*nb, zero, zero, work, n )
382 unfl = slamch(
'Safe minimum' )
385 ulp = slamch(
'Precision' )
386 smlnum = unfl*( n / ulp )
387 bignum = ( one-ulp ) / smlnum
396 work( j ) = work( j ) + abs( t( i, j ) )
429 ELSE IF( ki.EQ.1 )
THEN
432 ELSE IF( t( ki, ki-1 ).EQ.zero )
THEN
442 IF( .NOT.
SELECT( ki ) )
445 IF( .NOT.
SELECT( ki-1 ) )
455 $ wi = sqrt( abs( t( ki, ki-1 ) ) )*
456 $ sqrt( abs( t( ki-1, ki ) ) )
457 smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum )
464 work( ki + iv*n ) = one
469 work( k + iv*n ) = -t( k, ki )
476 DO 60 j = ki - 1, 1, -1
483 IF( t( j, j-1 ).NE.zero )
THEN
493 CALL slaln2( .false., 1, 1, smin, one, t( j, j ),
494 $ ldt, one, one, work( j+iv*n ), n, wr,
495 $ zero, x, 2, scale, xnorm, ierr )
500 IF( xnorm.GT.one )
THEN
501 IF( work( j ).GT.bignum / xnorm )
THEN
502 x( 1, 1 ) = x( 1, 1 ) / xnorm
503 scale = scale / xnorm
510 $
CALL sscal( ki, scale, work( 1+iv*n ), 1 )
511 work( j+iv*n ) = x( 1, 1 )
515 CALL saxpy( j-1, -x( 1, 1 ), t( 1, j ), 1,
516 $ work( 1+iv*n ), 1 )
522 CALL slaln2( .false., 2, 1, smin, one,
523 $ t( j-1, j-1 ), ldt, one, one,
524 $ work( j-1+iv*n ), n, wr, zero, x, 2,
525 $ scale, xnorm, ierr )
530 IF( xnorm.GT.one )
THEN
531 beta = max( work( j-1 ), work( j ) )
532 IF( beta.GT.bignum / xnorm )
THEN
533 x( 1, 1 ) = x( 1, 1 ) / xnorm
534 x( 2, 1 ) = x( 2, 1 ) / xnorm
535 scale = scale / xnorm
542 $
CALL sscal( ki, scale, work( 1+iv*n ), 1 )
543 work( j-1+iv*n ) = x( 1, 1 )
544 work( j +iv*n ) = x( 2, 1 )
548 CALL saxpy( j-2, -x( 1, 1 ), t( 1, j-1 ), 1,
549 $ work( 1+iv*n ), 1 )
550 CALL saxpy( j-2, -x( 2, 1 ), t( 1, j ), 1,
551 $ work( 1+iv*n ), 1 )
560 CALL scopy( ki, work( 1 + iv*n ), 1, vr( 1, is ), 1 )
562 ii = isamax( ki, vr( 1, is ), 1 )
563 remax = one / abs( vr( ii, is ) )
564 CALL sscal( ki, remax, vr( 1, is ), 1 )
570 ELSE IF( nb.EQ.1 )
THEN
574 $
CALL sgemv(
'N', n, ki-1, one, vr, ldvr,
575 $ work( 1 + iv*n ), 1, work( ki + iv*n ),
578 ii = isamax( n, vr( 1, ki ), 1 )
579 remax = one / abs( vr( ii, ki ) )
580 CALL sscal( n, remax, vr( 1, ki ), 1 )
587 work( k + iv*n ) = zero
601 IF( abs( t( ki-1, ki ) ).GE.abs( t( ki, ki-1 ) ) )
THEN
602 work( ki-1 + (iv-1)*n ) = one
603 work( ki + (iv )*n ) = wi / t( ki-1, ki )
605 work( ki-1 + (iv-1)*n ) = -wi / t( ki, ki-1 )
606 work( ki + (iv )*n ) = one
608 work( ki + (iv-1)*n ) = zero
609 work( ki-1 + (iv )*n ) = zero
614 work( k+(iv-1)*n ) = -work( ki-1+(iv-1)*n )*t(k,ki-1)
615 work( k+(iv )*n ) = -work( ki +(iv )*n )*t(k,ki )
622 DO 90 j = ki - 2, 1, -1
629 IF( t( j, j-1 ).NE.zero )
THEN
639 CALL slaln2( .false., 1, 2, smin, one, t( j, j ),
640 $ ldt, one, one, work( j+(iv-1)*n ), n,
641 $ wr, wi, x, 2, scale, xnorm, ierr )
646 IF( xnorm.GT.one )
THEN
647 IF( work( j ).GT.bignum / xnorm )
THEN
648 x( 1, 1 ) = x( 1, 1 ) / xnorm
649 x( 1, 2 ) = x( 1, 2 ) / xnorm
650 scale = scale / xnorm
656 IF( scale.NE.one )
THEN
657 CALL sscal( ki, scale, work( 1+(iv-1)*n ), 1 )
658 CALL sscal( ki, scale, work( 1+(iv )*n ), 1 )
660 work( j+(iv-1)*n ) = x( 1, 1 )
661 work( j+(iv )*n ) = x( 1, 2 )
665 CALL saxpy( j-1, -x( 1, 1 ), t( 1, j ), 1,
666 $ work( 1+(iv-1)*n ), 1 )
667 CALL saxpy( j-1, -x( 1, 2 ), t( 1, j ), 1,
668 $ work( 1+(iv )*n ), 1 )
674 CALL slaln2( .false., 2, 2, smin, one,
675 $ t( j-1, j-1 ), ldt, one, one,
676 $ work( j-1+(iv-1)*n ), n, wr, wi, x, 2,
677 $ scale, xnorm, ierr )
682 IF( xnorm.GT.one )
THEN
683 beta = max( work( j-1 ), work( j ) )
684 IF( beta.GT.bignum / xnorm )
THEN
686 x( 1, 1 ) = x( 1, 1 )*rec
687 x( 1, 2 ) = x( 1, 2 )*rec
688 x( 2, 1 ) = x( 2, 1 )*rec
689 x( 2, 2 ) = x( 2, 2 )*rec
696 IF( scale.NE.one )
THEN
697 CALL sscal( ki, scale, work( 1+(iv-1)*n ), 1 )
698 CALL sscal( ki, scale, work( 1+(iv )*n ), 1 )
700 work( j-1+(iv-1)*n ) = x( 1, 1 )
701 work( j +(iv-1)*n ) = x( 2, 1 )
702 work( j-1+(iv )*n ) = x( 1, 2 )
703 work( j +(iv )*n ) = x( 2, 2 )
707 CALL saxpy( j-2, -x( 1, 1 ), t( 1, j-1 ), 1,
708 $ work( 1+(iv-1)*n ), 1 )
709 CALL saxpy( j-2, -x( 2, 1 ), t( 1, j ), 1,
710 $ work( 1+(iv-1)*n ), 1 )
711 CALL saxpy( j-2, -x( 1, 2 ), t( 1, j-1 ), 1,
712 $ work( 1+(iv )*n ), 1 )
713 CALL saxpy( j-2, -x( 2, 2 ), t( 1, j ), 1,
714 $ work( 1+(iv )*n ), 1 )
723 CALL scopy( ki, work( 1+(iv-1)*n ), 1, vr(1,is-1), 1 )
724 CALL scopy( ki, work( 1+(iv )*n ), 1, vr(1,is ), 1 )
728 emax = max( emax, abs( vr( k, is-1 ) )+
729 $ abs( vr( k, is ) ) )
732 CALL sscal( ki, remax, vr( 1, is-1 ), 1 )
733 CALL sscal( ki, remax, vr( 1, is ), 1 )
740 ELSE IF( nb.EQ.1 )
THEN
744 CALL sgemv(
'N', n, ki-2, one, vr, ldvr,
745 $ work( 1 + (iv-1)*n ), 1,
746 $ work( ki-1 + (iv-1)*n ), vr(1,ki-1), 1)
747 CALL sgemv(
'N', n, ki-2, one, vr, ldvr,
748 $ work( 1 + (iv)*n ), 1,
749 $ work( ki + (iv)*n ), vr( 1, ki ), 1 )
751 CALL sscal( n, work(ki-1+(iv-1)*n), vr(1,ki-1), 1)
752 CALL sscal( n, work(ki +(iv )*n), vr(1,ki ), 1)
757 emax = max( emax, abs( vr( k, ki-1 ) )+
758 $ abs( vr( k, ki ) ) )
761 CALL sscal( n, remax, vr( 1, ki-1 ), 1 )
762 CALL sscal( n, remax, vr( 1, ki ), 1 )
769 work( k + (iv-1)*n ) = zero
770 work( k + (iv )*n ) = zero
772 iscomplex( iv-1 ) = -ip
792 IF( (iv.LE.2) .OR. (ki2.EQ.1) )
THEN
793 CALL sgemm(
'N',
'N', n, nb-iv+1, ki2+nb-iv, one,
795 $ work( 1 + (iv)*n ), n,
797 $ work( 1 + (nb+iv)*n ), n )
800 IF( iscomplex(k).EQ.0 )
THEN
802 ii = isamax( n, work( 1 + (nb+k)*n ), 1 )
803 remax = one / abs( work( ii + (nb+k)*n ) )
804 ELSE IF( iscomplex(k).EQ.1 )
THEN
809 $ abs( work( ii + (nb+k )*n ) )+
810 $ abs( work( ii + (nb+k+1)*n ) ) )
817 CALL sscal( n, remax, work( 1 + (nb+k)*n ), 1 )
819 CALL slacpy(
'F', n, nb-iv+1,
820 $ work( 1 + (nb+iv)*n ), n,
821 $ vr( 1, ki2 ), ldvr )
853 ELSE IF( ki.EQ.n )
THEN
856 ELSE IF( t( ki+1, ki ).EQ.zero )
THEN
865 IF( .NOT.
SELECT( ki ) )
874 $ wi = sqrt( abs( t( ki, ki+1 ) ) )*
875 $ sqrt( abs( t( ki+1, ki ) ) )
876 smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum )
883 work( ki + iv*n ) = one
888 work( k + iv*n ) = -t( ki, k )
905 IF( t( j+1, j ).NE.zero )
THEN
918 IF( work( j ).GT.vcrit )
THEN
920 CALL sscal( n-ki+1, rec, work( ki+iv*n ), 1 )
925 work( j+iv*n ) = work( j+iv*n ) -
926 $ sdot( j-ki-1, t( ki+1, j ), 1,
927 $ work( ki+1+iv*n ), 1 )
931 CALL slaln2( .false., 1, 1, smin, one, t( j, j ),
932 $ ldt, one, one, work( j+iv*n ), n, wr,
933 $ zero, x, 2, scale, xnorm, ierr )
938 $
CALL sscal( n-ki+1, scale, work( ki+iv*n ), 1 )
939 work( j+iv*n ) = x( 1, 1 )
940 vmax = max( abs( work( j+iv*n ) ), vmax )
941 vcrit = bignum / vmax
950 beta = max( work( j ), work( j+1 ) )
951 IF( beta.GT.vcrit )
THEN
953 CALL sscal( n-ki+1, rec, work( ki+iv*n ), 1 )
958 work( j+iv*n ) = work( j+iv*n ) -
959 $ sdot( j-ki-1, t( ki+1, j ), 1,
960 $ work( ki+1+iv*n ), 1 )
962 work( j+1+iv*n ) = work( j+1+iv*n ) -
963 $ sdot( j-ki-1, t( ki+1, j+1 ), 1,
964 $ work( ki+1+iv*n ), 1 )
970 CALL slaln2( .true., 2, 1, smin, one, t( j, j ),
971 $ ldt, one, one, work( j+iv*n ), n, wr,
972 $ zero, x, 2, scale, xnorm, ierr )
977 $
CALL sscal( n-ki+1, scale, work( ki+iv*n ), 1 )
978 work( j +iv*n ) = x( 1, 1 )
979 work( j+1+iv*n ) = x( 2, 1 )
981 vmax = max( abs( work( j +iv*n ) ),
982 $ abs( work( j+1+iv*n ) ), vmax )
983 vcrit = bignum / vmax
993 CALL scopy( n-ki+1, work( ki + iv*n ), 1,
996 ii = isamax( n-ki+1, vl( ki, is ), 1 ) + ki - 1
997 remax = one / abs( vl( ii, is ) )
998 CALL sscal( n-ki+1, remax, vl( ki, is ), 1 )
1000 DO 180 k = 1, ki - 1
1004 ELSE IF( nb.EQ.1 )
THEN
1008 $
CALL sgemv(
'N', n, n-ki, one,
1009 $ vl( 1, ki+1 ), ldvl,
1010 $ work( ki+1 + iv*n ), 1,
1011 $ work( ki + iv*n ), vl( 1, ki ), 1 )
1013 ii = isamax( n, vl( 1, ki ), 1 )
1014 remax = one / abs( vl( ii, ki ) )
1015 CALL sscal( n, remax, vl( 1, ki ), 1 )
1023 work( k + iv*n ) = zero
1025 iscomplex( iv ) = ip
1037 IF( abs( t( ki, ki+1 ) ).GE.abs( t( ki+1, ki ) ) )
THEN
1038 work( ki + (iv )*n ) = wi / t( ki, ki+1 )
1039 work( ki+1 + (iv+1)*n ) = one
1041 work( ki + (iv )*n ) = one
1042 work( ki+1 + (iv+1)*n ) = -wi / t( ki+1, ki )
1044 work( ki+1 + (iv )*n ) = zero
1045 work( ki + (iv+1)*n ) = zero
1049 DO 190 k = ki + 2, n
1050 work( k+(iv )*n ) = -work( ki +(iv )*n )*t(ki, k)
1051 work( k+(iv+1)*n ) = -work( ki+1+(iv+1)*n )*t(ki+1,k)
1061 DO 200 j = ki + 2, n
1068 IF( t( j+1, j ).NE.zero )
THEN
1081 IF( work( j ).GT.vcrit )
THEN
1083 CALL sscal( n-ki+1, rec, work(ki+(iv )*n), 1 )
1084 CALL sscal( n-ki+1, rec, work(ki+(iv+1)*n), 1 )
1089 work( j+(iv )*n ) = work( j+(iv)*n ) -
1090 $ sdot( j-ki-2, t( ki+2, j ), 1,
1091 $ work( ki+2+(iv)*n ), 1 )
1092 work( j+(iv+1)*n ) = work( j+(iv+1)*n ) -
1093 $ sdot( j-ki-2, t( ki+2, j ), 1,
1094 $ work( ki+2+(iv+1)*n ), 1 )
1098 CALL slaln2( .false., 1, 2, smin, one, t( j, j ),
1099 $ ldt, one, one, work( j+iv*n ), n, wr,
1100 $ -wi, x, 2, scale, xnorm, ierr )
1104 IF( scale.NE.one )
THEN
1105 CALL sscal( n-ki+1, scale, work(ki+(iv )*n), 1)
1106 CALL sscal( n-ki+1, scale, work(ki+(iv+1)*n), 1)
1108 work( j+(iv )*n ) = x( 1, 1 )
1109 work( j+(iv+1)*n ) = x( 1, 2 )
1110 vmax = max( abs( work( j+(iv )*n ) ),
1111 $ abs( work( j+(iv+1)*n ) ), vmax )
1112 vcrit = bignum / vmax
1121 beta = max( work( j ), work( j+1 ) )
1122 IF( beta.GT.vcrit )
THEN
1124 CALL sscal( n-ki+1, rec, work(ki+(iv )*n), 1 )
1125 CALL sscal( n-ki+1, rec, work(ki+(iv+1)*n), 1 )
1130 work( j +(iv )*n ) = work( j+(iv)*n ) -
1131 $ sdot( j-ki-2, t( ki+2, j ), 1,
1132 $ work( ki+2+(iv)*n ), 1 )
1134 work( j +(iv+1)*n ) = work( j+(iv+1)*n ) -
1135 $ sdot( j-ki-2, t( ki+2, j ), 1,
1136 $ work( ki+2+(iv+1)*n ), 1 )
1138 work( j+1+(iv )*n ) = work( j+1+(iv)*n ) -
1139 $ sdot( j-ki-2, t( ki+2, j+1 ), 1,
1140 $ work( ki+2+(iv)*n ), 1 )
1142 work( j+1+(iv+1)*n ) = work( j+1+(iv+1)*n ) -
1143 $ sdot( j-ki-2, t( ki+2, j+1 ), 1,
1144 $ work( ki+2+(iv+1)*n ), 1 )
1150 CALL slaln2( .true., 2, 2, smin, one, t( j, j ),
1151 $ ldt, one, one, work( j+iv*n ), n, wr,
1152 $ -wi, x, 2, scale, xnorm, ierr )
1156 IF( scale.NE.one )
THEN
1157 CALL sscal( n-ki+1, scale, work(ki+(iv )*n), 1)
1158 CALL sscal( n-ki+1, scale, work(ki+(iv+1)*n), 1)
1160 work( j +(iv )*n ) = x( 1, 1 )
1161 work( j +(iv+1)*n ) = x( 1, 2 )
1162 work( j+1+(iv )*n ) = x( 2, 1 )
1163 work( j+1+(iv+1)*n ) = x( 2, 2 )
1164 vmax = max( abs( x( 1, 1 ) ), abs( x( 1, 2 ) ),
1165 $ abs( x( 2, 1 ) ), abs( x( 2, 2 ) ),
1167 vcrit = bignum / vmax
1174 IF( .NOT.over )
THEN
1177 CALL scopy( n-ki+1, work( ki + (iv )*n ), 1,
1179 CALL scopy( n-ki+1, work( ki + (iv+1)*n ), 1,
1180 $ vl( ki, is+1 ), 1 )
1184 emax = max( emax, abs( vl( k, is ) )+
1185 $ abs( vl( k, is+1 ) ) )
1188 CALL sscal( n-ki+1, remax, vl( ki, is ), 1 )
1189 CALL sscal( n-ki+1, remax, vl( ki, is+1 ), 1 )
1191 DO 230 k = 1, ki - 1
1193 vl( k, is+1 ) = zero
1196 ELSE IF( nb.EQ.1 )
THEN
1199 IF( ki.LT.n-1 )
THEN
1200 CALL sgemv(
'N', n, n-ki-1, one,
1201 $ vl( 1, ki+2 ), ldvl,
1202 $ work( ki+2 + (iv)*n ), 1,
1203 $ work( ki + (iv)*n ),
1205 CALL sgemv(
'N', n, n-ki-1, one,
1206 $ vl( 1, ki+2 ), ldvl,
1207 $ work( ki+2 + (iv+1)*n ), 1,
1208 $ work( ki+1 + (iv+1)*n ),
1209 $ vl( 1, ki+1 ), 1 )
1211 CALL sscal( n, work(ki+ (iv )*n), vl(1, ki ), 1)
1212 CALL sscal( n, work(ki+1+(iv+1)*n), vl(1, ki+1), 1)
1217 emax = max( emax, abs( vl( k, ki ) )+
1218 $ abs( vl( k, ki+1 ) ) )
1221 CALL sscal( n, remax, vl( 1, ki ), 1 )
1222 CALL sscal( n, remax, vl( 1, ki+1 ), 1 )
1230 work( k + (iv )*n ) = zero
1231 work( k + (iv+1)*n ) = zero
1233 iscomplex( iv ) = ip
1234 iscomplex( iv+1 ) = -ip
1253 IF( (iv.GE.nb-1) .OR. (ki2.EQ.n) )
THEN
1254 CALL sgemm(
'N',
'N', n, iv, n-ki2+iv, one,
1255 $ vl( 1, ki2-iv+1 ), ldvl,
1256 $ work( ki2-iv+1 + (1)*n ), n,
1258 $ work( 1 + (nb+1)*n ), n )
1261 IF( iscomplex(k).EQ.0)
THEN
1263 ii = isamax( n, work( 1 + (nb+k)*n ), 1 )
1264 remax = one / abs( work( ii + (nb+k)*n ) )
1265 ELSE IF( iscomplex(k).EQ.1)
THEN
1270 $ abs( work( ii + (nb+k )*n ) )+
1271 $ abs( work( ii + (nb+k+1)*n ) ) )
1278 CALL sscal( n, remax, work( 1 + (nb+k)*n ), 1 )
1281 $ work( 1 + (nb+1)*n ), n,
1282 $ vl( 1, ki2-iv+1 ), ldvl )
subroutine slabad(SMALL, LARGE)
SLABAD
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 slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine xerbla(SRNAME, INFO)
XERBLA
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 strevc3(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, LWORK, INFO)
STREVC3
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine sscal(N, SA, SX, INCX)
SSCAL
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM