281 SUBROUTINE zhgeqz( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
282 $ ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK,
290 CHARACTER COMPQ, COMPZ, JOB
291 INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
294 DOUBLE PRECISION RWORK( * )
295 COMPLEX*16 ALPHA( * ), BETA( * ), H( LDH, * ),
296 $ q( ldq, * ), t( ldt, * ), work( * ),
303 COMPLEX*16 CZERO, CONE
304 PARAMETER ( CZERO = ( 0.0d+0, 0.0d+0 ),
305 $ cone = ( 1.0d+0, 0.0d+0 ) )
306 DOUBLE PRECISION ZERO, ONE
307 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0 )
308 DOUBLE PRECISION HALF
309 parameter( half = 0.5d+0 )
312 LOGICAL ILAZR2, ILAZRO, ILQ, ILSCHR, ILZ, LQUERY
313 INTEGER ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
314 $ ilastm, in, ischur, istart, j, jc, jch, jiter,
316 DOUBLE PRECISION ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL,
317 $ C, SAFMIN, TEMP, TEMP2, TEMPR, ULP
318 COMPLEX*16 ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2,
319 $ CTEMP3, ESHIFT, S, SHIFT, SIGNBC,
325 DOUBLE PRECISION DLAMCH, ZLANHS
326 EXTERNAL zladiv, lsame, dlamch, zlanhs
332 INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, min,
336 DOUBLE PRECISION ABS1
339 abs1( x ) = abs( dble( x ) ) + abs( dimag( x ) )
345 IF( lsame( job,
'E' ) )
THEN
348 ELSE IF( lsame( job,
'S' ) )
THEN
356 IF( lsame( compq,
'N' ) )
THEN
359 ELSE IF( lsame( compq,
'V' ) )
THEN
362 ELSE IF( lsame( compq,
'I' ) )
THEN
370 IF( lsame( compz,
'N' ) )
THEN
373 ELSE IF( lsame( compz,
'V' ) )
THEN
376 ELSE IF( lsame( compz,
'I' ) )
THEN
387 work( 1 ) = max( 1, n )
388 lquery = ( lwork.EQ.-1 )
389 IF( ischur.EQ.0 )
THEN
391 ELSE IF( icompq.EQ.0 )
THEN
393 ELSE IF( icompz.EQ.0 )
THEN
395 ELSE IF( n.LT.0 )
THEN
397 ELSE IF( ilo.LT.1 )
THEN
399 ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 )
THEN
401 ELSE IF( ldh.LT.n )
THEN
403 ELSE IF( ldt.LT.n )
THEN
405 ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) )
THEN
407 ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) )
THEN
409 ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery )
THEN
413 CALL xerbla(
'ZHGEQZ', -info )
415 ELSE IF( lquery )
THEN
423 work( 1 ) = dcmplx( 1 )
430 $
CALL zlaset(
'Full', n, n, czero, cone, q, ldq )
432 $
CALL zlaset(
'Full', n, n, czero, cone, z, ldz )
437 safmin = dlamch(
'S' )
438 ulp = dlamch(
'E' )*dlamch(
'B' )
439 anorm = zlanhs(
'F', in, h( ilo, ilo ), ldh, rwork )
440 bnorm = zlanhs(
'F', in, t( ilo, ilo ), ldt, rwork )
441 atol = max( safmin, ulp*anorm )
442 btol = max( safmin, ulp*bnorm )
443 ascale = one / max( safmin, anorm )
444 bscale = one / max( safmin, bnorm )
450 absb = abs( t( j, j ) )
451 IF( absb.GT.safmin )
THEN
452 signbc = dconjg( t( j, j ) / absb )
455 CALL zscal( j-1, signbc, t( 1, j ), 1 )
456 CALL zscal( j, signbc, h( 1, j ), 1 )
458 CALL zscal( 1, signbc, h( j, j ), 1 )
461 $
CALL zscal( n, signbc, z( 1, j ), 1 )
465 alpha( j ) = h( j, j )
466 beta( j ) = t( j, j )
499 maxit = 30*( ihi-ilo+1 )
501 DO 170 jiter = 1, maxit
516 IF( ilast.EQ.ilo )
THEN
519 IF( abs1( h( ilast, ilast-1 ) ).LE.max( safmin, ulp*(
520 $ abs1( h( ilast, ilast ) ) + abs1( h( ilast-1, ilast-1 )
522 h( ilast, ilast-1 ) = czero
527 IF( abs( t( ilast, ilast ) ).LE.btol )
THEN
528 t( ilast, ilast ) = czero
534 DO 40 j = ilast - 1, ilo, -1
541 IF( abs1( h( j, j-1 ) ).LE.max( safmin, ulp*(
542 $ abs1( h( j, j ) ) + abs1( h( j-1, j-1 ) )
553 IF( abs( t( j, j ) ).LT.btol )
THEN
559 IF( .NOT.ilazro )
THEN
560 IF( abs1( h( j, j-1 ) )*( ascale*abs1( h( j+1,
561 $ j ) ) ).LE.abs1( h( j, j ) )*( ascale*atol ) )
571 IF( ilazro .OR. ilazr2 )
THEN
572 DO 20 jch = j, ilast - 1
573 ctemp = h( jch, jch )
574 CALL zlartg( ctemp, h( jch+1, jch ), c, s,
576 h( jch+1, jch ) = czero
577 CALL zrot( ilastm-jch, h( jch, jch+1 ), ldh,
578 $ h( jch+1, jch+1 ), ldh, c, s )
579 CALL zrot( ilastm-jch, t( jch, jch+1 ), ldt,
580 $ t( jch+1, jch+1 ), ldt, c, s )
582 $
CALL zrot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
585 $ h( jch, jch-1 ) = h( jch, jch-1 )*c
587 IF( abs1( t( jch+1, jch+1 ) ).GE.btol )
THEN
588 IF( jch+1.GE.ilast )
THEN
595 t( jch+1, jch+1 ) = czero
603 DO 30 jch = j, ilast - 1
604 ctemp = t( jch, jch+1 )
605 CALL zlartg( ctemp, t( jch+1, jch+1 ), c, s,
607 t( jch+1, jch+1 ) = czero
608 IF( jch.LT.ilastm-1 )
609 $
CALL zrot( ilastm-jch-1, t( jch, jch+2 ), ldt,
610 $ t( jch+1, jch+2 ), ldt, c, s )
611 CALL zrot( ilastm-jch+2, h( jch, jch-1 ), ldh,
612 $ h( jch+1, jch-1 ), ldh, c, s )
614 $
CALL zrot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
616 ctemp = h( jch+1, jch )
617 CALL zlartg( ctemp, h( jch+1, jch-1 ), c, s,
619 h( jch+1, jch-1 ) = czero
620 CALL zrot( jch+1-ifrstm, h( ifrstm, jch ), 1,
621 $ h( ifrstm, jch-1 ), 1, c, s )
622 CALL zrot( jch-ifrstm, t( ifrstm, jch ), 1,
623 $ t( ifrstm, jch-1 ), 1, c, s )
625 $
CALL zrot( n, z( 1, jch ), 1, z( 1, jch-1 ), 1,
630 ELSE IF( ilazro )
THEN
651 ctemp = h( ilast, ilast )
652 CALL zlartg( ctemp, h( ilast, ilast-1 ), c, s,
653 $ h( ilast, ilast ) )
654 h( ilast, ilast-1 ) = czero
655 CALL zrot( ilast-ifrstm, h( ifrstm, ilast ), 1,
656 $ h( ifrstm, ilast-1 ), 1, c, s )
657 CALL zrot( ilast-ifrstm, t( ifrstm, ilast ), 1,
658 $ t( ifrstm, ilast-1 ), 1, c, s )
660 $
CALL zrot( n, z( 1, ilast ), 1, z( 1, ilast-1 ), 1, c, s )
665 absb = abs( t( ilast, ilast ) )
666 IF( absb.GT.safmin )
THEN
667 signbc = dconjg( t( ilast, ilast ) / absb )
668 t( ilast, ilast ) = absb
670 CALL zscal( ilast-ifrstm, signbc, t( ifrstm, ilast ), 1 )
671 CALL zscal( ilast+1-ifrstm, signbc, h( ifrstm, ilast ),
674 CALL zscal( 1, signbc, h( ilast, ilast ), 1 )
677 $
CALL zscal( n, signbc, z( 1, ilast ), 1 )
679 t( ilast, ilast ) = czero
681 alpha( ilast ) = h( ilast, ilast )
682 beta( ilast ) = t( ilast, ilast )
694 IF( .NOT.ilschr )
THEN
696 IF( ifrstm.GT.ilast )
708 IF( .NOT.ilschr )
THEN
718 IF( ( iiter / 10 )*10.NE.iiter )
THEN
727 u12 = ( bscale*t( ilast-1, ilast ) ) /
728 $ ( bscale*t( ilast, ilast ) )
729 ad11 = ( ascale*h( ilast-1, ilast-1 ) ) /
730 $ ( bscale*t( ilast-1, ilast-1 ) )
731 ad21 = ( ascale*h( ilast, ilast-1 ) ) /
732 $ ( bscale*t( ilast-1, ilast-1 ) )
733 ad12 = ( ascale*h( ilast-1, ilast ) ) /
734 $ ( bscale*t( ilast, ilast ) )
735 ad22 = ( ascale*h( ilast, ilast ) ) /
736 $ ( bscale*t( ilast, ilast ) )
737 abi22 = ad22 - u12*ad21
738 abi12 = ad12 - u12*ad11
741 ctemp = sqrt( abi12 )*sqrt( ad21 )
743 IF( ctemp.NE.zero )
THEN
744 x = half*( ad11-shift )
746 temp = max( temp, abs1( x ) )
747 y = temp*sqrt( ( x / temp )**2+( ctemp / temp )**2 )
748 IF( temp2.GT.zero )
THEN
749 IF( dble( x / temp2 )*dble( y )+
750 $ dimag( x / temp2 )*dimag( y ).LT.zero )y = -y
752 shift = shift - ctemp*zladiv( ctemp, ( x+y ) )
758 IF( ( iiter / 20 )*20.EQ.iiter .AND.
759 $ bscale*abs1(t( ilast, ilast )).GT.safmin )
THEN
760 eshift = eshift + ( ascale*h( ilast,
761 $ ilast ) )/( bscale*t( ilast, ilast ) )
763 eshift = eshift + ( ascale*h( ilast,
764 $ ilast-1 ) )/( bscale*t( ilast-1, ilast-1 ) )
771 DO 80 j = ilast - 1, ifirst + 1, -1
773 ctemp = ascale*h( j, j ) - shift*( bscale*t( j, j ) )
775 temp2 = ascale*abs1( h( j+1, j ) )
776 tempr = max( temp, temp2 )
777 IF( tempr.LT.one .AND. tempr.NE.zero )
THEN
779 temp2 = temp2 / tempr
781 IF( abs1( h( j, j-1 ) )*temp2.LE.temp*atol )
786 ctemp = ascale*h( ifirst, ifirst ) -
787 $ shift*( bscale*t( ifirst, ifirst ) )
794 ctemp2 = ascale*h( istart+1, istart )
795 CALL zlartg( ctemp, ctemp2, c, s, ctemp3 )
799 DO 150 j = istart, ilast - 1
800 IF( j.GT.istart )
THEN
802 CALL zlartg( ctemp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
803 h( j+1, j-1 ) = czero
806 DO 100 jc = j, ilastm
807 ctemp = c*h( j, jc ) + s*h( j+1, jc )
808 h( j+1, jc ) = -dconjg( s )*h( j, jc ) + c*h( j+1, jc )
810 ctemp2 = c*t( j, jc ) + s*t( j+1, jc )
811 t( j+1, jc ) = -dconjg( s )*t( j, jc ) + c*t( j+1, jc )
816 ctemp = c*q( jr, j ) + dconjg( s )*q( jr, j+1 )
817 q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
822 ctemp = t( j+1, j+1 )
823 CALL zlartg( ctemp, t( j+1, j ), c, s, t( j+1, j+1 ) )
826 DO 120 jr = ifrstm, min( j+2, ilast )
827 ctemp = c*h( jr, j+1 ) + s*h( jr, j )
828 h( jr, j ) = -dconjg( s )*h( jr, j+1 ) + c*h( jr, j )
831 DO 130 jr = ifrstm, j
832 ctemp = c*t( jr, j+1 ) + s*t( jr, j )
833 t( jr, j ) = -dconjg( s )*t( jr, j+1 ) + c*t( jr, j )
838 ctemp = c*z( jr, j+1 ) + s*z( jr, j )
839 z( jr, j ) = -dconjg( s )*z( jr, j+1 ) + c*z( jr, j )
861 DO 200 j = 1, ilo - 1
862 absb = abs( t( j, j ) )
863 IF( absb.GT.safmin )
THEN
864 signbc = dconjg( t( j, j ) / absb )
867 CALL zscal( j-1, signbc, t( 1, j ), 1 )
868 CALL zscal( j, signbc, h( 1, j ), 1 )
870 CALL zscal( 1, signbc, h( j, j ), 1 )
873 $
CALL zscal( n, signbc, z( 1, j ), 1 )
877 alpha( j ) = h( j, j )
878 beta( j ) = t( j, j )
888 work( 1 ) = dcmplx( n )
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
subroutine zhgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, INFO)
ZHGEQZ
subroutine zrot(N, CX, INCX, CY, INCY, C, S)
ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.