chgeqz.f.html

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

      ELSE
         IFRSTM = ILO
         ILASTM = IHI
      END IF
      IITER = 0
      ESHIFT = CZERO
      MAXIT = 30*( IHI-ILO+1 )
<span class="comment">*</span><span class="comment">
</span>      DO 170 JITER = 1, MAXIT
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Check for too many iterations.
</span><span class="comment">*</span><span class="comment">
</span>         IF( JITER.GT.MAXIT )
     $      GO TO 180
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Split the matrix if possible.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Two tests:
</span><span class="comment">*</span><span class="comment">           1: H(j,j-1)=0  or  j=ILO
</span><span class="comment">*</span><span class="comment">           2: T(j,j)=0
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Special case: j=ILAST
</span><span class="comment">*</span><span class="comment">
</span>         IF( ILAST.EQ.ILO ) THEN
            GO TO 60
         ELSE
            IF( ABS1( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN
               H( ILAST, ILAST-1 ) = CZERO
               GO TO 60
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN
            T( ILAST, ILAST ) = CZERO
            GO TO 50
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        General case: j&lt;ILAST
</span><span class="comment">*</span><span class="comment">
</span>         DO 40 J = ILAST - 1, ILO, -1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Test 1: for H(j,j-1)=0 or j=ILO
</span><span class="comment">*</span><span class="comment">
</span>            IF( J.EQ.ILO ) THEN
               ILAZRO = .TRUE.
            ELSE
               IF( ABS1( H( J, J-1 ) ).LE.ATOL ) THEN
                  H( J, J-1 ) = CZERO
                  ILAZRO = .TRUE.
               ELSE
                  ILAZRO = .FALSE.
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Test 2: for T(j,j)=0
</span><span class="comment">*</span><span class="comment">
</span>            IF( ABS( T( J, J ) ).LT.BTOL ) THEN
               T( J, J ) = CZERO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Test 1a: Check for 2 consecutive small subdiagonals in A
</span><span class="comment">*</span><span class="comment">
</span>               ILAZR2 = .FALSE.
               IF( .NOT.ILAZRO ) THEN
                  IF( ABS1( H( J, J-1 ) )*( ASCALE*ABS1( H( J+1,
     $                J ) ) ).LE.ABS1( H( J, J ) )*( ASCALE*ATOL ) )
     $                ILAZR2 = .TRUE.
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              If both tests pass (1 &amp; 2), i.e., the leading diagonal
</span><span class="comment">*</span><span class="comment">              element of B in the block is zero, split a 1x1 block off
</span><span class="comment">*</span><span class="comment">              at the top. (I.e., at the J-th row/column) The leading
</span><span class="comment">*</span><span class="comment">              diagonal element of the remainder can also be zero, so
</span><span class="comment">*</span><span class="comment">              this may have to be done repeatedly.
</span><span class="comment">*</span><span class="comment">
</span>               IF( ILAZRO .OR. ILAZR2 ) THEN
                  DO 20 JCH = J, ILAST - 1
                     CTEMP = H( JCH, JCH )
                     CALL <a name="CLARTG.451"></a><a href="clartg.f.html#CLARTG.1">CLARTG</a>( CTEMP, H( JCH+1, JCH ), C, S,
     $                            H( JCH, JCH ) )
                     H( JCH+1, JCH ) = CZERO
                     CALL <a name="CROT.454"></a><a href="crot.f.html#CROT.1">CROT</a>( ILASTM-JCH, H( JCH, JCH+1 ), LDH,
     $                          H( JCH+1, JCH+1 ), LDH, C, S )
                     CALL <a name="CROT.456"></a><a href="crot.f.html#CROT.1">CROT</a>( ILASTM-JCH, T( JCH, JCH+1 ), LDT,
     $                          T( JCH+1, JCH+1 ), LDT, C, S )
                     IF( ILQ )
     $                  CALL <a name="CROT.459"></a><a href="crot.f.html#CROT.1">CROT</a>( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
     $                             C, CONJG( S ) )
                     IF( ILAZR2 )
     $                  H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C
                     ILAZR2 = .FALSE.
                     IF( ABS1( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN
                        IF( JCH+1.GE.ILAST ) THEN
                           GO TO 60
                        ELSE
                           IFIRST = JCH + 1
                           GO TO 70
                        END IF
                     END IF
                     T( JCH+1, JCH+1 ) = CZERO
   20             CONTINUE
                  GO TO 50
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Only test 2 passed -- chase the zero to T(ILAST,ILAST)
</span><span class="comment">*</span><span class="comment">                 Then process as in the case T(ILAST,ILAST)=0
</span><span class="comment">*</span><span class="comment">
</span>                  DO 30 JCH = J, ILAST - 1
                     CTEMP = T( JCH, JCH+1 )
                     CALL <a name="CLARTG.482"></a><a href="clartg.f.html#CLARTG.1">CLARTG</a>( CTEMP, T( JCH+1, JCH+1 ), C, S,
     $                            T( JCH, JCH+1 ) )
                     T( JCH+1, JCH+1 ) = CZERO
                     IF( JCH.LT.ILASTM-1 )
     $                  CALL <a name="CROT.486"></a><a href="crot.f.html#CROT.1">CROT</a>( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT,
     $                             T( JCH+1, JCH+2 ), LDT, C, S )
                     CALL <a name="CROT.488"></a><a href="crot.f.html#CROT.1">CROT</a>( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH,
     $                          H( JCH+1, JCH-1 ), LDH, C, S )
                     IF( ILQ )
     $                  CALL <a name="CROT.491"></a><a href="crot.f.html#CROT.1">CROT</a>( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
     $                             C, CONJG( S ) )
                     CTEMP = H( JCH+1, JCH )
                     CALL <a name="CLARTG.494"></a><a href="clartg.f.html#CLARTG.1">CLARTG</a>( CTEMP, H( JCH+1, JCH-1 ), C, S,
     $                            H( JCH+1, JCH ) )
                     H( JCH+1, JCH-1 ) = CZERO
                     CALL <a name="CROT.497"></a><a href="crot.f.html#CROT.1">CROT</a>( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1,
     $                          H( IFRSTM, JCH-1 ), 1, C, S )
                     CALL <a name="CROT.499"></a><a href="crot.f.html#CROT.1">CROT</a>( JCH-IFRSTM, T( IFRSTM, JCH ), 1,
     $                          T( IFRSTM, JCH-1 ), 1, C, S )
                     IF( ILZ )
     $                  CALL <a name="CROT.502"></a><a href="crot.f.html#CROT.1">CROT</a>( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1,
     $                             C, S )
   30             CONTINUE
                  GO TO 50
               END IF
            ELSE IF( ILAZRO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Only test 1 passed -- work on J:ILAST
</span><span class="comment">*</span><span class="comment">
</span>               IFIRST = J
               GO TO 70
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Neither test passed -- try next J
</span><span class="comment">*</span><span class="comment">
</span>   40    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        (Drop-through is &quot;impossible&quot;)
</span><span class="comment">*</span><span class="comment">
</span>         INFO = 2*N + 1
         GO TO 210
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
</span><span class="comment">*</span><span class="comment">        1x1 block.
</span><span class="comment">*</span><span class="comment">
</span>   50    CONTINUE
         CTEMP = H( ILAST, ILAST )
         CALL <a name="CLARTG.529"></a><a href="clartg.f.html#CLARTG.1">CLARTG</a>( CTEMP, H( ILAST, ILAST-1 ), C, S,
     $                H( ILAST, ILAST ) )
         H( ILAST, ILAST-1 ) = CZERO
         CALL <a name="CROT.532"></a><a href="crot.f.html#CROT.1">CROT</a>( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1,
     $              H( IFRSTM, ILAST-1 ), 1, C, S )
         CALL <a name="CROT.534"></a><a href="crot.f.html#CROT.1">CROT</a>( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1,
     $              T( IFRSTM, ILAST-1 ), 1, C, S )
         IF( ILZ )
     $      CALL <a name="CROT.537"></a><a href="crot.f.html#CROT.1">CROT</a>( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHA and BETA
</span><span class="comment">*</span><span class="comment">
</span>   60    CONTINUE
         ABSB = ABS( T( ILAST, ILAST ) )
         IF( ABSB.GT.SAFMIN ) THEN
            SIGNBC = CONJG( T( ILAST, ILAST ) / ABSB )
            T( ILAST, ILAST ) = ABSB
            IF( ILSCHR ) THEN
               CALL CSCAL( ILAST-IFRSTM, SIGNBC, T( IFRSTM, ILAST ), 1 )
               CALL CSCAL( ILAST+1-IFRSTM, SIGNBC, H( IFRSTM, ILAST ),
     $                     1 )
            ELSE
               H( ILAST, ILAST ) = H( ILAST, ILAST )*SIGNBC
            END IF
            IF( ILZ )
     $         CALL CSCAL( N, SIGNBC, Z( 1, ILAST ), 1 )
         ELSE
            T( ILAST, ILAST ) = CZERO
         END IF
         ALPHA( ILAST ) = H( ILAST, ILAST )
         BETA( ILAST ) = T( ILAST, ILAST )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Go to next block -- exit if finished.
</span><span class="comment">*</span><span class="comment">
</span>         ILAST = ILAST - 1
         IF( ILAST.LT.ILO )
     $      GO TO 190
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Reset counters
</span><span class="comment">*</span><span class="comment">
</span>         IITER = 0