ssterf.f.html

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

     $      GO TO 90
<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.179"></a><a href="slae2.f.html#SLAE2.1">SLAE2</a> to compute its
</span><span class="comment">*</span><span class="comment">        eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span>         IF( M.EQ.L+1 ) THEN
            RTE = SQRT( E( L ) )
            CALL <a name="SLAE2.184"></a><a href="slae2.f.html#SLAE2.1">SLAE2</a>( D( L ), RTE, D( L+1 ), RT1, RT2 )
            D( L ) = RT1
            D( L+1 ) = RT2
            E( L ) = ZERO
            L = L + 2
            IF( L.LE.LEND )
     $         GO TO 50
            GO TO 150
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( JTOT.EQ.NMAXIT )
     $      GO TO 150
         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>         RTE = SQRT( E( L ) )
         SIGMA = ( D( L+1 )-P ) / ( TWO*RTE )
         R = <a name="SLAPY2.202"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( SIGMA, ONE )
         SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
<span class="comment">*</span><span class="comment">
</span>         C = ONE
         S = ZERO
         GAMMA = D( M ) - SIGMA
         P = GAMMA*GAMMA
<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>         DO 80 I = M - 1, L, -1
            BB = E( I )
            R = P + BB
            IF( I.NE.M-1 )
     $         E( I+1 ) = S*R
            OLDC = C
            C = P / R
            S = BB / R
            OLDGAM = GAMMA
            ALPHA = D( I )
            GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
            D( I+1 ) = OLDGAM + ( ALPHA-GAMMA )
            IF( C.NE.ZERO ) THEN
               P = ( GAMMA*GAMMA ) / C
            ELSE
               P = OLDC*BB
            END IF
   80    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         E( L ) = S*P
         D( L ) = SIGMA + GAMMA
         GO TO 50
<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>   90    CONTINUE
         D( L ) = P
<span class="comment">*</span><span class="comment">
</span>         L = L + 1
         IF( L.LE.LEND )
     $      GO TO 50
         GO TO 150
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        QR Iteration
</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>  100    CONTINUE
         DO 110 M = L, LEND + 1, -1
            IF( ABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) )
     $         GO TO 120
  110    CONTINUE
         M = LEND
<span class="comment">*</span><span class="comment">
</span>  120    CONTINUE
         IF( M.GT.LEND )
     $      E( M-1 ) = ZERO
         P = D( L )
         IF( M.EQ.L )
     $      GO TO 140
<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.265"></a><a href="slae2.f.html#SLAE2.1">SLAE2</a> to compute its
</span><span class="comment">*</span><span class="comment">        eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span>         IF( M.EQ.L-1 ) THEN
            RTE = SQRT( E( L-1 ) )
            CALL <a name="SLAE2.270"></a><a href="slae2.f.html#SLAE2.1">SLAE2</a>( D( L ), RTE, D( L-1 ), RT1, RT2 )
            D( L ) = RT1
            D( L-1 ) = RT2
            E( L-1 ) = ZERO
            L = L - 2
            IF( L.GE.LEND )
     $         GO TO 100
            GO TO 150
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( JTOT.EQ.NMAXIT )
     $      GO TO 150
         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>         RTE = SQRT( E( L-1 ) )
         SIGMA = ( D( L-1 )-P ) / ( TWO*RTE )
         R = <a name="SLAPY2.288"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( SIGMA, ONE )
         SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
<span class="comment">*</span><span class="comment">
</span>         C = ONE
         S = ZERO
         GAMMA = D( M ) - SIGMA
         P = GAMMA*GAMMA
<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>         DO 130 I = M, L - 1
            BB = E( I )
            R = P + BB
            IF( I.NE.M )
     $         E( I-1 ) = S*R
            OLDC = C
            C = P / R
            S = BB / R
            OLDGAM = GAMMA
            ALPHA = D( I+1 )
            GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
            D( I ) = OLDGAM + ( ALPHA-GAMMA )
            IF( C.NE.ZERO ) THEN
               P = ( GAMMA*GAMMA ) / C
            ELSE
               P = OLDC*BB
            END IF
  130    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         E( L-1 ) = S*P
         D( L ) = SIGMA + GAMMA
         GO TO 100
<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>  140    CONTINUE
         D( L ) = P
<span class="comment">*</span><span class="comment">
</span>         L = L - 1
         IF( L.GE.LEND )
     $      GO TO 100
         GO TO 150
<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>  150 CONTINUE
      IF( ISCALE.EQ.1 )
     $   CALL <a name="SLASCL.337"></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 )
      IF( ISCALE.EQ.2 )
     $   CALL <a name="SLASCL.340"></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 )
<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 160 I = 1, N - 1
         IF( E( I ).NE.ZERO )
     $      INFO = INFO + 1
  160 CONTINUE
      GO TO 180
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Sort eigenvalues in increasing order.
</span><span class="comment">*</span><span class="comment">
</span>  170 CONTINUE
      CALL <a name="SLASRT.357"></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>  180 CONTINUE
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="SSTERF.362"></a><a href="ssterf.f.html#SSTERF.1">SSTERF</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>