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,
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 )
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 h(
j,
j ) = h(
j,
j )*signbc
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 h( ilast, ilast ) = h( ilast, ilast )*signbc
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 h(
j,
j ) = h(
j,
j )*signbc
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 zhgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, INFO)
ZHGEQZ
double precision function zlanhs(NORM, N, A, LDA, WORK)
ZLANHS returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine zlartg(F, G, CS, SN, R)
ZLARTG generates a plane rotation with real cosine and complex sine.
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
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...
double precision function dlamch(CMACH)
DLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
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