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
262 EXTERNAL lsame, slamch, cladiv
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 )
374 ulp = slamch(
'Epsilon' )*slamch(
'Base' )
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
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine ctgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
CTGEVC