313 CHARACTER compq, compz, job
314 INTEGER ihi, ilo, info, ldh, ldq, ldt, ldz, lwork, n
317 REAL alphai( * ), alphar( * ), beta( * ),
318 $ h( ldh, * ), q( ldq, * ), t( ldt, * ),
319 $ work( * ), z( ldz, * )
326 REAL half, zero, one, safety
327 parameter ( half = 0.5e+0, zero = 0.0e+0, one = 1.0e+0,
331 LOGICAL ilazr2, ilazro, ilpivt, ilq, ilschr, ilz,
333 INTEGER icompq, icompz, ifirst, ifrstm, iiter, ilast,
334 $ ilastm, in, ischur, istart, j, jc, jch, jiter,
336 REAL a11, a12, a1i, a1r, a21, a22, a2i, a2r, ad11,
337 $ ad11l, ad12, ad12l, ad21, ad21l, ad22, ad22l,
338 $ ad32l, an, anorm, ascale, atol, b11, b1a, b1i,
339 $ b1r, b22, b2a, b2i, b2r, bn, bnorm, bscale,
340 $ btol, c, c11i, c11r, c12, c21, c22i, c22r, cl,
341 $ cq, cr, cz, eshift, s, s1, s1inv, s2, safmax,
342 $ safmin, scale, sl, sqi, sqr, sr, szi, szr, t1,
343 $ tau, temp, temp2, tempi, tempr, u1, u12, u12l,
344 $ u2, ulp, vs, w11, w12, w21, w22, wabs, wi, wr,
360 INTRINSIC abs, max, min,
REAL, sqrt
366 IF(
lsame( job,
'E' ) )
THEN
369 ELSE IF(
lsame( job,
'S' ) )
THEN
376 IF(
lsame( compq,
'N' ) )
THEN
379 ELSE IF(
lsame( compq,
'V' ) )
THEN
382 ELSE IF(
lsame( compq,
'I' ) )
THEN
389 IF(
lsame( compz,
'N' ) )
THEN
392 ELSE IF(
lsame( compz,
'V' ) )
THEN
395 ELSE IF(
lsame( compz,
'I' ) )
THEN
405 work( 1 ) = max( 1, n )
406 lquery = ( lwork.EQ.-1 )
407 IF( ischur.EQ.0 )
THEN
409 ELSE IF( icompq.EQ.0 )
THEN
411 ELSE IF( icompz.EQ.0 )
THEN
413 ELSE IF( n.LT.0 )
THEN
415 ELSE IF( ilo.LT.1 )
THEN
417 ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 )
THEN
419 ELSE IF( ldh.LT.n )
THEN
421 ELSE IF( ldt.LT.n )
THEN
423 ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) )
THEN
425 ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) )
THEN
427 ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery )
THEN
431 CALL xerbla(
'SHGEQZ', -info )
433 ELSE IF( lquery )
THEN
440 work( 1 ) =
REAL( 1 )
447 $
CALL slaset(
'Full', n, n, zero, one, q, ldq )
449 $
CALL slaset(
'Full', n, n, zero, one, z, ldz )
455 safmax = one / safmin
457 anorm =
slanhs(
'F', in, h( ilo, ilo ), ldh, work )
458 bnorm =
slanhs(
'F', in, t( ilo, ilo ), ldt, work )
459 atol = max( safmin, ulp*anorm )
460 btol = max( safmin, ulp*bnorm )
461 ascale = one / max( safmin, anorm )
462 bscale = one / max( safmin, bnorm )
467 IF( t( j, j ).LT.zero )
THEN
470 h( jr, j ) = -h( jr, j )
471 t( jr, j ) = -t( jr, j )
474 h( j, j ) = -h( j, j )
475 t( j, j ) = -t( j, j )
479 z( jr, j ) = -z( jr, j )
483 alphar( j ) = h( j, j )
485 beta( j ) = t( j, j )
518 maxit = 30*( ihi-ilo+1 )
520 DO 360 jiter = 1, maxit
528 IF( ilast.EQ.ilo )
THEN
534 IF( abs( h( ilast, ilast-1 ) ).LE.atol )
THEN
535 h( ilast, ilast-1 ) = zero
540 IF( abs( t( ilast, ilast ) ).LE.btol )
THEN
541 t( ilast, ilast ) = zero
547 DO 60 j = ilast - 1, ilo, -1
554 IF( abs( h( j, j-1 ) ).LE.atol )
THEN
564 IF( abs( t( j, j ) ).LT.btol )
THEN
570 IF( .NOT.ilazro )
THEN
571 temp = abs( h( j, j-1 ) )
572 temp2 = abs( h( j, j ) )
573 tempr = max( temp, temp2 )
574 IF( tempr.LT.one .AND. tempr.NE.zero )
THEN
576 temp2 = temp2 / tempr
578 IF( temp*( ascale*abs( h( j+1, j ) ) ).LE.temp2*
579 $ ( ascale*atol ) )ilazr2 = .true.
588 IF( ilazro .OR. ilazr2 )
THEN
589 DO 40 jch = j, ilast - 1
591 CALL slartg( temp, h( jch+1, jch ), c, s,
593 h( jch+1, jch ) = zero
594 CALL srot( ilastm-jch, h( jch, jch+1 ), ldh,
595 $ h( jch+1, jch+1 ), ldh, c, s )
596 CALL srot( ilastm-jch, t( jch, jch+1 ), ldt,
597 $ t( jch+1, jch+1 ), ldt, c, s )
599 $
CALL srot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
602 $ h( jch, jch-1 ) = h( jch, jch-1 )*c
604 IF( abs( t( jch+1, jch+1 ) ).GE.btol )
THEN
605 IF( jch+1.GE.ilast )
THEN
612 t( jch+1, jch+1 ) = zero
620 DO 50 jch = j, ilast - 1
621 temp = t( jch, jch+1 )
622 CALL slartg( temp, t( jch+1, jch+1 ), c, s,
624 t( jch+1, jch+1 ) = zero
625 IF( jch.LT.ilastm-1 )
626 $
CALL srot( ilastm-jch-1, t( jch, jch+2 ), ldt,
627 $ t( jch+1, jch+2 ), ldt, c, s )
628 CALL srot( ilastm-jch+2, h( jch, jch-1 ), ldh,
629 $ h( jch+1, jch-1 ), ldh, c, s )
631 $
CALL srot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
633 temp = h( jch+1, jch )
634 CALL slartg( temp, h( jch+1, jch-1 ), c, s,
636 h( jch+1, jch-1 ) = zero
637 CALL srot( jch+1-ifrstm, h( ifrstm, jch ), 1,
638 $ h( ifrstm, jch-1 ), 1, c, s )
639 CALL srot( jch-ifrstm, t( ifrstm, jch ), 1,
640 $ t( ifrstm, jch-1 ), 1, c, s )
642 $
CALL srot( n, z( 1, jch ), 1, z( 1, jch-1 ), 1,
647 ELSE IF( ilazro )
THEN
668 temp = h( ilast, ilast )
669 CALL slartg( temp, h( ilast, ilast-1 ), c, s,
670 $ h( ilast, ilast ) )
671 h( ilast, ilast-1 ) = zero
672 CALL srot( ilast-ifrstm, h( ifrstm, ilast ), 1,
673 $ h( ifrstm, ilast-1 ), 1, c, s )
674 CALL srot( ilast-ifrstm, t( ifrstm, ilast ), 1,
675 $ t( ifrstm, ilast-1 ), 1, c, s )
677 $
CALL srot( n, z( 1, ilast ), 1, z( 1, ilast-1 ), 1, c, s )
683 IF( t( ilast, ilast ).LT.zero )
THEN
685 DO 90 j = ifrstm, ilast
686 h( j, ilast ) = -h( j, ilast )
687 t( j, ilast ) = -t( j, ilast )
690 h( ilast, ilast ) = -h( ilast, ilast )
691 t( ilast, ilast ) = -t( ilast, ilast )
695 z( j, ilast ) = -z( j, ilast )
699 alphar( ilast ) = h( ilast, ilast )
700 alphai( ilast ) = zero
701 beta( ilast ) = t( ilast, ilast )
713 IF( .NOT.ilschr )
THEN
715 IF( ifrstm.GT.ilast )
727 IF( .NOT.ilschr )
THEN
737 IF( ( iiter / 10 )*10.EQ.iiter )
THEN
742 IF( (
REAL( maxit )*safmin )*abs( h( ilast, ilast-1 ) ).LT.
743 $ abs( t( ilast-1, ilast-1 ) ) )
THEN
744 eshift = h( ilast, ilast-1 ) /
745 $ t( ilast-1, ilast-1 )
747 eshift = eshift + one / ( safmin*
REAL( MAXIT ) )
758 CALL slag2( h( ilast-1, ilast-1 ), ldh,
759 $ t( ilast-1, ilast-1 ), ldt, safmin*safety, s1,
762 IF ( abs( (wr/s1)*t( ilast, ilast ) - h( ilast, ilast ) )
763 $ .GT. abs( (wr2/s2)*t( ilast, ilast )
764 $ - h( ilast, ilast ) ) )
THEN
772 temp = max( s1, safmin*max( one, abs( wr ), abs( wi ) ) )
779 temp = min( ascale, one )*( half*safmax )
780 IF( s1.GT.temp )
THEN
786 temp = min( bscale, one )*( half*safmax )
787 IF( abs( wr ).GT.temp )
788 $ scale = min( scale, temp / abs( wr ) )
794 DO 120 j = ilast - 1, ifirst + 1, -1
796 temp = abs( s1*h( j, j-1 ) )
797 temp2 = abs( s1*h( j, j )-wr*t( j, j ) )
798 tempr = max( temp, temp2 )
799 IF( tempr.LT.one .AND. tempr.NE.zero )
THEN
801 temp2 = temp2 / tempr
803 IF( abs( ( ascale*h( j+1, j ) )*temp ).LE.( ascale*atol )*
814 temp = s1*h( istart, istart ) - wr*t( istart, istart )
815 temp2 = s1*h( istart+1, istart )
816 CALL slartg( temp, temp2, c, s, tempr )
820 DO 190 j = istart, ilast - 1
821 IF( j.GT.istart )
THEN
823 CALL slartg( temp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
827 DO 140 jc = j, ilastm
828 temp = c*h( j, jc ) + s*h( j+1, jc )
829 h( j+1, jc ) = -s*h( j, jc ) + c*h( j+1, jc )
831 temp2 = c*t( j, jc ) + s*t( j+1, jc )
832 t( j+1, jc ) = -s*t( j, jc ) + c*t( j+1, jc )
837 temp = c*q( jr, j ) + s*q( jr, j+1 )
838 q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
844 CALL slartg( temp, t( j+1, j ), c, s, t( j+1, j+1 ) )
847 DO 160 jr = ifrstm, min( j+2, ilast )
848 temp = c*h( jr, j+1 ) + s*h( jr, j )
849 h( jr, j ) = -s*h( jr, j+1 ) + c*h( jr, j )
852 DO 170 jr = ifrstm, j
853 temp = c*t( jr, j+1 ) + s*t( jr, j )
854 t( jr, j ) = -s*t( jr, j+1 ) + c*t( jr, j )
859 temp = c*z( jr, j+1 ) + s*z( jr, j )
860 z( jr, j ) = -s*z( jr, j+1 ) + c*z( jr, j )
876 IF( ifirst+1.EQ.ilast )
THEN
886 CALL slasv2( t( ilast-1, ilast-1 ), t( ilast-1, ilast ),
887 $ t( ilast, ilast ), b22, b11, sr, cr, sl, cl )
889 IF( b11.LT.zero )
THEN
896 CALL srot( ilastm+1-ifirst, h( ilast-1, ilast-1 ), ldh,
897 $ h( ilast, ilast-1 ), ldh, cl, sl )
898 CALL srot( ilast+1-ifrstm, h( ifrstm, ilast-1 ), 1,
899 $ h( ifrstm, ilast ), 1, cr, sr )
901 IF( ilast.LT.ilastm )
902 $
CALL srot( ilastm-ilast, t( ilast-1, ilast+1 ), ldt,
903 $ t( ilast, ilast+1 ), ldt, cl, sl )
904 IF( ifrstm.LT.ilast-1 )
905 $
CALL srot( ifirst-ifrstm, t( ifrstm, ilast-1 ), 1,
906 $ t( ifrstm, ilast ), 1, cr, sr )
909 $
CALL srot( n, q( 1, ilast-1 ), 1, q( 1, ilast ), 1, cl,
912 $
CALL srot( n, z( 1, ilast-1 ), 1, z( 1, ilast ), 1, cr,
915 t( ilast-1, ilast-1 ) = b11
916 t( ilast-1, ilast ) = zero
917 t( ilast, ilast-1 ) = zero
918 t( ilast, ilast ) = b22
922 IF( b22.LT.zero )
THEN
923 DO 210 j = ifrstm, ilast
924 h( j, ilast ) = -h( j, ilast )
925 t( j, ilast ) = -t( j, ilast )
930 z( j, ilast ) = -z( j, ilast )
940 CALL slag2( h( ilast-1, ilast-1 ), ldh,
941 $ t( ilast-1, ilast-1 ), ldt, safmin*safety, s1,
942 $ temp, wr, temp2, wi )
953 a11 = h( ilast-1, ilast-1 )
954 a21 = h( ilast, ilast-1 )
955 a12 = h( ilast-1, ilast )
956 a22 = h( ilast, ilast )
964 c11r = s1*a11 - wr*b11
968 c22r = s1*a22 - wr*b22
971 IF( abs( c11r )+abs( c11i )+abs( c12 ).GT.abs( c21 )+
972 $ abs( c22r )+abs( c22i ) )
THEN
973 t1 =
slapy3( c12, c11r, c11i )
979 IF( cz.LE.safmin )
THEN
988 szr = -c21*tempr / t1
999 an = abs( a11 ) + abs( a12 ) + abs( a21 ) + abs( a22 )
1000 bn = abs( b11 ) + abs( b22 )
1001 wabs = abs( wr ) + abs( wi )
1002 IF( s1*an.GT.wabs*bn )
THEN
1007 a1r = cz*a11 + szr*a12
1009 a2r = cz*a21 + szr*a22
1012 IF( cq.LE.safmin )
THEN
1019 sqr = tempr*a2r + tempi*a2i
1020 sqi = tempi*a2r - tempr*a2i
1023 t1 =
slapy3( cq, sqr, sqi )
1030 tempr = sqr*szr - sqi*szi
1031 tempi = sqr*szi + sqi*szr
1032 b1r = cq*cz*b11 + tempr*b22
1035 b2r = cq*cz*b22 + tempr*b11
1041 beta( ilast-1 ) = b1a
1043 alphar( ilast-1 ) = ( wr*b1a )*s1inv
1044 alphai( ilast-1 ) = ( wi*b1a )*s1inv
1045 alphar( ilast ) = ( wr*b2a )*s1inv
1046 alphai( ilast ) = -( wi*b2a )*s1inv
1058 IF( .NOT.ilschr )
THEN
1060 IF( ifrstm.GT.ilast )
1078 ad11 = ( ascale*h( ilast-1, ilast-1 ) ) /
1079 $ ( bscale*t( ilast-1, ilast-1 ) )
1080 ad21 = ( ascale*h( ilast, ilast-1 ) ) /
1081 $ ( bscale*t( ilast-1, ilast-1 ) )
1082 ad12 = ( ascale*h( ilast-1, ilast ) ) /
1083 $ ( bscale*t( ilast, ilast ) )
1084 ad22 = ( ascale*h( ilast, ilast ) ) /
1085 $ ( bscale*t( ilast, ilast ) )
1086 u12 = t( ilast-1, ilast ) / t( ilast, ilast )
1087 ad11l = ( ascale*h( ifirst, ifirst ) ) /
1088 $ ( bscale*t( ifirst, ifirst ) )
1089 ad21l = ( ascale*h( ifirst+1, ifirst ) ) /
1090 $ ( bscale*t( ifirst, ifirst ) )
1091 ad12l = ( ascale*h( ifirst, ifirst+1 ) ) /
1092 $ ( bscale*t( ifirst+1, ifirst+1 ) )
1093 ad22l = ( ascale*h( ifirst+1, ifirst+1 ) ) /
1094 $ ( bscale*t( ifirst+1, ifirst+1 ) )
1095 ad32l = ( ascale*h( ifirst+2, ifirst+1 ) ) /
1096 $ ( bscale*t( ifirst+1, ifirst+1 ) )
1097 u12l = t( ifirst, ifirst+1 ) / t( ifirst+1, ifirst+1 )
1099 v( 1 ) = ( ad11-ad11l )*( ad22-ad11l ) - ad12*ad21 +
1100 $ ad21*u12*ad11l + ( ad12l-ad11l*u12l )*ad21l
1101 v( 2 ) = ( ( ad22l-ad11l )-ad21l*u12l-( ad11-ad11l )-
1102 $ ( ad22-ad11l )+ad21*u12 )*ad21l
1103 v( 3 ) = ad32l*ad21l
1107 CALL slarfg( 3, v( 1 ), v( 2 ), 1, tau )
1112 DO 290 j = istart, ilast - 2
1118 IF( j.GT.istart )
THEN
1119 v( 1 ) = h( j, j-1 )
1120 v( 2 ) = h( j+1, j-1 )
1121 v( 3 ) = h( j+2, j-1 )
1123 CALL slarfg( 3, h( j, j-1 ), v( 2 ), 1, tau )
1125 h( j+1, j-1 ) = zero
1126 h( j+2, j-1 ) = zero
1129 DO 230 jc = j, ilastm
1130 temp = tau*( h( j, jc )+v( 2 )*h( j+1, jc )+v( 3 )*
1132 h( j, jc ) = h( j, jc ) - temp
1133 h( j+1, jc ) = h( j+1, jc ) - temp*v( 2 )
1134 h( j+2, jc ) = h( j+2, jc ) - temp*v( 3 )
1135 temp2 = tau*( t( j, jc )+v( 2 )*t( j+1, jc )+v( 3 )*
1137 t( j, jc ) = t( j, jc ) - temp2
1138 t( j+1, jc ) = t( j+1, jc ) - temp2*v( 2 )
1139 t( j+2, jc ) = t( j+2, jc ) - temp2*v( 3 )
1143 temp = tau*( q( jr, j )+v( 2 )*q( jr, j+1 )+v( 3 )*
1145 q( jr, j ) = q( jr, j ) - temp
1146 q( jr, j+1 ) = q( jr, j+1 ) - temp*v( 2 )
1147 q( jr, j+2 ) = q( jr, j+2 ) - temp*v( 3 )
1156 temp = max( abs( t( j+1, j+1 ) ), abs( t( j+1, j+2 ) ) )
1157 temp2 = max( abs( t( j+2, j+1 ) ), abs( t( j+2, j+2 ) ) )
1158 IF( max( temp, temp2 ).LT.safmin )
THEN
1163 ELSE IF( temp.GE.temp2 )
THEN
1181 IF( abs( w12 ).GT.abs( w11 ) )
THEN
1195 w22 = w22 - temp*w12
1201 IF( abs( w22 ).LT.safmin )
THEN
1207 IF( abs( w22 ).LT.abs( u2 ) )
1208 $ scale = abs( w22 / u2 )
1209 IF( abs( w11 ).LT.abs( u1 ) )
1210 $ scale = min( scale, abs( w11 / u1 ) )
1214 u2 = ( scale*u2 ) / w22
1215 u1 = ( scale*u1-w12*u2 ) / w11
1226 t1 = sqrt( scale**2+u1**2+u2**2 )
1227 tau = one + scale / t1
1228 vs = -one / ( scale+t1 )
1235 DO 260 jr = ifrstm, min( j+3, ilast )
1236 temp = tau*( h( jr, j )+v( 2 )*h( jr, j+1 )+v( 3 )*
1238 h( jr, j ) = h( jr, j ) - temp
1239 h( jr, j+1 ) = h( jr, j+1 ) - temp*v( 2 )
1240 h( jr, j+2 ) = h( jr, j+2 ) - temp*v( 3 )
1242 DO 270 jr = ifrstm, j + 2
1243 temp = tau*( t( jr, j )+v( 2 )*t( jr, j+1 )+v( 3 )*
1245 t( jr, j ) = t( jr, j ) - temp
1246 t( jr, j+1 ) = t( jr, j+1 ) - temp*v( 2 )
1247 t( jr, j+2 ) = t( jr, j+2 ) - temp*v( 3 )
1251 temp = tau*( z( jr, j )+v( 2 )*z( jr, j+1 )+v( 3 )*
1253 z( jr, j ) = z( jr, j ) - temp
1254 z( jr, j+1 ) = z( jr, j+1 ) - temp*v( 2 )
1255 z( jr, j+2 ) = z( jr, j+2 ) - temp*v( 3 )
1268 CALL slartg( temp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
1269 h( j+1, j-1 ) = zero
1271 DO 300 jc = j, ilastm
1272 temp = c*h( j, jc ) + s*h( j+1, jc )
1273 h( j+1, jc ) = -s*h( j, jc ) + c*h( j+1, jc )
1275 temp2 = c*t( j, jc ) + s*t( j+1, jc )
1276 t( j+1, jc ) = -s*t( j, jc ) + c*t( j+1, jc )
1281 temp = c*q( jr, j ) + s*q( jr, j+1 )
1282 q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
1289 temp = t( j+1, j+1 )
1290 CALL slartg( temp, t( j+1, j ), c, s, t( j+1, j+1 ) )
1293 DO 320 jr = ifrstm, ilast
1294 temp = c*h( jr, j+1 ) + s*h( jr, j )
1295 h( jr, j ) = -s*h( jr, j+1 ) + c*h( jr, j )
1298 DO 330 jr = ifrstm, ilast - 1
1299 temp = c*t( jr, j+1 ) + s*t( jr, j )
1300 t( jr, j ) = -s*t( jr, j+1 ) + c*t( jr, j )
1305 temp = c*z( jr, j+1 ) + s*z( jr, j )
1306 z( jr, j ) = -s*z( jr, j+1 ) + c*z( jr, j )
1333 DO 410 j = 1, ilo - 1
1334 IF( t( j, j ).LT.zero )
THEN
1337 h( jr, j ) = -h( jr, j )
1338 t( jr, j ) = -t( jr, j )
1341 h( j, j ) = -h( j, j )
1342 t( j, j ) = -t( j, j )
1346 z( jr, j ) = -z( jr, j )
1350 alphar( j ) = h( j, j )
1352 beta( j ) = t( j, j )
1362 work( 1 ) =
REAL( n )
subroutine slarfg(N, ALPHA, X, INCX, TAU)
SLARFG generates an elementary reflector (Householder matrix).
real function slanhs(NORM, N, A, LDA, WORK)
SLANHS returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
real function slapy3(X, Y, Z)
SLAPY3 returns sqrt(x2+y2+z2).
real function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine slasv2(F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL)
SLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
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 ...
real function slamch(CMACH)
SLAMCH
logical function lsame(CA, CB)
LSAME