ssteqr.f.html

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

</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Look for small superdiagonal element.
</span><span class="comment">*</span><span class="comment">
</span>   90    CONTINUE
         IF( L.NE.LEND ) THEN
            LENDP1 = LEND + 1
            DO 100 M = L, LENDP1, -1
               TST = ABS( E( M-1 ) )**2
               IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+
     $             SAFMIN )GO TO 110
  100       CONTINUE
         END IF
<span class="comment">*</span><span class="comment">
</span>         M = LEND
<span class="comment">*</span><span class="comment">
</span>  110    CONTINUE
         IF( M.GT.LEND )
     $      E( M-1 ) = ZERO
         P = D( L )
         IF( M.EQ.L )
     $      GO TO 130
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        If remaining matrix is 2-by-2, use <a name="SLAE2.354"></a><a href="slae2.f.html#SLAE2.1">SLAE2</a> or <a name="SLAEV2.354"></a><a href="slaev2.f.html#SLAEV2.1">SLAEV2</a>
</span><span class="comment">*</span><span class="comment">        to compute its eigensystem.
</span><span class="comment">*</span><span class="comment">
</span>         IF( M.EQ.L-1 ) THEN
            IF( ICOMPZ.GT.0 ) THEN
               CALL <a name="SLAEV2.359"></a><a href="slaev2.f.html#SLAEV2.1">SLAEV2</a>( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S )
               WORK( M ) = C
               WORK( N-1+M ) = S
               CALL <a name="SLASR.362"></a><a href="slasr.f.html#SLASR.1">SLASR</a>( <span class="string">'R'</span>, <span class="string">'V'</span>, <span class="string">'F'</span>, N, 2, WORK( M ),
     $                     WORK( N-1+M ), Z( 1, L-1 ), LDZ )
            ELSE
               CALL <a name="SLAE2.365"></a><a href="slae2.f.html#SLAE2.1">SLAE2</a>( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 )
            END IF
            D( L-1 ) = RT1
            D( L ) = RT2
            E( L-1 ) = ZERO
            L = L - 2
            IF( L.GE.LEND )
     $         GO TO 90
            GO TO 140
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( JTOT.EQ.NMAXIT )
     $      GO TO 140
         JTOT = JTOT + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Form shift.
</span><span class="comment">*</span><span class="comment">
</span>         G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) )
         R = <a name="SLAPY2.383"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( G, ONE )
         G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) )
<span class="comment">*</span><span class="comment">
</span>         S = ONE
         C = ONE
         P = ZERO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Inner loop
</span><span class="comment">*</span><span class="comment">
</span>         LM1 = L - 1
         DO 120 I = M, LM1
            F = S*E( I )
            B = C*E( I )
            CALL <a name="SLARTG.396"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( G, F, C, S, R )
            IF( I.NE.M )
     $         E( I-1 ) = R
            G = D( I ) - P
            R = ( D( I+1 )-G )*S + TWO*C*B
            P = S*R
            D( I ) = G + P
            G = C*R - B
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           If eigenvectors are desired, then save rotations.
</span><span class="comment">*</span><span class="comment">
</span>            IF( ICOMPZ.GT.0 ) THEN
               WORK( I ) = C
               WORK( N-1+I ) = S
            END IF
<span class="comment">*</span><span class="comment">
</span>  120    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        If eigenvectors are desired, then apply saved rotations.
</span><span class="comment">*</span><span class="comment">
</span>         IF( ICOMPZ.GT.0 ) THEN
            MM = L - M + 1
            CALL <a name="SLASR.418"></a><a href="slasr.f.html#SLASR.1">SLASR</a>( <span class="string">'R'</span>, <span class="string">'V'</span>, <span class="string">'F'</span>, N, MM, WORK( M ), WORK( N-1+M ),
     $                  Z( 1, M ), LDZ )
         END IF
<span class="comment">*</span><span class="comment">
</span>         D( L ) = D( L ) - P
         E( LM1 ) = G
         GO TO 90
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Eigenvalue found.
</span><span class="comment">*</span><span class="comment">
</span>  130    CONTINUE
         D( L ) = P
<span class="comment">*</span><span class="comment">
</span>         L = L - 1
         IF( L.GE.LEND )
     $      GO TO 90
         GO TO 140
<span class="comment">*</span><span class="comment">
</span>      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Undo scaling if necessary
</span><span class="comment">*</span><span class="comment">
</span>  140 CONTINUE
      IF( ISCALE.EQ.1 ) THEN
         CALL <a name="SLASCL.442"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
     $                D( LSV ), N, INFO )
         CALL <a name="SLASCL.444"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ),
     $                N, INFO )
      ELSE IF( ISCALE.EQ.2 ) THEN
         CALL <a name="SLASCL.447"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
     $                D( LSV ), N, INFO )
         CALL <a name="SLASCL.449"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ),
     $                N, INFO )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Check for no convergence to an eigenvalue after a total
</span><span class="comment">*</span><span class="comment">     of N*MAXIT iterations.
</span><span class="comment">*</span><span class="comment">
</span>      IF( JTOT.LT.NMAXIT )
     $   GO TO 10
      DO 150 I = 1, N - 1
         IF( E( I ).NE.ZERO )
     $      INFO = INFO + 1
  150 CONTINUE
      GO TO 190
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Order eigenvalues and eigenvectors.
</span><span class="comment">*</span><span class="comment">
</span>  160 CONTINUE
      IF( ICOMPZ.EQ.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Use Quick Sort
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SLASRT.471"></a><a href="slasrt.f.html#SLASRT.1">SLASRT</a>( <span class="string">'I'</span>, N, D, INFO )
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Use Selection Sort to minimize swaps of eigenvectors
</span><span class="comment">*</span><span class="comment">
</span>         DO 180 II = 2, N
            I = II - 1
            K = I
            P = D( I )
            DO 170 J = II, N
               IF( D( J ).LT.P ) THEN
                  K = J
                  P = D( J )
               END IF
  170       CONTINUE
            IF( K.NE.I ) THEN
               D( K ) = D( I )
               D( I ) = P
               CALL SSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
            END IF
  180    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span>  190 CONTINUE
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="SSTEQR.498"></a><a href="ssteqr.f.html#SSTEQR.1">SSTEQR</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>