215 SUBROUTINE ztgevc( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
216 $ LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
223 CHARACTER HOWMNY, SIDE
224 INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
228 DOUBLE PRECISION RWORK( * )
229 COMPLEX*16 P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
230 $ vr( ldvr, * ), work( * )
237 DOUBLE PRECISION ZERO, ONE
238 parameter( zero = 0.0d+0, one = 1.0d+0 )
239 COMPLEX*16 CZERO, CONE
240 parameter( czero = ( 0.0d+0, 0.0d+0 ),
241 $ cone = ( 1.0d+0, 0.0d+0 ) )
244 LOGICAL COMPL, COMPR, ILALL, ILBACK, ILBBAD, ILCOMP,
246 INTEGER I, IBEG, IEIG, IEND, IHWMNY, IM, ISIDE, ISRC,
248 DOUBLE PRECISION ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG,
249 $ bignum, bnorm, bscale, dmin, safmin, sbeta,
250 $ scale, small, temp, ulp, xmax
251 COMPLEX*16 BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X
255 DOUBLE PRECISION DLAMCH
257 EXTERNAL lsame, dlamch, zladiv
263 INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, min
266 DOUBLE PRECISION ABS1
269 abs1( x ) = abs( dble( x ) ) + abs( dimag( x ) )
275 IF( lsame( howmny,
'A' ) )
THEN
279 ELSE IF( lsame( howmny,
'S' ) )
THEN
283 ELSE IF( lsame( howmny,
'B' ) )
THEN
291 IF( lsame( side,
'R' ) )
THEN
295 ELSE IF( lsame( side,
'L' ) )
THEN
299 ELSE IF( lsame( side,
'B' ) )
THEN
308 IF( iside.LT.0 )
THEN
310 ELSE IF( ihwmny.LT.0 )
THEN
312 ELSE IF( n.LT.0 )
THEN
314 ELSE IF( lds.LT.max( 1, n ) )
THEN
316 ELSE IF( ldp.LT.max( 1, n ) )
THEN
320 CALL xerbla(
'ZTGEVC', -info )
326 IF( .NOT.ilall )
THEN
340 IF( dimag( p( j, j ) ).NE.zero )
346 ELSE IF( compl .AND. ldvl.LT.n .OR. ldvl.LT.1 )
THEN
348 ELSE IF( compr .AND. ldvr.LT.n .OR. ldvr.LT.1 )
THEN
350 ELSE IF( mm.LT.im )
THEN
354 CALL xerbla(
'ZTGEVC', -info )
366 safmin = dlamch(
'Safe minimum' )
368 ulp = dlamch(
'Epsilon' )*dlamch(
'Base' )
369 small = safmin*n / ulp
371 bignum = one / ( safmin*n )
377 anorm = abs1( s( 1, 1 ) )
378 bnorm = abs1( p( 1, 1 ) )
385 rwork( j ) = rwork( j ) + abs1( s( i, j ) )
386 rwork( n+j ) = rwork( n+j ) + abs1( p( i, j ) )
388 anorm = max( anorm, rwork( j )+abs1( s( j, j ) ) )
389 bnorm = max( bnorm, rwork( n+j )+abs1( p( j, j ) ) )
392 ascale = one / max( anorm, safmin )
393 bscale = one / max( bnorm, safmin )
406 ilcomp =
SELECT( je )
411 IF( abs1( s( je, je ) ).LE.safmin .AND.
412 $ abs( dble( p( je, je ) ) ).LE.safmin )
THEN
417 vl( jr, ieig ) = czero
419 vl( ieig, ieig ) = cone
428 temp = one / max( abs1( s( je, je ) )*ascale,
429 $ abs( dble( p( je, je ) ) )*bscale, safmin )
430 salpha = ( temp*s( je, je ) )*ascale
431 sbeta = ( temp*dble( p( je, je ) ) )*bscale
432 acoeff = sbeta*ascale
433 bcoeff = salpha*bscale
437 lsa = abs( sbeta ).GE.safmin .AND. abs( acoeff ).LT.small
438 lsb = abs1( salpha ).GE.safmin .AND. abs1( bcoeff ).LT.
443 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
445 $ scale = max( scale, ( small / abs1( salpha ) )*
446 $ min( bnorm, big ) )
447 IF( lsa .OR. lsb )
THEN
448 scale = min( scale, one /
449 $ ( safmin*max( one, abs( acoeff ),
450 $ abs1( bcoeff ) ) ) )
452 acoeff = ascale*( scale*sbeta )
454 acoeff = scale*acoeff
457 bcoeff = bscale*( scale*salpha )
459 bcoeff = scale*bcoeff
463 acoefa = abs( acoeff )
464 bcoefa = abs1( bcoeff )
470 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
487 IF( acoefa*rwork( j )+bcoefa*rwork( n+j ).GT.bignum*
490 work( jr ) = temp*work( jr )
498 suma = suma + dconjg( s( jr, j ) )*work( jr )
499 sumb = sumb + dconjg( p( jr, j ) )*work( jr )
501 sum = acoeff*suma - dconjg( bcoeff )*sumb
507 d = dconjg( acoeff*s( j, j )-bcoeff*p( j, j ) )
508 IF( abs1( d ).LE.dmin )
511 IF( abs1( d ).LT.one )
THEN
512 IF( abs1( sum ).GE.bignum*abs1( d ) )
THEN
513 temp = one / abs1( sum )
515 work( jr ) = temp*work( jr )
521 work( j ) = zladiv( -sum, d )
522 xmax = max( xmax, abs1( work( j ) ) )
528 CALL zgemv(
'N', n, n+1-je, cone, vl( 1, je ),
530 $ work( je ), 1, czero, work( n+1 ), 1 )
542 xmax = max( xmax, abs1( work( ( isrc-1 )*n+jr ) ) )
545 IF( xmax.GT.safmin )
THEN
548 vl( jr, ieig ) = temp*work( ( isrc-1 )*n+jr )
554 DO 130 jr = 1, ibeg - 1
555 vl( jr, ieig ) = czero
573 ilcomp =
SELECT( je )
578 IF( abs1( s( je, je ) ).LE.safmin .AND.
579 $ abs( dble( p( je, je ) ) ).LE.safmin )
THEN
584 vr( jr, ieig ) = czero
586 vr( ieig, ieig ) = cone
595 temp = one / max( abs1( s( je, je ) )*ascale,
596 $ abs( dble( p( je, je ) ) )*bscale, safmin )
597 salpha = ( temp*s( je, je ) )*ascale
598 sbeta = ( temp*dble( p( je, je ) ) )*bscale
599 acoeff = sbeta*ascale
600 bcoeff = salpha*bscale
604 lsa = abs( sbeta ).GE.safmin .AND. abs( acoeff ).LT.small
605 lsb = abs1( salpha ).GE.safmin .AND. abs1( bcoeff ).LT.
610 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
612 $ scale = max( scale, ( small / abs1( salpha ) )*
613 $ min( bnorm, big ) )
614 IF( lsa .OR. lsb )
THEN
615 scale = min( scale, one /
616 $ ( safmin*max( one, abs( acoeff ),
617 $ abs1( bcoeff ) ) ) )
619 acoeff = ascale*( scale*sbeta )
621 acoeff = scale*acoeff
624 bcoeff = bscale*( scale*salpha )
626 bcoeff = scale*bcoeff
630 acoefa = abs( acoeff )
631 bcoefa = abs1( bcoeff )
637 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
644 DO 170 jr = 1, je - 1
645 work( jr ) = acoeff*s( jr, je ) - bcoeff*p( jr, je )
649 DO 210 j = je - 1, 1, -1
654 d = acoeff*s( j, j ) - bcoeff*p( j, j )
655 IF( abs1( d ).LE.dmin )
658 IF( abs1( d ).LT.one )
THEN
659 IF( abs1( work( j ) ).GE.bignum*abs1( d ) )
THEN
660 temp = one / abs1( work( j ) )
662 work( jr ) = temp*work( jr )
667 work( j ) = zladiv( -work( j ), d )
673 IF( abs1( work( j ) ).GT.one )
THEN
674 temp = one / abs1( work( j ) )
675 IF( acoefa*rwork( j )+bcoefa*rwork( n+j ).GE.
678 work( jr ) = temp*work( jr )
683 ca = acoeff*work( j )
684 cb = bcoeff*work( j )
686 work( jr ) = work( jr ) + ca*s( jr, j ) -
695 CALL zgemv(
'N', n, je, cone, vr, ldvr, work, 1,
696 $ czero, work( n+1 ), 1 )
708 xmax = max( xmax, abs1( work( ( isrc-1 )*n+jr ) ) )
711 IF( xmax.GT.safmin )
THEN
714 vr( jr, ieig ) = temp*work( ( isrc-1 )*n+jr )
720 DO 240 jr = iend + 1, n
721 vr( jr, ieig ) = czero