246 SUBROUTINE ctrevc3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
247 $ LDVR, MM, M, WORK, LWORK, RWORK, LRWORK, INFO)
256 CHARACTER HOWMNY, SIDE
257 INTEGER INFO, LDT, LDVL, LDVR, LWORK, LRWORK, M, MM, N
262 COMPLEX T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
270 parameter( zero = 0.0e+0, one = 1.0e+0 )
272 parameter( czero = ( 0.0e+0, 0.0e+0 ),
273 $ cone = ( 1.0e+0, 0.0e+0 ) )
275 parameter( nbmin = 8, nbmax = 128 )
278 LOGICAL ALLV, BOTHV, LEFTV, LQUERY, OVER, RIGHTV, SOMEV
279 INTEGER I, II, IS, J, K, KI, IV, MAXWRK, NB
280 REAL OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL
285 INTEGER ILAENV, ICAMAX
287 EXTERNAL lsame, ilaenv, icamax, slamch, scasum
294 INTRINSIC abs, real, cmplx, conjg, aimag, max
300 cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
306 bothv = lsame( side,
'B' )
307 rightv = lsame( side,
'R' ) .OR. bothv
308 leftv = lsame( side,
'L' ) .OR. bothv
310 allv = lsame( howmny,
'A' )
311 over = lsame( howmny,
'B' )
312 somev = lsame( howmny,
'S' )
328 nb = ilaenv( 1,
'CTREVC', side // howmny, n, -1, -1, -1 )
332 lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 )
333 IF( .NOT.rightv .AND. .NOT.leftv )
THEN
335 ELSE IF( .NOT.allv .AND. .NOT.over .AND. .NOT.somev )
THEN
337 ELSE IF( n.LT.0 )
THEN
339 ELSE IF( ldt.LT.max( 1, n ) )
THEN
341 ELSE IF( ldvl.LT.1 .OR. ( leftv .AND. ldvl.LT.n ) )
THEN
343 ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) )
THEN
345 ELSE IF( mm.LT.m )
THEN
347 ELSE IF( lwork.LT.max( 1, 2*n ) .AND. .NOT.lquery )
THEN
349 ELSE IF ( lrwork.LT.max( 1, n ) .AND. .NOT.lquery )
THEN
353 CALL xerbla(
'CTREVC3', -info )
355 ELSE IF( lquery )
THEN
367 IF( over .AND. lwork .GE. n + 2*n*nbmin )
THEN
368 nb = (lwork - n) / (2*n)
369 nb = min( nb, nbmax )
370 CALL claset(
'F', n, 1+2*nb, czero, czero, work, n )
377 unfl = slamch(
'Safe minimum' )
380 ulp = slamch(
'Precision' )
381 smlnum = unfl*( n / ulp )
386 work( i ) = t( i, i )
394 rwork( j ) = scasum( j-1, t( 1, j ), 1 )
410 IF( .NOT.
SELECT( ki ) )
413 smin = max( ulp*( cabs1( t( ki, ki ) ) ), smlnum )
418 work( ki + iv*n ) = cone
423 work( k + iv*n ) = -t( k, ki )
430 t( k, k ) = t( k, k ) - t( ki, ki )
431 IF( cabs1( t( k, k ) ).LT.smin )
436 CALL clatrs(
'Upper',
'No transpose',
'Non-unit',
'Y',
437 $ ki-1, t, ldt, work( 1 + iv*n ), scale,
439 work( ki + iv*n ) = scale
447 CALL ccopy( ki, work( 1 + iv*n ), 1, vr( 1, is ), 1 )
449 ii = icamax( ki, vr( 1, is ), 1 )
450 remax = one / cabs1( vr( ii, is ) )
451 CALL csscal( ki, remax, vr( 1, is ), 1 )
457 ELSE IF( nb.EQ.1 )
THEN
461 $
CALL cgemv(
'N', n, ki-1, cone, vr, ldvr,
462 $ work( 1 + iv*n ), 1, cmplx( scale ),
465 ii = icamax( n, vr( 1, ki ), 1 )
466 remax = one / cabs1( vr( ii, ki ) )
467 CALL csscal( n, remax, vr( 1, ki ), 1 )
474 work( k + iv*n ) = czero
480 IF( (iv.EQ.1) .OR. (ki.EQ.1) )
THEN
481 CALL cgemm(
'N',
'N', n, nb-iv+1, ki+nb-iv, cone,
483 $ work( 1 + (iv)*n ), n,
485 $ work( 1 + (nb+iv)*n ), n )
488 ii = icamax( n, work( 1 + (nb+k)*n ), 1 )
489 remax = one / cabs1( work( ii + (nb+k)*n ) )
490 CALL csscal( n, remax, work( 1 + (nb+k)*n ), 1 )
492 CALL clacpy(
'F', n, nb-iv+1,
493 $ work( 1 + (nb+iv)*n ), n,
494 $ vr( 1, ki ), ldvr )
504 t( k, k ) = work( k )
525 IF( .NOT.
SELECT( ki ) )
528 smin = max( ulp*( cabs1( t( ki, ki ) ) ), smlnum )
533 work( ki + iv*n ) = cone
538 work( k + iv*n ) = -conjg( t( ki, k ) )
545 t( k, k ) = t( k, k ) - t( ki, ki )
546 IF( cabs1( t( k, k ) ).LT.smin )
551 CALL clatrs(
'Upper',
'Conjugate transpose',
'Non-unit',
552 $
'Y', n-ki, t( ki+1, ki+1 ), ldt,
553 $ work( ki+1 + iv*n ), scale, rwork, info )
554 work( ki + iv*n ) = scale
562 CALL ccopy( n-ki+1, work( ki + iv*n ), 1, vl(ki,is), 1 )
564 ii = icamax( n-ki+1, vl( ki, is ), 1 ) + ki - 1
565 remax = one / cabs1( vl( ii, is ) )
566 CALL csscal( n-ki+1, remax, vl( ki, is ), 1 )
572 ELSE IF( nb.EQ.1 )
THEN
576 $
CALL cgemv(
'N', n, n-ki, cone, vl( 1, ki+1 ), ldvl,
577 $ work( ki+1 + iv*n ), 1, cmplx( scale ),
580 ii = icamax( n, vl( 1, ki ), 1 )
581 remax = one / cabs1( vl( ii, ki ) )
582 CALL csscal( n, remax, vl( 1, ki ), 1 )
590 work( k + iv*n ) = czero
596 IF( (iv.EQ.nb) .OR. (ki.EQ.n) )
THEN
597 CALL cgemm(
'N',
'N', n, iv, n-ki+iv, cone,
598 $ vl( 1, ki-iv+1 ), ldvl,
599 $ work( ki-iv+1 + (1)*n ), n,
601 $ work( 1 + (nb+1)*n ), n )
604 ii = icamax( n, work( 1 + (nb+k)*n ), 1 )
605 remax = one / cabs1( work( ii + (nb+k)*n ) )
606 CALL csscal( n, remax, work( 1 + (nb+k)*n ), 1 )
609 $ work( 1 + (nb+1)*n ), n,
610 $ vl( 1, ki-iv+1 ), ldvl )
620 t( k, k ) = work( k )