slahqr.f.html

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

</span>            H21S = H( M+1, M )
            S = ABS( H( M, M )-RT2R ) + ABS( RT2I ) + ABS( H21S )
            H21S = H( M+1, M ) / S
            V( 1 ) = H21S*H( M, M+1 ) + ( H( M, M )-RT1R )*
     $               ( ( H( M, M )-RT2R ) / S ) - RT1I*( RT2I / S )
            V( 2 ) = H21S*( H( M, M )+H( M+1, M+1 )-RT1R-RT2R )
            V( 3 ) = H21S*H( M+2, M+1 )
            S = ABS( V( 1 ) ) + ABS( V( 2 ) ) + ABS( V( 3 ) )
            V( 1 ) = V( 1 ) / S
            V( 2 ) = V( 2 ) / S
            V( 3 ) = V( 3 ) / S
            IF( M.EQ.L )
     $         GO TO 60
            IF( ABS( H( M, M-1 ) )*( ABS( V( 2 ) )+ABS( V( 3 ) ) ).LE.
     $          ULP*ABS( V( 1 ) )*( ABS( H( M-1, M-1 ) )+ABS( H( M,
     $          M ) )+ABS( H( M+1, M+1 ) ) ) )GO TO 60
   50    CONTINUE
   60    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Double-shift QR step
</span><span class="comment">*</span><span class="comment">
</span>         DO 130 K = M, I - 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           The first iteration of this loop determines a reflection G
</span><span class="comment">*</span><span class="comment">           from the vector V and applies it from left and right to H,
</span><span class="comment">*</span><span class="comment">           thus creating a nonzero bulge below the subdiagonal.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Each subsequent iteration determines a reflection G to
</span><span class="comment">*</span><span class="comment">           restore the Hessenberg form in the (K-1)th column, and thus
</span><span class="comment">*</span><span class="comment">           chases the bulge one step toward the bottom of the active
</span><span class="comment">*</span><span class="comment">           submatrix. NR is the order of G.
</span><span class="comment">*</span><span class="comment">
</span>            NR = MIN( 3, I-K+1 )
            IF( K.GT.M )
     $         CALL SCOPY( NR, H( K, K-1 ), 1, V, 1 )
            CALL <a name="SLARFG.369"></a><a href="slarfg.f.html#SLARFG.1">SLARFG</a>( NR, V( 1 ), V( 2 ), 1, T1 )
            IF( K.GT.M ) THEN
               H( K, K-1 ) = V( 1 )
               H( K+1, K-1 ) = ZERO
               IF( K.LT.I-1 )
     $            H( K+2, K-1 ) = ZERO
            ELSE IF( M.GT.L ) THEN
               H( K, K-1 ) = -H( K, K-1 )
            END IF
            V2 = V( 2 )
            T2 = T1*V2
            IF( NR.EQ.3 ) THEN
               V3 = V( 3 )
               T3 = T1*V3
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Apply G from the left to transform the rows of the matrix
</span><span class="comment">*</span><span class="comment">              in columns K to I2.
</span><span class="comment">*</span><span class="comment">
</span>               DO 70 J = K, I2
                  SUM = H( K, J ) + V2*H( K+1, J ) + V3*H( K+2, J )
                  H( K, J ) = H( K, J ) - SUM*T1
                  H( K+1, J ) = H( K+1, J ) - SUM*T2
                  H( K+2, J ) = H( K+2, J ) - SUM*T3
   70          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Apply G from the right to transform the columns of the
</span><span class="comment">*</span><span class="comment">              matrix in rows I1 to min(K+3,I).
</span><span class="comment">*</span><span class="comment">
</span>               DO 80 J = I1, MIN( K+3, I )
                  SUM = H( J, K ) + V2*H( J, K+1 ) + V3*H( J, K+2 )
                  H( J, K ) = H( J, K ) - SUM*T1
                  H( J, K+1 ) = H( J, K+1 ) - SUM*T2
                  H( J, K+2 ) = H( J, K+2 ) - SUM*T3
   80          CONTINUE
