chgeqz.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 783 行 · 第 1/4 页

</span><span class="comment">*</span><span class="comment">                     i=INFO-N+1,...,N should be correct.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  We assume that complex ABS works as long as its value is less than
</span><span class="comment">*</span><span class="comment">  overflow.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      COMPLEX            CZERO, CONE
      PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
     $                   CONE = ( 1.0E+0, 0.0E+0 ) )
      REAL               ZERO, ONE
      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
      REAL               HALF
      PARAMETER          ( HALF = 0.5E+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            ILAZR2, ILAZRO, ILQ, ILSCHR, ILZ, LQUERY
      INTEGER            ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
     $                   ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER,
     $                   JR, MAXIT
      REAL               ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL,
     $                   C, SAFMIN, TEMP, TEMP2, TEMPR, ULP
      COMPLEX            ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2,
     $                   CTEMP3, ESHIFT, RTDISC, S, SHIFT, SIGNBC, T1,
     $                   U12, X
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.209"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      REAL               <a name="CLANHS.210"></a><a href="clanhs.f.html#CLANHS.1">CLANHS</a>, <a name="SLAMCH.210"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
      EXTERNAL           <a name="LSAME.211"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="CLANHS.211"></a><a href="clanhs.f.html#CLANHS.1">CLANHS</a>, <a name="SLAMCH.211"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="CLARTG.214"></a><a href="clartg.f.html#CLARTG.1">CLARTG</a>, <a name="CLASET.214"></a><a href="claset.f.html#CLASET.1">CLASET</a>, <a name="CROT.214"></a><a href="crot.f.html#CROT.1">CROT</a>, CSCAL, <a name="XERBLA.214"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL, SQRT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Functions ..
</span>      REAL               ABS1
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Function definitions ..
</span>      ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Decode JOB, COMPQ, COMPZ
</span><span class="comment">*</span><span class="comment">
</span>      IF( <a name="LSAME.229"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'E'</span> ) ) THEN
         ILSCHR = .FALSE.
         ISCHUR = 1
      ELSE IF( <a name="LSAME.232"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'S'</span> ) ) THEN
         ILSCHR = .TRUE.
         ISCHUR = 2
      ELSE
         ISCHUR = 0
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( <a name="LSAME.239"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPQ, <span class="string">'N'</span> ) ) THEN
         ILQ = .FALSE.
         ICOMPQ = 1
      ELSE IF( <a name="LSAME.242"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPQ, <span class="string">'V'</span> ) ) THEN
         ILQ = .TRUE.
         ICOMPQ = 2
      ELSE IF( <a name="LSAME.245"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPQ, <span class="string">'I'</span> ) ) THEN
         ILQ = .TRUE.
         ICOMPQ = 3
      ELSE
         ICOMPQ = 0
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( <a name="LSAME.252"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPZ, <span class="string">'N'</span> ) ) THEN
         ILZ = .FALSE.
         ICOMPZ = 1
      ELSE IF( <a name="LSAME.255"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPZ, <span class="string">'V'</span> ) ) THEN
         ILZ = .TRUE.
         ICOMPZ = 2
      ELSE IF( <a name="LSAME.258"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPZ, <span class="string">'I'</span> ) ) THEN
         ILZ = .TRUE.
         ICOMPZ = 3
      ELSE
         ICOMPZ = 0
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Check Argument Values
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      WORK( 1 ) = MAX( 1, N )
      LQUERY = ( LWORK.EQ.-1 )
      IF( ISCHUR.EQ.0 ) THEN
         INFO = -1
      ELSE IF( ICOMPQ.EQ.0 ) THEN
         INFO = -2
      ELSE IF( ICOMPZ.EQ.0 ) THEN
         INFO = -3
      ELSE IF( N.LT.0 ) THEN
         INFO = -4
      ELSE IF( ILO.LT.1 ) THEN
         INFO = -5
      ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
         INFO = -6
      ELSE IF( LDH.LT.N ) THEN
         INFO = -8
      ELSE IF( LDT.LT.N ) THEN
         INFO = -10
      ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN
         INFO = -14
      ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN
         INFO = -16
      ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
         INFO = -18
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.294"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CHGEQZ.294"></a><a href="chgeqz.f.html#CHGEQZ.1">CHGEQZ</a>'</span>, -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">c</span><span class="comment">     WORK( 1 ) = CMPLX( 1 )
</span>      IF( N.LE.0 ) THEN
         WORK( 1 ) = CMPLX( 1 )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initialize Q and Z
