283 SUBROUTINE chgeqz( 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
298 COMPLEX alpha( * ), beta( * ), h( ldh, * ),
299 $ q( ldq, * ), t( ldt, * ), work( * ),
307 parameter( czero = ( 0.0e+0, 0.0e+0 ),
308 $ cone = ( 1.0e+0, 0.0e+0 ) )
310 parameter( zero = 0.0e+0, one = 1.0e+0 )
312 parameter( half = 0.5e+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 REAL absb, anorm, ascale, atol, bnorm, bscale, btol,
320 $ c, safmin, temp, temp2, tempr, ulp
321 COMPLEX abi22, ad11, ad12, ad21, ad22, ctemp, ctemp2,
322 $ ctemp3, eshift, rtdisc, s, shift, signbc, t1,
334 INTRINSIC abs, aimag, cmplx, conjg, max, min,
REAL, sqrt
340 abs1(
x ) = abs(
REAL( X ) ) + abs( aimag(
x ) )
346 IF(
lsame( job,
'E' ) )
THEN
349 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
369 IF(
lsame( compz,
'N' ) )
THEN
372 ELSE IF(
lsame( compz,
'V' ) )
THEN
375 ELSE IF(
lsame( compz,
'I' ) )
THEN
385 work( 1 ) = max( 1, n )
386 lquery = ( lwork.EQ.-1 )
387 IF( ischur.EQ.0 )
THEN
389 ELSE IF( icompq.EQ.0 )
THEN
391 ELSE IF( icompz.EQ.0 )
THEN
393 ELSE IF( n.LT.0 )
THEN
395 ELSE IF( ilo.LT.1 )
THEN
397 ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 )
THEN
399 ELSE IF( ldh.LT.n )
THEN
401 ELSE IF( ldt.LT.n )
THEN
403 ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) )
THEN
405 ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) )
THEN
407 ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery )
THEN
411 CALL
xerbla(
'CHGEQZ', -info )
413 ELSE IF( lquery )
THEN
421 work( 1 ) = cmplx( 1 )
428 $ CALL
claset(
'Full', n, n, czero, cone, q, ldq )
430 $ CALL
claset(
'Full', n, n, czero, cone, z, ldz )
437 anorm =
clanhs(
'F', in, h( ilo, ilo ), ldh, rwork )
438 bnorm =
clanhs(
'F', in, t( ilo, ilo ), ldt, rwork )
439 atol = max( safmin, ulp*anorm )
440 btol = max( safmin, ulp*bnorm )
441 ascale = one / max( safmin, anorm )
442 bscale = one / max( safmin, bnorm )
448 absb = abs( t(
j,
j ) )
449 IF( absb.GT.safmin )
THEN
450 signbc = conjg( t(
j,
j ) / absb )
453 CALL
cscal(
j-1, signbc, t( 1,
j ), 1 )
454 CALL
cscal(
j, signbc, h( 1,
j ), 1 )
456 h(
j,
j ) = h(
j,
j )*signbc
459 $ CALL
cscal( n, signbc, z( 1,
j ), 1 )
463 alpha(
j ) = h(
j,
j )
464 beta(
j ) = t(
j,
j )
497 maxit = 30*( ihi-ilo+1 )
499 DO 170 jiter = 1, maxit
514 IF( ilast.EQ.ilo )
THEN
517 IF( abs1( h( ilast, ilast-1 ) ).LE.atol )
THEN
518 h( ilast, ilast-1 ) = czero
523 IF( abs( t( ilast, ilast ) ).LE.btol )
THEN
524 t( ilast, ilast ) = czero
530 DO 40
j = ilast - 1, ilo, -1
537 IF( abs1( h(
j,
j-1 ) ).LE.atol )
THEN
547 IF( abs( t(
j,
j ) ).LT.btol )
THEN
553 IF( .NOT.ilazro )
THEN
554 IF( abs1( h(
j,
j-1 ) )*( ascale*abs1( h(
j+1,
555 $
j ) ) ).LE.abs1( h(
j,
j ) )*( ascale*atol ) )
565 IF( ilazro .OR. ilazr2 )
THEN
566 DO 20 jch =
j, ilast - 1
567 ctemp = h( jch, jch )
568 CALL
clartg( ctemp, h( jch+1, jch ), c, s,
570 h( jch+1, jch ) = czero
571 CALL
crot( ilastm-jch, h( jch, jch+1 ), ldh,
572 $ h( jch+1, jch+1 ), ldh, c, s )
573 CALL
crot( ilastm-jch, t( jch, jch+1 ), ldt,
574 $ t( jch+1, jch+1 ), ldt, c, s )
576 $ CALL
crot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
579 $ h( jch, jch-1 ) = h( jch, jch-1 )*c
581 IF( abs1( t( jch+1, jch+1 ) ).GE.btol )
THEN
582 IF( jch+1.GE.ilast )
THEN
589 t( jch+1, jch+1 ) = czero
597 DO 30 jch =
j, ilast - 1
598 ctemp = t( jch, jch+1 )
599 CALL
clartg( ctemp, t( jch+1, jch+1 ), c, s,
601 t( jch+1, jch+1 ) = czero
602 IF( jch.LT.ilastm-1 )
603 $ CALL
crot( ilastm-jch-1, t( jch, jch+2 ), ldt,
604 $ t( jch+1, jch+2 ), ldt, c, s )
605 CALL
crot( ilastm-jch+2, h( jch, jch-1 ), ldh,
606 $ h( jch+1, jch-1 ), ldh, c, s )
608 $ CALL
crot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
610 ctemp = h( jch+1, jch )
611 CALL
clartg( ctemp, h( jch+1, jch-1 ), c, s,
613 h( jch+1, jch-1 ) = czero
614 CALL
crot( jch+1-ifrstm, h( ifrstm, jch ), 1,
615 $ h( ifrstm, jch-1 ), 1, c, s )
616 CALL
crot( jch-ifrstm, t( ifrstm, jch ), 1,
617 $ t( ifrstm, jch-1 ), 1, c, s )
619 $ CALL
crot( n, z( 1, jch ), 1, z( 1, jch-1 ), 1,
624 ELSE IF( ilazro )
THEN
645 ctemp = h( ilast, ilast )
646 CALL
clartg( ctemp, h( ilast, ilast-1 ), c, s,
647 $ h( ilast, ilast ) )
648 h( ilast, ilast-1 ) = czero
649 CALL
crot( ilast-ifrstm, h( ifrstm, ilast ), 1,
650 $ h( ifrstm, ilast-1 ), 1, c, s )
651 CALL
crot( ilast-ifrstm, t( ifrstm, ilast ), 1,
652 $ t( ifrstm, ilast-1 ), 1, c, s )
654 $ CALL
crot( n, z( 1, ilast ), 1, z( 1, ilast-1 ), 1, c, s )
659 absb = abs( t( ilast, ilast ) )
660 IF( absb.GT.safmin )
THEN
661 signbc = conjg( t( ilast, ilast ) / absb )
662 t( ilast, ilast ) = absb
664 CALL
cscal( ilast-ifrstm, signbc, t( ifrstm, ilast ), 1 )
665 CALL
cscal( ilast+1-ifrstm, signbc, h( ifrstm, ilast ),
668 h( ilast, ilast ) = h( ilast, ilast )*signbc
671 $ CALL
cscal( n, signbc, z( 1, ilast ), 1 )
673 t( ilast, ilast ) = czero
675 alpha( ilast ) = h( ilast, ilast )
676 beta( ilast ) = t( ilast, ilast )
688 IF( .NOT.ilschr )
THEN
690 IF( ifrstm.GT.ilast )
702 IF( .NOT.ilschr )
THEN
712 IF( ( iiter / 10 )*10.NE.iiter )
THEN
721 u12 = ( bscale*t( ilast-1, ilast ) ) /
722 $ ( bscale*t( ilast, ilast ) )
723 ad11 = ( ascale*h( ilast-1, ilast-1 ) ) /
724 $ ( bscale*t( ilast-1, ilast-1 ) )
725 ad21 = ( ascale*h( ilast, ilast-1 ) ) /
726 $ ( bscale*t( ilast-1, ilast-1 ) )
727 ad12 = ( ascale*h( ilast-1, ilast ) ) /
728 $ ( bscale*t( ilast, ilast ) )
729 ad22 = ( ascale*h( ilast, ilast ) ) /
730 $ ( bscale*t( ilast, ilast ) )
731 abi22 = ad22 - u12*ad21
733 t1 = half*( ad11+abi22 )
734 rtdisc = sqrt( t1**2+ad12*ad21-ad11*ad22 )
735 temp =
REAL( t1-abi22 )*
REAL( RTDISC ) +
736 $ aimag( t1-abi22 )*aimag( rtdisc )
737 IF( temp.LE.zero )
THEN
746 eshift = eshift + (ascale*h(ilast,ilast-1))/
747 $ (bscale*t(ilast-1,ilast-1))
753 DO 80
j = ilast - 1, ifirst + 1, -1
755 ctemp = ascale*h(
j,
j ) - shift*( bscale*t(
j,
j ) )
757 temp2 = ascale*abs1( h(
j+1,
j ) )
758 tempr = max( temp, temp2 )
759 IF( tempr.LT.one .AND. tempr.NE.zero )
THEN
761 temp2 = temp2 / tempr
763 IF( abs1( h(
j,
j-1 ) )*temp2.LE.temp*atol )
768 ctemp = ascale*h( ifirst, ifirst ) -
769 $ shift*( bscale*t( ifirst, ifirst ) )
776 ctemp2 = ascale*h( istart+1, istart )
777 CALL
clartg( ctemp, ctemp2, c, s, ctemp3 )
781 DO 150
j = istart, ilast - 1
782 IF(
j.GT.istart )
THEN
784 CALL
clartg( ctemp, h(
j+1,
j-1 ), c, s, h(
j,
j-1 ) )
785 h(
j+1,
j-1 ) = czero
788 DO 100 jc =
j, ilastm
789 ctemp = c*h(
j, jc ) + s*h(
j+1, jc )
790 h(
j+1, jc ) = -conjg( s )*h(
j, jc ) + c*h(
j+1, jc )
792 ctemp2 = c*t(
j, jc ) + s*t(
j+1, jc )
793 t(
j+1, jc ) = -conjg( s )*t(
j, jc ) + c*t(
j+1, jc )
798 ctemp = c*q( jr,
j ) + conjg( s )*q( jr,
j+1 )
799 q( jr,
j+1 ) = -s*q( jr,
j ) + c*q( jr,
j+1 )
804 ctemp = t(
j+1,
j+1 )
805 CALL
clartg( ctemp, t(
j+1,
j ), c, s, t(
j+1,
j+1 ) )
808 DO 120 jr = ifrstm, min(
j+2, ilast )
809 ctemp = c*h( jr,
j+1 ) + s*h( jr,
j )
810 h( jr,
j ) = -conjg( s )*h( jr,
j+1 ) + c*h( jr,
j )
813 DO 130 jr = ifrstm,
j
814 ctemp = c*t( jr,
j+1 ) + s*t( jr,
j )
815 t( jr,
j ) = -conjg( s )*t( jr,
j+1 ) + c*t( jr,
j )
820 ctemp = c*z( jr,
j+1 ) + s*z( jr,
j )
821 z( jr,
j ) = -conjg( s )*z( jr,
j+1 ) + c*z( jr,
j )
843 DO 200
j = 1, ilo - 1
844 absb = abs( t(
j,
j ) )
845 IF( absb.GT.safmin )
THEN
846 signbc = conjg( t(
j,
j ) / absb )
849 CALL
cscal(
j-1, signbc, t( 1,
j ), 1 )
850 CALL
cscal(
j, signbc, h( 1,
j ), 1 )
852 h(
j,
j ) = h(
j,
j )*signbc
855 $ CALL
cscal( n, signbc, z( 1,
j ), 1 )
859 alpha(
j ) = h(
j,
j )
860 beta(
j ) = t(
j,
j )
870 work( 1 ) = cmplx( n )
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine crot(N, CX, INCX, CY, INCY, C, S)
CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors...
LOGICAL function lsame(CA, CB)
LSAME
subroutine cscal(N, CA, CX, INCX)
CSCAL
REAL function slamch(CMACH)
SLAMCH
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine clartg(F, G, CS, SN, R)
CLARTG generates a plane rotation with real cosine and complex sine.
REAL function clanhs(NORM, N, A, LDA, WORK)
CLANHS returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine chgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, INFO)
CHGEQZ