<span class="comment">*</span><span class="comment">
</span>               IF( WANTZ ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Accumulate transformations in the matrix Z
</span><span class="comment">*</span><span class="comment">
</span>                  DO 90 J = ILOZ, IHIZ
                     SUM = Z( J, K ) + V2*Z( J, K+1 ) + V3*Z( J, K+2 )
                     Z( J, K ) = Z( J, K ) - SUM*T1
                     Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2
                     Z( J, K+2 ) = Z( J, K+2 ) - SUM*T3
   90             CONTINUE
               END IF
            ELSE IF( NR.EQ.2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Apply G from the left to transform the rows of the matrix
</span><span class="comment">*</span><span class="comment">              in columns K to I2.
</span><span class="comment">*</span><span class="comment">
</span>               DO 100 J = K, I2
                  SUM = H( K, J ) + V2*H( K+1, J )
                  H( K, J ) = H( K, J ) - SUM*T1
                  H( K+1, J ) = H( K+1, J ) - SUM*T2
  100          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Apply G from the right to transform the columns of the
</span><span class="comment">*</span><span class="comment">              matrix in rows I1 to min(K+3,I).
</span><span class="comment">*</span><span class="comment">
</span>               DO 110 J = I1, I
                  SUM = H( J, K ) + V2*H( J, K+1 )
                  H( J, K ) = H( J, K ) - SUM*T1
                  H( J, K+1 ) = H( J, K+1 ) - SUM*T2
  110          CONTINUE
<span class="comment">*</span><span class="comment">
</span>               IF( WANTZ ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Accumulate transformations in the matrix Z
</span><span class="comment">*</span><span class="comment">
</span>                  DO 120 J = ILOZ, IHIZ
                     SUM = Z( J, K ) + V2*Z( J, K+1 )
                     Z( J, K ) = Z( J, K ) - SUM*T1
                     Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2
  120             CONTINUE
               END IF
            END IF
  130    CONTINUE
<span class="comment">*</span><span class="comment">
</span>  140 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Failure to converge in remaining number of iterations
</span><span class="comment">*</span><span class="comment">
</span>      INFO = I
      RETURN
<span class="comment">*</span><span class="comment">
</span>  150 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      IF( L.EQ.I ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        H(I,I-1) is negligible: one eigenvalue has converged.
</span><span class="comment">*</span><span class="comment">
</span>         WR( I ) = H( I, I )
         WI( I ) = ZERO
      ELSE IF( L.EQ.I-1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Transform the 2-by-2 submatrix to standard Schur form,
</span><span class="comment">*</span><span class="comment">        and compute and store the eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SLANV2.470"></a><a href="slanv2.f.html#SLANV2.1">SLANV2</a>( H( I-1, I-1 ), H( I-1, I ), H( I, I-1 ),
     $                H( I, I ), WR( I-1 ), WI( I-1 ), WR( I ), WI( I ),
     $                CS, SN )
<span class="comment">*</span><span class="comment">
</span>         IF( WANTT ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Apply the transformation to the rest of H.
</span><span class="comment">*</span><span class="comment">
</span>            IF( I2.GT.I )
     $         CALL SROT( I2-I, H( I-1, I+1 ), LDH, H( I, I+1 ), LDH,
     $                    CS, SN )
            CALL SROT( I-I1-1, H( I1, I-1 ), 1, H( I1, I ), 1, CS, SN )
         END IF
         IF( WANTZ ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Apply the transformation to Z.
</span><span class="comment">*</span><span class="comment">
</span>            CALL SROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, SN )
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     return to start of the main loop with new value of I.
</span><span class="comment">*</span><span class="comment">
</span>      I = L - 1
      GO TO 20
<span class="comment">*</span><span class="comment">
</span>  160 CONTINUE
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="SLAHQR.499"></a><a href="slahqr.f.html#SLAHQR.1">SLAHQR</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>