293 SUBROUTINE stgevc( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
294 $ LDVL, VR, LDVR, MM, M, WORK, INFO )
301 CHARACTER HOWMNY, SIDE
302 INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
306 REAL P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
307 $ vr( ldvr, * ), work( * )
314 REAL ZERO, ONE, SAFETY
315 parameter( zero = 0.0e+0, one = 1.0e+0,
319 LOGICAL COMPL, COMPR, IL2BY2, ILABAD, ILALL, ILBACK,
320 $ ilbbad, ilcomp, ilcplx, lsa, lsb
321 INTEGER I, IBEG, IEIG, IEND, IHWMNY, IINFO, IM, ISIDE,
322 $ j, ja, jc, je, jr, jw, na, nw
323 REAL ACOEF, ACOEFA, ANORM, ASCALE, BCOEFA, BCOEFI,
324 $ bcoefr, big, bignum, bnorm, bscale, cim2a,
325 $ cim2b, cimaga, cimagb, cre2a, cre2b, creala,
326 $ crealb, dmin, safmin, salfar, sbeta, scale,
327 $ small, temp, temp2, temp2i, temp2r, ulp, xmax,
331 REAL BDIAG( 2 ), SUM( 2, 2 ), SUMS( 2, 2 ),
337 EXTERNAL lsame, slamch
343 INTRINSIC abs, max, min
349 IF( lsame( howmny,
'A' ) )
THEN
353 ELSE IF( lsame( howmny,
'S' ) )
THEN
357 ELSE IF( lsame( howmny,
'B' ) )
THEN
366 IF( lsame( side,
'R' ) )
THEN
370 ELSE IF( lsame( side,
'L' ) )
THEN
374 ELSE IF( lsame( side,
'B' ) )
THEN
383 IF( iside.LT.0 )
THEN
385 ELSE IF( ihwmny.LT.0 )
THEN
387 ELSE IF( n.LT.0 )
THEN
389 ELSE IF( lds.LT.max( 1, n ) )
THEN
391 ELSE IF( ldp.LT.max( 1, n ) )
THEN
395 CALL xerbla(
'STGEVC', -info )
401 IF( .NOT.ilall )
THEN
410 IF( s( j+1, j ).NE.zero )
414 IF(
SELECT( j ) .OR.
SELECT( j+1 ) )
430 IF( s( j+1, j ).NE.zero )
THEN
431 IF( p( j, j ).EQ.zero .OR. p( j+1, j+1 ).EQ.zero .OR.
432 $ p( j, j+1 ).NE.zero )ilbbad = .true.
434 IF( s( j+2, j+1 ).NE.zero )
442 ELSE IF( ilbbad )
THEN
444 ELSE IF( compl .AND. ldvl.LT.n .OR. ldvl.LT.1 )
THEN
446 ELSE IF( compr .AND. ldvr.LT.n .OR. ldvr.LT.1 )
THEN
448 ELSE IF( mm.LT.im )
THEN
452 CALL xerbla(
'STGEVC', -info )
464 safmin = slamch(
'Safe minimum' )
466 ulp = slamch(
'Epsilon' )*slamch(
'Base' )
467 small = safmin*n / ulp
469 bignum = one / ( safmin*n )
476 anorm = abs( s( 1, 1 ) )
478 $ anorm = anorm + abs( s( 2, 1 ) )
479 bnorm = abs( p( 1, 1 ) )
486 IF( s( j, j-1 ).EQ.zero )
THEN
492 temp = temp + abs( s( i, j ) )
493 temp2 = temp2 + abs( p( i, j ) )
497 DO 40 i = iend + 1, min( j+1, n )
498 temp = temp + abs( s( i, j ) )
499 temp2 = temp2 + abs( p( i, j ) )
501 anorm = max( anorm, temp )
502 bnorm = max( bnorm, temp2 )
505 ascale = one / max( anorm, safmin )
506 bscale = one / max( bnorm, safmin )
529 IF( s( je+1, je ).NE.zero )
THEN
536 ELSE IF( ilcplx )
THEN
537 ilcomp =
SELECT( je ) .OR.
SELECT( je+1 )
539 ilcomp =
SELECT( je )
547 IF( .NOT.ilcplx )
THEN
548 IF( abs( s( je, je ) ).LE.safmin .AND.
549 $ abs( p( je, je ) ).LE.safmin )
THEN
555 vl( jr, ieig ) = zero
557 vl( ieig, ieig ) = one
565 work( 2*n+jr ) = zero
572 IF( .NOT.ilcplx )
THEN
576 temp = one / max( abs( s( je, je ) )*ascale,
577 $ abs( p( je, je ) )*bscale, safmin )
578 salfar = ( temp*s( je, je ) )*ascale
579 sbeta = ( temp*p( je, je ) )*bscale
581 bcoefr = salfar*bscale
587 lsa = abs( sbeta ).GE.safmin .AND. abs( acoef ).LT.small
588 lsb = abs( salfar ).GE.safmin .AND. abs( bcoefr ).LT.
591 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
593 $ scale = max( scale, ( small / abs( salfar ) )*
594 $ min( bnorm, big ) )
595 IF( lsa .OR. lsb )
THEN
596 scale = min( scale, one /
597 $ ( safmin*max( one, abs( acoef ),
598 $ abs( bcoefr ) ) ) )
600 acoef = ascale*( scale*sbeta )
605 bcoefr = bscale*( scale*salfar )
607 bcoefr = scale*bcoefr
610 acoefa = abs( acoef )
611 bcoefa = abs( bcoefr )
621 CALL slag2( s( je, je ), lds, p( je, je ), ldp,
622 $ safmin*safety, acoef, temp, bcoefr, temp2,
625 IF( bcoefi.EQ.zero )
THEN
632 acoefa = abs( acoef )
633 bcoefa = abs( bcoefr ) + abs( bcoefi )
635 IF( acoefa*ulp.LT.safmin .AND. acoefa.GE.safmin )
636 $ scale = ( safmin / ulp ) / acoefa
637 IF( bcoefa*ulp.LT.safmin .AND. bcoefa.GE.safmin )
638 $ scale = max( scale, ( safmin / ulp ) / bcoefa )
639 IF( safmin*acoefa.GT.ascale )
640 $ scale = ascale / ( safmin*acoefa )
641 IF( safmin*bcoefa.GT.bscale )
642 $ scale = min( scale, bscale / ( safmin*bcoefa ) )
643 IF( scale.NE.one )
THEN
645 acoefa = abs( acoef )
646 bcoefr = scale*bcoefr
647 bcoefi = scale*bcoefi
648 bcoefa = abs( bcoefr ) + abs( bcoefi )
653 temp = acoef*s( je+1, je )
654 temp2r = acoef*s( je, je ) - bcoefr*p( je, je )
655 temp2i = -bcoefi*p( je, je )
656 IF( abs( temp ).GT.abs( temp2r )+abs( temp2i ) )
THEN
658 work( 3*n+je ) = zero
659 work( 2*n+je+1 ) = -temp2r / temp
660 work( 3*n+je+1 ) = -temp2i / temp
662 work( 2*n+je+1 ) = one
663 work( 3*n+je+1 ) = zero
664 temp = acoef*s( je, je+1 )
665 work( 2*n+je ) = ( bcoefr*p( je+1, je+1 )-acoef*
666 $ s( je+1, je+1 ) ) / temp
667 work( 3*n+je ) = bcoefi*p( je+1, je+1 ) / temp
669 xmax = max( abs( work( 2*n+je ) )+abs( work( 3*n+je ) ),
670 $ abs( work( 2*n+je+1 ) )+abs( work( 3*n+je+1 ) ) )
673 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
683 DO 160 j = je + nw, n
690 bdiag( 1 ) = p( j, j )
692 IF( s( j+1, j ).NE.zero )
THEN
694 bdiag( 2 ) = p( j+1, j+1 )
701 xscale = one / max( one, xmax )
702 temp = max( work( j ), work( n+j ),
703 $ acoefa*work( j )+bcoefa*work( n+j ) )
705 $ temp = max( temp, work( j+1 ), work( n+j+1 ),
706 $ acoefa*work( j+1 )+bcoefa*work( n+j+1 ) )
707 IF( temp.GT.bignum*xscale )
THEN
710 work( ( jw+2 )*n+jr ) = xscale*
711 $ work( ( jw+2 )*n+jr )
735 sums( ja, jw ) = zero
736 sump( ja, jw ) = zero
738 DO 100 jr = je, j - 1
739 sums( ja, jw ) = sums( ja, jw ) +
741 $ work( ( jw+1 )*n+jr )
742 sump( ja, jw ) = sump( ja, jw ) +
744 $ work( ( jw+1 )*n+jr )
751 sum( ja, 1 ) = -acoef*sums( ja, 1 ) +
752 $ bcoefr*sump( ja, 1 ) -
753 $ bcoefi*sump( ja, 2 )
754 sum( ja, 2 ) = -acoef*sums( ja, 2 ) +
755 $ bcoefr*sump( ja, 2 ) +
756 $ bcoefi*sump( ja, 1 )
758 sum( ja, 1 ) = -acoef*sums( ja, 1 ) +
759 $ bcoefr*sump( ja, 1 )
767 CALL slaln2( .true., na, nw, dmin, acoef, s( j, j ), lds,
768 $ bdiag( 1 ), bdiag( 2 ), sum, 2, bcoefr,
769 $ bcoefi, work( 2*n+j ), n, scale, temp,
771 IF( scale.LT.one )
THEN
772 DO 150 jw = 0, nw - 1
773 DO 140 jr = je, j - 1
774 work( ( jw+2 )*n+jr ) = scale*
775 $ work( ( jw+2 )*n+jr )
780 xmax = max( xmax, temp )
788 DO 170 jw = 0, nw - 1
789 CALL sgemv(
'N', n, n+1-je, one, vl( 1, je ), ldvl,
790 $ work( ( jw+2 )*n+je ), 1, zero,
791 $ work( ( jw+4 )*n+1 ), 1 )
793 CALL slacpy(
' ', n, nw, work( 4*n+1 ), n, vl( 1, je ),
797 CALL slacpy(
' ', n, nw, work( 2*n+1 ), n, vl( 1, ieig ),
807 xmax = max( xmax, abs( vl( j, ieig ) )+
808 $ abs( vl( j, ieig+1 ) ) )
812 xmax = max( xmax, abs( vl( j, ieig ) ) )
816 IF( xmax.GT.safmin )
THEN
819 DO 210 jw = 0, nw - 1
821 vl( jr, ieig+jw ) = xscale*vl( jr, ieig+jw )
854 IF( s( je, je-1 ).NE.zero )
THEN
861 ELSE IF( ilcplx )
THEN
862 ilcomp =
SELECT( je ) .OR.
SELECT( je-1 )
864 ilcomp =
SELECT( je )
872 IF( .NOT.ilcplx )
THEN
873 IF( abs( s( je, je ) ).LE.safmin .AND.
874 $ abs( p( je, je ) ).LE.safmin )
THEN
880 vr( jr, ieig ) = zero
882 vr( ieig, ieig ) = one
889 DO 250 jw = 0, nw - 1
891 work( ( jw+2 )*n+jr ) = zero
899 IF( .NOT.ilcplx )
THEN
903 temp = one / max( abs( s( je, je ) )*ascale,
904 $ abs( p( je, je ) )*bscale, safmin )
905 salfar = ( temp*s( je, je ) )*ascale
906 sbeta = ( temp*p( je, je ) )*bscale
908 bcoefr = salfar*bscale
914 lsa = abs( sbeta ).GE.safmin .AND. abs( acoef ).LT.small
915 lsb = abs( salfar ).GE.safmin .AND. abs( bcoefr ).LT.
918 $ scale = ( small / abs( sbeta ) )*min( anorm, big )
920 $ scale = max( scale, ( small / abs( salfar ) )*
921 $ min( bnorm, big ) )
922 IF( lsa .OR. lsb )
THEN
923 scale = min( scale, one /
924 $ ( safmin*max( one, abs( acoef ),
925 $ abs( bcoefr ) ) ) )
927 acoef = ascale*( scale*sbeta )
932 bcoefr = bscale*( scale*salfar )
934 bcoefr = scale*bcoefr
937 acoefa = abs( acoef )
938 bcoefa = abs( bcoefr )
948 DO 260 jr = 1, je - 1
949 work( 2*n+jr ) = bcoefr*p( jr, je ) -
956 CALL slag2( s( je-1, je-1 ), lds, p( je-1, je-1 ), ldp,
957 $ safmin*safety, acoef, temp, bcoefr, temp2,
959 IF( bcoefi.EQ.zero )
THEN
966 acoefa = abs( acoef )
967 bcoefa = abs( bcoefr ) + abs( bcoefi )
969 IF( acoefa*ulp.LT.safmin .AND. acoefa.GE.safmin )
970 $ scale = ( safmin / ulp ) / acoefa
971 IF( bcoefa*ulp.LT.safmin .AND. bcoefa.GE.safmin )
972 $ scale = max( scale, ( safmin / ulp ) / bcoefa )
973 IF( safmin*acoefa.GT.ascale )
974 $ scale = ascale / ( safmin*acoefa )
975 IF( safmin*bcoefa.GT.bscale )
976 $ scale = min( scale, bscale / ( safmin*bcoefa ) )
977 IF( scale.NE.one )
THEN
979 acoefa = abs( acoef )
980 bcoefr = scale*bcoefr
981 bcoefi = scale*bcoefi
982 bcoefa = abs( bcoefr ) + abs( bcoefi )
988 temp = acoef*s( je, je-1 )
989 temp2r = acoef*s( je, je ) - bcoefr*p( je, je )
990 temp2i = -bcoefi*p( je, je )
991 IF( abs( temp ).GE.abs( temp2r )+abs( temp2i ) )
THEN
993 work( 3*n+je ) = zero
994 work( 2*n+je-1 ) = -temp2r / temp
995 work( 3*n+je-1 ) = -temp2i / temp
997 work( 2*n+je-1 ) = one
998 work( 3*n+je-1 ) = zero
999 temp = acoef*s( je-1, je )
1000 work( 2*n+je ) = ( bcoefr*p( je-1, je-1 )-acoef*
1001 $ s( je-1, je-1 ) ) / temp
1002 work( 3*n+je ) = bcoefi*p( je-1, je-1 ) / temp
1005 xmax = max( abs( work( 2*n+je ) )+abs( work( 3*n+je ) ),
1006 $ abs( work( 2*n+je-1 ) )+abs( work( 3*n+je-1 ) ) )
1011 creala = acoef*work( 2*n+je-1 )
1012 cimaga = acoef*work( 3*n+je-1 )
1013 crealb = bcoefr*work( 2*n+je-1 ) -
1014 $ bcoefi*work( 3*n+je-1 )
1015 cimagb = bcoefi*work( 2*n+je-1 ) +
1016 $ bcoefr*work( 3*n+je-1 )
1017 cre2a = acoef*work( 2*n+je )
1018 cim2a = acoef*work( 3*n+je )
1019 cre2b = bcoefr*work( 2*n+je ) - bcoefi*work( 3*n+je )
1020 cim2b = bcoefi*work( 2*n+je ) + bcoefr*work( 3*n+je )
1021 DO 270 jr = 1, je - 2
1022 work( 2*n+jr ) = -creala*s( jr, je-1 ) +
1023 $ crealb*p( jr, je-1 ) -
1024 $ cre2a*s( jr, je ) + cre2b*p( jr, je )
1025 work( 3*n+jr ) = -cimaga*s( jr, je-1 ) +
1026 $ cimagb*p( jr, je-1 ) -
1027 $ cim2a*s( jr, je ) + cim2b*p( jr, je )
1031 dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
1036 DO 370 j = je - nw, 1, -1
1041 IF( .NOT.il2by2 .AND. j.GT.1 )
THEN
1042 IF( s( j, j-1 ).NE.zero )
THEN
1047 bdiag( 1 ) = p( j, j )
1050 bdiag( 2 ) = p( j+1, j+1 )
1057 CALL slaln2( .false., na, nw, dmin, acoef, s( j, j ),
1058 $ lds, bdiag( 1 ), bdiag( 2 ), work( 2*n+j ),
1059 $ n, bcoefr, bcoefi, sum, 2, scale, temp,
1061 IF( scale.LT.one )
THEN
1063 DO 290 jw = 0, nw - 1
1065 work( ( jw+2 )*n+jr ) = scale*
1066 $ work( ( jw+2 )*n+jr )
1070 xmax = max( scale*xmax, temp )
1074 work( ( jw+1 )*n+j+ja-1 ) = sum( ja, jw )
1084 xscale = one / max( one, xmax )
1085 temp = acoefa*work( j ) + bcoefa*work( n+j )
1087 $ temp = max( temp, acoefa*work( j+1 )+bcoefa*
1089 temp = max( temp, acoefa, bcoefa )
1090 IF( temp.GT.bignum*xscale )
THEN
1092 DO 330 jw = 0, nw - 1
1094 work( ( jw+2 )*n+jr ) = xscale*
1095 $ work( ( jw+2 )*n+jr )
1108 creala = acoef*work( 2*n+j+ja-1 )
1109 cimaga = acoef*work( 3*n+j+ja-1 )
1110 crealb = bcoefr*work( 2*n+j+ja-1 ) -
1111 $ bcoefi*work( 3*n+j+ja-1 )
1112 cimagb = bcoefi*work( 2*n+j+ja-1 ) +
1113 $ bcoefr*work( 3*n+j+ja-1 )
1114 DO 340 jr = 1, j - 1
1115 work( 2*n+jr ) = work( 2*n+jr ) -
1116 $ creala*s( jr, j+ja-1 ) +
1117 $ crealb*p( jr, j+ja-1 )
1118 work( 3*n+jr ) = work( 3*n+jr ) -
1119 $ cimaga*s( jr, j+ja-1 ) +
1120 $ cimagb*p( jr, j+ja-1 )
1123 creala = acoef*work( 2*n+j+ja-1 )
1124 crealb = bcoefr*work( 2*n+j+ja-1 )
1125 DO 350 jr = 1, j - 1
1126 work( 2*n+jr ) = work( 2*n+jr ) -
1127 $ creala*s( jr, j+ja-1 ) +
1128 $ crealb*p( jr, j+ja-1 )
1143 DO 410 jw = 0, nw - 1
1145 work( ( jw+4 )*n+jr ) = work( ( jw+2 )*n+1 )*
1155 work( ( jw+4 )*n+jr ) = work( ( jw+4 )*n+jr ) +
1156 $ work( ( jw+2 )*n+jc )*vr( jr, jc )
1161 DO 430 jw = 0, nw - 1
1163 vr( jr, ieig+jw ) = work( ( jw+4 )*n+jr )
1169 DO 450 jw = 0, nw - 1
1171 vr( jr, ieig+jw ) = work( ( jw+2 )*n+jr )
1183 xmax = max( xmax, abs( vr( j, ieig ) )+
1184 $ abs( vr( j, ieig+1 ) ) )
1188 xmax = max( xmax, abs( vr( j, ieig ) ) )
1192 IF( xmax.GT.safmin )
THEN
1194 DO 490 jw = 0, nw - 1
1196 vr( jr, ieig+jw ) = xscale*vr( jr, ieig+jw )
subroutine xerbla(srname, info)
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
subroutine slag2(a, lda, b, ldb, safmin, scale1, scale2, wr1, wr2, wi)
SLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
subroutine slaln2(ltrans, na, nw, smin, ca, a, lda, d1, d2, b, ldb, wr, wi, x, ldx, scale, xnorm, info)
SLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.
subroutine stgevc(side, howmny, select, n, s, lds, p, ldp, vl, ldvl, vr, ldvr, mm, m, work, info)
STGEVC