239 SUBROUTINE dtrevc3( 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 DOUBLE PRECISION T( ldt, * ), VL( ldvl, * ), VR( ldvr, * ),
261 DOUBLE PRECISION ZERO, ONE
262 parameter ( zero = 0.0d+0, one = 1.0d+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 DOUBLE PRECISION BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
272 $ smin, smlnum, ulp, unfl, vcrit, vmax, wi, wr,
277 INTEGER IDAMAX, ILAENV
278 DOUBLE PRECISION DDOT, DLAMCH
279 EXTERNAL lsame, idamax, ilaenv, ddot, dlamch
286 INTRINSIC abs, max, sqrt
289 DOUBLE PRECISION X( 2, 2 )
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,
'DTREVC', 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(
'DTREVC3', -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 dlaset(
'F', n, 1+2*nb, zero, zero, work, n )
387 unfl = dlamch(
'Safe minimum' )
390 ulp = dlamch(
'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 dlaln2( .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 dscal( ki, scale, work( 1+iv*n ), 1 )
516 work( j+iv*n ) = x( 1, 1 )
520 CALL daxpy( j-1, -x( 1, 1 ), t( 1, j ), 1,
521 $ work( 1+iv*n ), 1 )
527 CALL dlaln2( .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 dscal( 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 daxpy( j-2, -x( 1, 1 ), t( 1, j-1 ), 1,
554 $ work( 1+iv*n ), 1 )
555 CALL daxpy( j-2, -x( 2, 1 ), t( 1, j ), 1,
556 $ work( 1+iv*n ), 1 )
565 CALL dcopy( ki, work( 1 + iv*n ), 1, vr( 1, is ), 1 )
567 ii = idamax( ki, vr( 1, is ), 1 )
568 remax = one / abs( vr( ii, is ) )
569 CALL dscal( ki, remax, vr( 1, is ), 1 )
575 ELSE IF( nb.EQ.1 )
THEN
579 $
CALL dgemv(
'N', n, ki-1, one, vr, ldvr,
580 $ work( 1 + iv*n ), 1, work( ki + iv*n ),
583 ii = idamax( n, vr( 1, ki ), 1 )
584 remax = one / abs( vr( ii, ki ) )
585 CALL dscal( 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 dlaln2( .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 dscal( ki, scale, work( 1+(iv-1)*n ), 1 )
663 CALL dscal( 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 daxpy( j-1, -x( 1, 1 ), t( 1, j ), 1,
671 $ work( 1+(iv-1)*n ), 1 )
672 CALL daxpy( j-1, -x( 1, 2 ), t( 1, j ), 1,
673 $ work( 1+(iv )*n ), 1 )
679 CALL dlaln2( .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 dscal( ki, scale, work( 1+(iv-1)*n ), 1 )
703 CALL dscal( 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 daxpy( j-2, -x( 1, 1 ), t( 1, j-1 ), 1,
713 $ work( 1+(iv-1)*n ), 1 )
714 CALL daxpy( j-2, -x( 2, 1 ), t( 1, j ), 1,
715 $ work( 1+(iv-1)*n ), 1 )
716 CALL daxpy( j-2, -x( 1, 2 ), t( 1, j-1 ), 1,
717 $ work( 1+(iv )*n ), 1 )
718 CALL daxpy( j-2, -x( 2, 2 ), t( 1, j ), 1,
719 $ work( 1+(iv )*n ), 1 )
728 CALL dcopy( ki, work( 1+(iv-1)*n ), 1, vr(1,is-1), 1 )
729 CALL dcopy( 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 dscal( ki, remax, vr( 1, is-1 ), 1 )
738 CALL dscal( ki, remax, vr( 1, is ), 1 )
745 ELSE IF( nb.EQ.1 )
THEN
749 CALL dgemv(
'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 dgemv(
'N', n, ki-2, one, vr, ldvr,
753 $ work( 1 + (iv)*n ), 1,
754 $ work( ki + (iv)*n ), vr( 1, ki ), 1 )
756 CALL dscal( n, work(ki-1+(iv-1)*n), vr(1,ki-1), 1)
757 CALL dscal( 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 dscal( n, remax, vr( 1, ki-1 ), 1 )
767 CALL dscal( 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 dgemm(
'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 = idamax( 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 dscal( n, remax, work( 1 + (nb+k)*n ), 1 )
824 CALL dlacpy(
'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 dscal( n-ki+1, rec, work( ki+iv*n ), 1 )
930 work( j+iv*n ) = work( j+iv*n ) -
931 $ ddot( j-ki-1, t( ki+1, j ), 1,
932 $ work( ki+1+iv*n ), 1 )
936 CALL dlaln2( .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 dscal( 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 dscal( n-ki+1, rec, work( ki+iv*n ), 1 )
963 work( j+iv*n ) = work( j+iv*n ) -
964 $ ddot( 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 $ ddot( j-ki-1, t( ki+1, j+1 ), 1,
969 $ work( ki+1+iv*n ), 1 )
975 CALL dlaln2( .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 dscal( 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 dcopy( n-ki+1, work( ki + iv*n ), 1,
1001 ii = idamax( n-ki+1, vl( ki, is ), 1 ) + ki - 1
1002 remax = one / abs( vl( ii, is ) )
1003 CALL dscal( n-ki+1, remax, vl( ki, is ), 1 )
1005 DO 180 k = 1, ki - 1
1009 ELSE IF( nb.EQ.1 )
THEN
1013 $
CALL dgemv(
'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 = idamax( n, vl( 1, ki ), 1 )
1019 remax = one / abs( vl( ii, ki ) )
1020 CALL dscal( 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 dscal( n-ki+1, rec, work(ki+(iv )*n), 1 )
1089 CALL dscal( n-ki+1, rec, work(ki+(iv+1)*n), 1 )
1094 work( j+(iv )*n ) = work( j+(iv)*n ) -
1095 $ ddot( 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 $ ddot( j-ki-2, t( ki+2, j ), 1,
1099 $ work( ki+2+(iv+1)*n ), 1 )
1103 CALL dlaln2( .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 dscal( n-ki+1, scale, work(ki+(iv )*n), 1)
1111 CALL dscal( 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 dscal( n-ki+1, rec, work(ki+(iv )*n), 1 )
1130 CALL dscal( n-ki+1, rec, work(ki+(iv+1)*n), 1 )
1135 work( j +(iv )*n ) = work( j+(iv)*n ) -
1136 $ ddot( 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 $ ddot( 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 $ ddot( 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 $ ddot( j-ki-2, t( ki+2, j+1 ), 1,
1149 $ work( ki+2+(iv+1)*n ), 1 )
1155 CALL dlaln2( .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 dscal( n-ki+1, scale, work(ki+(iv )*n), 1)
1163 CALL dscal( 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 dcopy( n-ki+1, work( ki + (iv )*n ), 1,
1184 CALL dcopy( 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 dscal( n-ki+1, remax, vl( ki, is ), 1 )
1194 CALL dscal( 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 dgemv(
'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 dgemv(
'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 dscal( n, work(ki+ (iv )*n), vl(1, ki ), 1)
1217 CALL dscal( 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 dscal( n, remax, vl( 1, ki ), 1 )
1227 CALL dscal( 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 dgemm(
'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 = idamax( 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 dscal( n, remax, work( 1 + (nb+k)*n ), 1 )
1286 $ work( 1 + (nb+1)*n ), n,
1287 $ vl( 1, ki2-iv+1 ), ldvl )
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dtrevc3(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, LWORK, INFO)
DTREVC3
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
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.
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine dlabad(SMALL, LARGE)
DLABAD