283 SUBROUTINE zhgeqz( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
284 $ alpha, beta, q, ldq, z, ldz, work, lwork,
293 CHARACTER COMPQ, COMPZ, JOB
294 INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
297 DOUBLE PRECISION RWORK( * )
298 COMPLEX*16 ALPHA( * ), BETA( * ), H( ldh, * ),
299 $ q( ldq, * ), t( ldt, * ), work( * ),
306 COMPLEX*16 CZERO, CONE
307 parameter ( czero = ( 0.0d+0, 0.0d+0 ),
308 $ cone = ( 1.0d+0, 0.0d+0 ) )
309 DOUBLE PRECISION ZERO, ONE
310 parameter ( zero = 0.0d+0, one = 1.0d+0 )
311 DOUBLE PRECISION HALF
312 parameter ( half = 0.5d+0 )
315 LOGICAL ILAZR2, ILAZRO, ILQ, ILSCHR, ILZ, LQUERY
316 INTEGER ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
317 $ ilastm, in, ischur, istart, j, jc, jch, jiter,
319 DOUBLE PRECISION ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL,
320 $ c, safmin, temp, temp2, tempr, ulp
321 COMPLEX*16 ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2,
322 $ ctemp3, eshift, rtdisc, s, shift, signbc, t1,
327 DOUBLE PRECISION DLAMCH, ZLANHS
328 EXTERNAL lsame, dlamch, zlanhs
334 INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, min,
338 DOUBLE PRECISION ABS1
341 abs1( x ) = abs( dble( x ) ) + abs( dimag( x ) )
347 IF( lsame( job,
'E' ) )
THEN
350 ELSE IF( lsame( job,
'S' ) )
THEN
357 IF( lsame( compq,
'N' ) )
THEN
360 ELSE IF( lsame( compq,
'V' ) )
THEN
363 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
386 work( 1 ) = max( 1, n )
387 lquery = ( lwork.EQ.-1 )
388 IF( ischur.EQ.0 )
THEN
390 ELSE IF( icompq.EQ.0 )
THEN
392 ELSE IF( icompz.EQ.0 )
THEN
394 ELSE IF( n.LT.0 )
THEN
396 ELSE IF( ilo.LT.1 )
THEN
398 ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 )
THEN
400 ELSE IF( ldh.LT.n )
THEN
402 ELSE IF( ldt.LT.n )
THEN
404 ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) )
THEN
406 ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) )
THEN
408 ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery )
THEN
412 CALL xerbla(
'ZHGEQZ', -info )
414 ELSE IF( lquery )
THEN
422 work( 1 ) = dcmplx( 1 )
429 $
CALL zlaset(
'Full', n, n, czero, cone, q, ldq )
431 $
CALL zlaset(
'Full', n, n, czero, cone, z, ldz )
436 safmin = dlamch(
'S' )
437 ulp = dlamch(
'E' )*dlamch(
'B' )
438 anorm = zlanhs(
'F', in, h( ilo, ilo ), ldh, rwork )
439 bnorm = zlanhs(
'F', in, t( ilo, ilo ), ldt, rwork )
440 atol = max( safmin, ulp*anorm )
441 btol = max( safmin, ulp*bnorm )
442 ascale = one / max( safmin, anorm )
443 bscale = one / max( safmin, bnorm )
449 absb = abs( t( j, j ) )
450 IF( absb.GT.safmin )
THEN
451 signbc = dconjg( t( j, j ) / absb )
454 CALL zscal( j-1, signbc, t( 1, j ), 1 )
455 CALL zscal( j, signbc, h( 1, j ), 1 )
457 CALL zscal( 1, signbc, h( j, j ), 1 )
460 $
CALL zscal( n, signbc, z( 1, j ), 1 )
464 alpha( j ) = h( j, j )
465 beta( j ) = t( j, j )
498 maxit = 30*( ihi-ilo+1 )
500 DO 170 jiter = 1, maxit
515 IF( ilast.EQ.ilo )
THEN
518 IF( abs1( h( ilast, ilast-1 ) ).LE.atol )
THEN
519 h( ilast, ilast-1 ) = czero
524 IF( abs( t( ilast, ilast ) ).LE.btol )
THEN
525 t( ilast, ilast ) = czero
531 DO 40 j = ilast - 1, ilo, -1
538 IF( abs1( h( j, j-1 ) ).LE.atol )
THEN
548 IF( abs( t( j, j ) ).LT.btol )
THEN
554 IF( .NOT.ilazro )
THEN
555 IF( abs1( h( j, j-1 ) )*( ascale*abs1( h( j+1,
556 $ j ) ) ).LE.abs1( h( j, j ) )*( ascale*atol ) )
566 IF( ilazro .OR. ilazr2 )
THEN
567 DO 20 jch = j, ilast - 1
568 ctemp = h( jch, jch )
569 CALL zlartg( ctemp, h( jch+1, jch ), c, s,
571 h( jch+1, jch ) = czero
572 CALL zrot( ilastm-jch, h( jch, jch+1 ), ldh,
573 $ h( jch+1, jch+1 ), ldh, c, s )
574 CALL zrot( ilastm-jch, t( jch, jch+1 ), ldt,
575 $ t( jch+1, jch+1 ), ldt, c, s )
577 $
CALL zrot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
580 $ h( jch, jch-1 ) = h( jch, jch-1 )*c
582 IF( abs1( t( jch+1, jch+1 ) ).GE.btol )
THEN
583 IF( jch+1.GE.ilast )
THEN
590 t( jch+1, jch+1 ) = czero
598 DO 30 jch = j, ilast - 1
599 ctemp = t( jch, jch+1 )
600 CALL zlartg( ctemp, t( jch+1, jch+1 ), c, s,
602 t( jch+1, jch+1 ) = czero
603 IF( jch.LT.ilastm-1 )
604 $
CALL zrot( ilastm-jch-1, t( jch, jch+2 ), ldt,
605 $ t( jch+1, jch+2 ), ldt, c, s )
606 CALL zrot( ilastm-jch+2, h( jch, jch-1 ), ldh,
607 $ h( jch+1, jch-1 ), ldh, c, s )
609 $
CALL zrot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
611 ctemp = h( jch+1, jch )
612 CALL zlartg( ctemp, h( jch+1, jch-1 ), c, s,
614 h( jch+1, jch-1 ) = czero
615 CALL zrot( jch+1-ifrstm, h( ifrstm, jch ), 1,
616 $ h( ifrstm, jch-1 ), 1, c, s )
617 CALL zrot( jch-ifrstm, t( ifrstm, jch ), 1,
618 $ t( ifrstm, jch-1 ), 1, c, s )
620 $
CALL zrot( n, z( 1, jch ), 1, z( 1, jch-1 ), 1,
625 ELSE IF( ilazro )
THEN
646 ctemp = h( ilast, ilast )
647 CALL zlartg( ctemp, h( ilast, ilast-1 ), c, s,
648 $ h( ilast, ilast ) )
649 h( ilast, ilast-1 ) = czero
650 CALL zrot( ilast-ifrstm, h( ifrstm, ilast ), 1,
651 $ h( ifrstm, ilast-1 ), 1, c, s )
652 CALL zrot( ilast-ifrstm, t( ifrstm, ilast ), 1,
653 $ t( ifrstm, ilast-1 ), 1, c, s )
655 $
CALL zrot( n, z( 1, ilast ), 1, z( 1, ilast-1 ), 1, c, s )
660 absb = abs( t( ilast, ilast ) )
661 IF( absb.GT.safmin )
THEN
662 signbc = dconjg( t( ilast, ilast ) / absb )
663 t( ilast, ilast ) = absb
665 CALL zscal( ilast-ifrstm, signbc, t( ifrstm, ilast ), 1 )
666 CALL zscal( ilast+1-ifrstm, signbc, h( ifrstm, ilast ),
669 CALL zscal( 1, signbc, h( ilast, ilast ), 1 )
672 $
CALL zscal( n, signbc, z( 1, ilast ), 1 )
674 t( ilast, ilast ) = czero
676 alpha( ilast ) = h( ilast, ilast )
677 beta( ilast ) = t( ilast, ilast )
689 IF( .NOT.ilschr )
THEN
691 IF( ifrstm.GT.ilast )
703 IF( .NOT.ilschr )
THEN
713 IF( ( iiter / 10 )*10.NE.iiter )
THEN
722 u12 = ( bscale*t( ilast-1, ilast ) ) /
723 $ ( bscale*t( ilast, ilast ) )
724 ad11 = ( ascale*h( ilast-1, ilast-1 ) ) /
725 $ ( bscale*t( ilast-1, ilast-1 ) )
726 ad21 = ( ascale*h( ilast, ilast-1 ) ) /
727 $ ( bscale*t( ilast-1, ilast-1 ) )
728 ad12 = ( ascale*h( ilast-1, ilast ) ) /
729 $ ( bscale*t( ilast, ilast ) )
730 ad22 = ( ascale*h( ilast, ilast ) ) /
731 $ ( bscale*t( ilast, ilast ) )
732 abi22 = ad22 - u12*ad21
734 t1 = half*( ad11+abi22 )
735 rtdisc = sqrt( t1**2+ad12*ad21-ad11*ad22 )
736 temp = dble( t1-abi22 )*dble( rtdisc ) +
737 $ dimag( t1-abi22 )*dimag( rtdisc )
738 IF( temp.LE.zero )
THEN
747 eshift = eshift + (ascale*h(ilast,ilast-1))/
748 $ (bscale*t(ilast-1,ilast-1))
754 DO 80 j = ilast - 1, ifirst + 1, -1
756 ctemp = ascale*h( j, j ) - shift*( bscale*t( j, j ) )
758 temp2 = ascale*abs1( h( j+1, j ) )
759 tempr = max( temp, temp2 )
760 IF( tempr.LT.one .AND. tempr.NE.zero )
THEN
762 temp2 = temp2 / tempr
764 IF( abs1( h( j, j-1 ) )*temp2.LE.temp*atol )
769 ctemp = ascale*h( ifirst, ifirst ) -
770 $ shift*( bscale*t( ifirst, ifirst ) )
777 ctemp2 = ascale*h( istart+1, istart )
778 CALL zlartg( ctemp, ctemp2, c, s, ctemp3 )
782 DO 150 j = istart, ilast - 1
783 IF( j.GT.istart )
THEN
785 CALL zlartg( ctemp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
786 h( j+1, j-1 ) = czero
789 DO 100 jc = j, ilastm
790 ctemp = c*h( j, jc ) + s*h( j+1, jc )
791 h( j+1, jc ) = -dconjg( s )*h( j, jc ) + c*h( j+1, jc )
793 ctemp2 = c*t( j, jc ) + s*t( j+1, jc )
794 t( j+1, jc ) = -dconjg( s )*t( j, jc ) + c*t( j+1, jc )
799 ctemp = c*q( jr, j ) + dconjg( s )*q( jr, j+1 )
800 q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
805 ctemp = t( j+1, j+1 )
806 CALL zlartg( ctemp, t( j+1, j ), c, s, t( j+1, j+1 ) )
809 DO 120 jr = ifrstm, min( j+2, ilast )
810 ctemp = c*h( jr, j+1 ) + s*h( jr, j )
811 h( jr, j ) = -dconjg( s )*h( jr, j+1 ) + c*h( jr, j )
814 DO 130 jr = ifrstm, j
815 ctemp = c*t( jr, j+1 ) + s*t( jr, j )
816 t( jr, j ) = -dconjg( s )*t( jr, j+1 ) + c*t( jr, j )
821 ctemp = c*z( jr, j+1 ) + s*z( jr, j )
822 z( jr, j ) = -dconjg( s )*z( jr, j+1 ) + c*z( jr, j )
844 DO 200 j = 1, ilo - 1
845 absb = abs( t( j, j ) )
846 IF( absb.GT.safmin )
THEN
847 signbc = dconjg( t( j, j ) / absb )
850 CALL zscal( j-1, signbc, t( 1, j ), 1 )
851 CALL zscal( j, signbc, h( 1, j ), 1 )
853 CALL zscal( 1, signbc, h( j, j ), 1 )
856 $
CALL zscal( n, signbc, z( 1, j ), 1 )
860 alpha( j ) = h( j, j )
861 beta( j ) = t( j, j )
871 work( 1 ) = dcmplx( n )
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...
subroutine xerbla(SRNAME, INFO)
XERBLA
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 zscal(N, ZA, ZX, INCX)
ZSCAL
subroutine zlartg(F, G, CS, SN, R)
ZLARTG generates a plane rotation with real cosine and complex sine.