219 SUBROUTINE ctgevc( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
220 $ ldvl, vr, ldvr, mm, m, work, rwork, info )
228 CHARACTER howmny, side
229 INTEGER info, ldp, lds, ldvl, ldvr, m, mm, n
234 COMPLEX p( ldp, * ), s( lds, * ), vl( ldvl, * ),
235 $ vr( ldvr, * ), work( * )
243 parameter( zero = 0.0e+0, one = 1.0e+0 )
245 parameter( czero = ( 0.0e+0, 0.0e+0 ),
246 $ cone = ( 1.0e+0, 0.0e+0 ) )
249 LOGICAL compl, compr, ilall, ilback, ilbbad, ilcomp,
251 INTEGER i, ibeg, ieig, iend, ihwmny, im, iside, isrc,
253 REAL acoefa, acoeff, anorm, ascale, bcoefa, big,
254 $ bignum, bnorm, bscale, dmin, safmin, sbeta,
255 $ scale, small, temp, ulp, xmax
256 COMPLEX bcoeff, ca, cb, d, salpha, sum, suma, sumb, x
268 INTRINSIC abs, aimag, cmplx, conjg, max, min, real
274 abs1( x ) = abs(
REAL( X ) ) + abs( aimag( x ) )
280 IF(
lsame( howmny,
'A' ) )
THEN
284 ELSE IF(
lsame( howmny,
'S' ) )
THEN
288 ELSE IF(
lsame( howmny,
'B' ) )
THEN
296 IF(
lsame( side,
'R' ) )
THEN
300 ELSE IF(
lsame( side,
'L' ) )
THEN
304 ELSE IF(
lsame( side,
'B' ) )
THEN
313 IF( iside.LT.0 )
THEN
315 ELSE IF( ihwmny.LT.0 )
THEN
317 ELSE IF( n.LT.0 )
THEN
319 ELSE IF( lds.LT.max( 1, n ) )
THEN
321 ELSE IF( ldp.LT.max( 1, n ) )
THEN
325 CALL
xerbla(
'CTGEVC', -info )
331 IF( .NOT.ilall )
THEN
345 IF( aimag( p( j, j ) ).NE.zero )
351 ELSE IF( compl .AND. ldvl.LT.n .OR. ldvl.LT.1 )
THEN
353 ELSE IF( compr .AND. ldvr.LT.n .OR. ldvr.LT.1 )
THEN
355 ELSE IF( mm.LT.im )
THEN
359 CALL
xerbla(
'CTGEVC', -info )
371 safmin =
slamch(
'Safe minimum' )
373 CALL
slabad( safmin, big )
375 small = safmin*n / ulp
377 bignum = one / ( safmin*n )
383 anorm = abs1( s( 1, 1 ) )
384 bnorm = abs1( p( 1, 1 ) )
391 rwork( j ) = rwork( j ) + abs1( s( i, j ) )
392 rwork( n+j ) = rwork( n+j ) + abs1( p( i, j ) )
394 anorm = max( anorm, rwork( j )+abs1( s( j, j ) ) )
395 bnorm = max( bnorm, rwork( n+j )+abs1( p( j, j ) ) )
398 ascale = one / max( anorm, safmin )
399 bscale = one / max( bnorm, safmin )
412 ilcomp =
SELECT( je )
417 IF( abs1( s( je, je ) ).LE.safmin .AND.
418 $ abs(
REAL( P( JE, JE ) ) ).LE.safmin ) then
423 vl( jr, ieig ) = czero
425 vl( ieig, ieig ) = cone
434 temp = one / max( abs1( s( je, je ) )*ascale,
435 $ abs(
REAL( P( JE, JE ) ) )*bscale, safmin )
436 salpha = ( temp*s( je, je ) )*ascale
437 sbeta = ( temp*
REAL( P( JE, JE ) ) )*bscale
438 acoeff = sbeta*ascale
439 bcoeff = salpha*bscale
443 lsa = abs( sbeta ).GE.safmin .AND. abs( acoeff ).LT.small
444 lsb = abs1( salpha ).GE.safmin .AND. abs1( bcoeff ).LT.
449 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
451 $ scale = max( scale, ( small / abs1( salpha ) )*
452 $ min( bnorm, big ) )
453 IF( lsa .OR. lsb )
THEN
454 scale = min( scale, one /
455 $ ( safmin*max( one, abs( acoeff ),
456 $ abs1( bcoeff ) ) ) )
458 acoeff = ascale*( scale*sbeta )
460 acoeff = scale*acoeff
463 bcoeff = bscale*( scale*salpha )
465 bcoeff = scale*bcoeff
469 acoefa = abs( acoeff )
470 bcoefa = abs1( bcoeff )
476 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
493 IF( acoefa*rwork( j )+bcoefa*rwork( n+j ).GT.bignum*
496 work( jr ) = temp*work( jr )
504 suma = suma + conjg( s( jr, j ) )*work( jr )
505 sumb = sumb + conjg( p( jr, j ) )*work( jr )
507 sum = acoeff*suma - conjg( bcoeff )*sumb
513 d = conjg( acoeff*s( j, j )-bcoeff*p( j, j ) )
514 IF( abs1( d ).LE.dmin )
517 IF( abs1( d ).LT.one )
THEN
518 IF( abs1( sum ).GE.bignum*abs1( d ) )
THEN
519 temp = one / abs1( sum )
521 work( jr ) = temp*work( jr )
527 work( j ) =
cladiv( -sum, d )
528 xmax = max( xmax, abs1( work( j ) ) )
534 CALL
cgemv(
'N', n, n+1-je, cone, vl( 1, je ), ldvl,
535 $ work( je ), 1, czero, work( n+1 ), 1 )
547 xmax = max( xmax, abs1( work( ( isrc-1 )*n+jr ) ) )
550 IF( xmax.GT.safmin )
THEN
553 vl( jr, ieig ) = temp*work( ( isrc-1 )*n+jr )
559 DO 130 jr = 1, ibeg - 1
560 vl( jr, ieig ) = czero
578 ilcomp =
SELECT( je )
583 IF( abs1( s( je, je ) ).LE.safmin .AND.
584 $ abs(
REAL( P( JE, JE ) ) ).LE.safmin ) then
589 vr( jr, ieig ) = czero
591 vr( ieig, ieig ) = cone
600 temp = one / max( abs1( s( je, je ) )*ascale,
601 $ abs(
REAL( P( JE, JE ) ) )*bscale, safmin )
602 salpha = ( temp*s( je, je ) )*ascale
603 sbeta = ( temp*
REAL( P( JE, JE ) ) )*bscale
604 acoeff = sbeta*ascale
605 bcoeff = salpha*bscale
609 lsa = abs( sbeta ).GE.safmin .AND. abs( acoeff ).LT.small
610 lsb = abs1( salpha ).GE.safmin .AND. abs1( bcoeff ).LT.
615 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
617 $ scale = max( scale, ( small / abs1( salpha ) )*
618 $ min( bnorm, big ) )
619 IF( lsa .OR. lsb )
THEN
620 scale = min( scale, one /
621 $ ( safmin*max( one, abs( acoeff ),
622 $ abs1( bcoeff ) ) ) )
624 acoeff = ascale*( scale*sbeta )
626 acoeff = scale*acoeff
629 bcoeff = bscale*( scale*salpha )
631 bcoeff = scale*bcoeff
635 acoefa = abs( acoeff )
636 bcoefa = abs1( bcoeff )
642 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
649 DO 170 jr = 1, je - 1
650 work( jr ) = acoeff*s( jr, je ) - bcoeff*p( jr, je )
654 DO 210 j = je - 1, 1, -1
659 d = acoeff*s( j, j ) - bcoeff*p( j, j )
660 IF( abs1( d ).LE.dmin )
663 IF( abs1( d ).LT.one )
THEN
664 IF( abs1( work( j ) ).GE.bignum*abs1( d ) )
THEN
665 temp = one / abs1( work( j ) )
667 work( jr ) = temp*work( jr )
672 work( j ) =
cladiv( -work( j ), d )
678 IF( abs1( work( j ) ).GT.one )
THEN
679 temp = one / abs1( work( j ) )
680 IF( acoefa*rwork( j )+bcoefa*rwork( n+j ).GE.
683 work( jr ) = temp*work( jr )
688 ca = acoeff*work( j )
689 cb = bcoeff*work( j )
691 work( jr ) = work( jr ) + ca*s( jr, j ) -
700 CALL
cgemv(
'N', n, je, cone, vr, ldvr, work, 1,
701 $ czero, work( n+1 ), 1 )
713 xmax = max( xmax, abs1( work( ( isrc-1 )*n+jr ) ) )
716 IF( xmax.GT.safmin )
THEN
719 vr( jr, ieig ) = temp*work( ( isrc-1 )*n+jr )
725 DO 240 jr = iend + 1, n
726 vr( jr, ieig ) = czero