</span><span class="comment">*</span><span class="comment">
</span>      IF( ICOMPQ.EQ.3 )
     $   CALL <a name="CLASET.311"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'Full'</span>, N, N, CZERO, CONE, Q, LDQ )
      IF( ICOMPZ.EQ.3 )
     $   CALL <a name="CLASET.313"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'Full'</span>, N, N, CZERO, CONE, Z, LDZ )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Machine Constants
</span><span class="comment">*</span><span class="comment">
</span>      IN = IHI + 1 - ILO
      SAFMIN = <a name="SLAMCH.318"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'S'</span> )
      ULP = <a name="SLAMCH.319"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'E'</span> )*<a name="SLAMCH.319"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'B'</span> )
      ANORM = <a name="CLANHS.320"></a><a href="clanhs.f.html#CLANHS.1">CLANHS</a>( <span class="string">'F'</span>, IN, H( ILO, ILO ), LDH, RWORK )
      BNORM = <a name="CLANHS.321"></a><a href="clanhs.f.html#CLANHS.1">CLANHS</a>( <span class="string">'F'</span>, IN, T( ILO, ILO ), LDT, RWORK )
      ATOL = MAX( SAFMIN, ULP*ANORM )
      BTOL = MAX( SAFMIN, ULP*BNORM )
      ASCALE = ONE / MAX( SAFMIN, ANORM )
      BSCALE = ONE / MAX( SAFMIN, BNORM )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set Eigenvalues IHI+1:N
</span><span class="comment">*</span><span class="comment">
</span>      DO 10 J = IHI + 1, N
         ABSB = ABS( T( J, J ) )
         IF( ABSB.GT.SAFMIN ) THEN
            SIGNBC = CONJG( T( J, J ) / ABSB )
            T( J, J ) = ABSB
            IF( ILSCHR ) THEN
               CALL CSCAL( J-1, SIGNBC, T( 1, J ), 1 )
               CALL CSCAL( J, SIGNBC, H( 1, J ), 1 )
            ELSE
               H( J, J ) = H( J, J )*SIGNBC
            END IF
            IF( ILZ )
     $         CALL CSCAL( N, SIGNBC, Z( 1, J ), 1 )
         ELSE
            T( J, J ) = CZERO
         END IF
         ALPHA( J ) = H( J, J )
         BETA( J ) = T( J, J )
   10 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If IHI &lt; ILO, skip QZ steps
</span><span class="comment">*</span><span class="comment">
</span>      IF( IHI.LT.ILO )
     $   GO TO 190
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     MAIN QZ ITERATION LOOP
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initialize dynamic indices
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Eigenvalues ILAST+1:N have been found.
</span><span class="comment">*</span><span class="comment">        Column operations modify rows IFRSTM:whatever
</span><span class="comment">*</span><span class="comment">        Row operations modify columns whatever:ILASTM
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If only eigenvalues are being computed, then
</span><span class="comment">*</span><span class="comment">        IFRSTM is the row of the last splitting row above row ILAST;
</span><span class="comment">*</span><span class="comment">        this is always at least ILO.
</span><span class="comment">*</span><span class="comment">     IITER counts iterations since the last eigenvalue was found,
</span><span class="comment">*</span><span class="comment">        to tell when to use an extraordinary shift.
</span><span class="comment">*</span><span class="comment">     MAXIT is the maximum number of QZ sweeps allowed.
</span><span class="comment">*</span><span class="comment">
</span>      ILAST = IHI
      IF( ILSCHR ) THEN
         IFRSTM = 1
         ILASTM = N