217 SUBROUTINE ctgevc( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
218 $ LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
225 CHARACTER HOWMNY, SIDE
226 INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
231 COMPLEX P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
232 $ vr( ldvr, * ), work( * )
240 parameter( zero = 0.0e+0, one = 1.0e+0 )
242 parameter( czero = ( 0.0e+0, 0.0e+0 ),
243 $ cone = ( 1.0e+0, 0.0e+0 ) )
246 LOGICAL COMPL, COMPR, ILALL, ILBACK, ILBBAD, ILCOMP,
248 INTEGER I, IBEG, IEIG, IEND, IHWMNY, IM, ISIDE, ISRC,
250 REAL ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG,
251 $ bignum, bnorm, bscale, dmin, safmin, sbeta,
252 $ scale, small, temp, ulp, xmax
253 COMPLEX BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X
259 EXTERNAL lsame, slamch, cladiv
265 INTRINSIC abs, aimag, cmplx, conjg, max, min, real
271 abs1( x ) = abs( real( x ) ) + abs( aimag( x ) )
277 IF( lsame( howmny,
'A' ) )
THEN
281 ELSE IF( lsame( howmny,
'S' ) )
THEN
285 ELSE IF( lsame( howmny,
'B' ) )
THEN
293 IF( lsame( side,
'R' ) )
THEN
297 ELSE IF( lsame( side,
'L' ) )
THEN
301 ELSE IF( lsame( side,
'B' ) )
THEN
310 IF( iside.LT.0 )
THEN
312 ELSE IF( ihwmny.LT.0 )
THEN
314 ELSE IF( n.LT.0 )
THEN
316 ELSE IF( lds.LT.max( 1, n ) )
THEN
318 ELSE IF( ldp.LT.max( 1, n ) )
THEN
322 CALL xerbla(
'CTGEVC', -info )
328 IF( .NOT.ilall )
THEN
342 IF( aimag( p( j, j ) ).NE.zero )
348 ELSE IF( compl .AND. ldvl.LT.n .OR. ldvl.LT.1 )
THEN
350 ELSE IF( compr .AND. ldvr.LT.n .OR. ldvr.LT.1 )
THEN
352 ELSE IF( mm.LT.im )
THEN
356 CALL xerbla(
'CTGEVC', -info )
368 safmin = slamch(
'Safe minimum' )
370 CALL slabad( safmin, big )
371 ulp = slamch(
'Epsilon' )*slamch(
'Base' )
372 small = safmin*n / ulp
374 bignum = one / ( safmin*n )
380 anorm = abs1( s( 1, 1 ) )
381 bnorm = abs1( p( 1, 1 ) )
388 rwork( j ) = rwork( j ) + abs1( s( i, j ) )
389 rwork( n+j ) = rwork( n+j ) + abs1( p( i, j ) )
391 anorm = max( anorm, rwork( j )+abs1( s( j, j ) ) )
392 bnorm = max( bnorm, rwork( n+j )+abs1( p( j, j ) ) )
395 ascale = one / max( anorm, safmin )
396 bscale = one / max( bnorm, safmin )
409 ilcomp =
SELECT( je )
414 IF( abs1( s( je, je ) ).LE.safmin .AND.
415 $ abs( real( p( je, je ) ) ).LE.safmin )
THEN
420 vl( jr, ieig ) = czero
422 vl( ieig, ieig ) = cone
431 temp = one / max( abs1( s( je, je ) )*ascale,
432 $ abs( real( p( je, je ) ) )*bscale, safmin )
433 salpha = ( temp*s( je, je ) )*ascale
434 sbeta = ( temp*real( p( je, je ) ) )*bscale
435 acoeff = sbeta*ascale
436 bcoeff = salpha*bscale
440 lsa = abs( sbeta ).GE.safmin .AND. abs( acoeff ).LT.small
441 lsb = abs1( salpha ).GE.safmin .AND. abs1( bcoeff ).LT.
446 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
448 $ scale = max( scale, ( small / abs1( salpha ) )*
449 $ min( bnorm, big ) )
450 IF( lsa .OR. lsb )
THEN
451 scale = min( scale, one /
452 $ ( safmin*max( one, abs( acoeff ),
453 $ abs1( bcoeff ) ) ) )
455 acoeff = ascale*( scale*sbeta )
457 acoeff = scale*acoeff
460 bcoeff = bscale*( scale*salpha )
462 bcoeff = scale*bcoeff
466 acoefa = abs( acoeff )
467 bcoefa = abs1( bcoeff )
473 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
490 IF( acoefa*rwork( j )+bcoefa*rwork( n+j ).GT.bignum*
493 work( jr ) = temp*work( jr )
501 suma = suma + conjg( s( jr, j ) )*work( jr )
502 sumb = sumb + conjg( p( jr, j ) )*work( jr )
504 sum = acoeff*suma - conjg( bcoeff )*sumb
510 d = conjg( acoeff*s( j, j )-bcoeff*p( j, j ) )
511 IF( abs1( d ).LE.dmin )
514 IF( abs1( d ).LT.one )
THEN
515 IF( abs1( sum ).GE.bignum*abs1( d ) )
THEN
516 temp = one / abs1( sum )
518 work( jr ) = temp*work( jr )
524 work( j ) = cladiv( -sum, d )
525 xmax = max( xmax, abs1( work( j ) ) )
531 CALL cgemv(
'N', n, n+1-je, cone, vl( 1, je ), ldvl,
532 $ work( je ), 1, czero, work( n+1 ), 1 )
544 xmax = max( xmax, abs1( work( ( isrc-1 )*n+jr ) ) )
547 IF( xmax.GT.safmin )
THEN
550 vl( jr, ieig ) = temp*work( ( isrc-1 )*n+jr )
556 DO 130 jr = 1, ibeg - 1
557 vl( jr, ieig ) = czero
575 ilcomp =
SELECT( je )
580 IF( abs1( s( je, je ) ).LE.safmin .AND.
581 $ abs( real( p( je, je ) ) ).LE.safmin )
THEN
586 vr( jr, ieig ) = czero
588 vr( ieig, ieig ) = cone
597 temp = one / max( abs1( s( je, je ) )*ascale,
598 $ abs( real( p( je, je ) ) )*bscale, safmin )
599 salpha = ( temp*s( je, je ) )*ascale
600 sbeta = ( temp*real( p( je, je ) ) )*bscale
601 acoeff = sbeta*ascale
602 bcoeff = salpha*bscale
606 lsa = abs( sbeta ).GE.safmin .AND. abs( acoeff ).LT.small
607 lsb = abs1( salpha ).GE.safmin .AND. abs1( bcoeff ).LT.
612 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
614 $ scale = max( scale, ( small / abs1( salpha ) )*
615 $ min( bnorm, big ) )
616 IF( lsa .OR. lsb )
THEN
617 scale = min( scale, one /
618 $ ( safmin*max( one, abs( acoeff ),
619 $ abs1( bcoeff ) ) ) )
621 acoeff = ascale*( scale*sbeta )
623 acoeff = scale*acoeff
626 bcoeff = bscale*( scale*salpha )
628 bcoeff = scale*bcoeff
632 acoefa = abs( acoeff )
633 bcoefa = abs1( bcoeff )
639 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
646 DO 170 jr = 1, je - 1
647 work( jr ) = acoeff*s( jr, je ) - bcoeff*p( jr, je )
651 DO 210 j = je - 1, 1, -1
656 d = acoeff*s( j, j ) - bcoeff*p( j, j )
657 IF( abs1( d ).LE.dmin )
660 IF( abs1( d ).LT.one )
THEN
661 IF( abs1( work( j ) ).GE.bignum*abs1( d ) )
THEN
662 temp = one / abs1( work( j ) )
664 work( jr ) = temp*work( jr )
669 work( j ) = cladiv( -work( j ), d )
675 IF( abs1( work( j ) ).GT.one )
THEN
676 temp = one / abs1( work( j ) )
677 IF( acoefa*rwork( j )+bcoefa*rwork( n+j ).GE.
680 work( jr ) = temp*work( jr )
685 ca = acoeff*work( j )
686 cb = bcoeff*work( j )
688 work( jr ) = work( jr ) + ca*s( jr, j ) -
697 CALL cgemv(
'N', n, je, cone, vr, ldvr, work, 1,
698 $ czero, work( n+1 ), 1 )
710 xmax = max( xmax, abs1( work( ( isrc-1 )*n+jr ) ) )
713 IF( xmax.GT.safmin )
THEN
716 vr( jr, ieig ) = temp*work( ( isrc-1 )*n+jr )
722 DO 240 jr = iend + 1, n
723 vr( jr, ieig ) = czero
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
subroutine ctgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
CTGEVC