ssteqr.f.html

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

</span>      NMAXIT = N*MAXIT
      JTOT = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Determine where the matrix splits and choose QL or QR iteration
</span><span class="comment">*</span><span class="comment">     for each block, according to whether top or bottom diagonal
</span><span class="comment">*</span><span class="comment">     element is smaller.
</span><span class="comment">*</span><span class="comment">
</span>      L1 = 1
      NM1 = N - 1
<span class="comment">*</span><span class="comment">
</span>   10 CONTINUE
      IF( L1.GT.N )
     $   GO TO 160
      IF( L1.GT.1 )
     $   E( L1-1 ) = ZERO
      IF( L1.LE.NM1 ) THEN
         DO 20 M = L1, NM1
            TST = ABS( E( M ) )
            IF( TST.EQ.ZERO )
     $         GO TO 30
            IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
     $          1 ) ) ) )*EPS ) THEN
               E( M ) = ZERO
               GO TO 30
            END IF
   20    CONTINUE
      END IF
      M = N
<span class="comment">*</span><span class="comment">
</span>   30 CONTINUE
      L = L1
      LSV = L
      LEND = M
      LENDSV = LEND
      L1 = M + 1
      IF( LEND.EQ.L )
     $   GO TO 10
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Scale submatrix in rows and columns L to LEND
</span><span class="comment">*</span><span class="comment">
</span>      ANORM = <a name="SLANST.197"></a><a href="slanst.f.html#SLANST.1">SLANST</a>( <span class="string">'I'</span>, LEND-L+1, D( L ), E( L ) )
      ISCALE = 0
      IF( ANORM.EQ.ZERO )
     $   GO TO 10
      IF( ANORM.GT.SSFMAX ) THEN
         ISCALE = 1
         CALL <a name="SLASCL.203"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
     $                INFO )
         CALL <a name="SLASCL.205"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
     $                INFO )
      ELSE IF( ANORM.LT.SSFMIN ) THEN
         ISCALE = 2
         CALL <a name="SLASCL.209"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
     $                INFO )
         CALL <a name="SLASCL.211"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
     $                INFO )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Choose between QL and QR iteration
</span><span class="comment">*</span><span class="comment">
</span>      IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
         LEND = LSV
         L = LENDSV
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( LEND.GT.L ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        QL Iteration
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Look for small subdiagonal element.
</span><span class="comment">*</span><span class="comment">
</span>   40    CONTINUE
         IF( L.NE.LEND ) THEN
            LENDM1 = LEND - 1
            DO 50 M = L, LENDM1
               TST = ABS( E( M ) )**2
               IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M+1 ) )+
     $             SAFMIN )GO TO 60
   50       CONTINUE
         END IF
<span class="comment">*</span><span class="comment">
</span>         M = LEND
<span class="comment">*</span><span class="comment">
</span>   60    CONTINUE
         IF( M.LT.LEND )
     $      E( M ) = ZERO
         P = D( L )
         IF( M.EQ.L )
     $      GO TO 80
<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.247"></a><a href="slae2.f.html#SLAE2.1">SLAE2</a> or <a name="SLAEV2.247"></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.252"></a><a href="slaev2.f.html#SLAEV2.1">SLAEV2</a>( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S )
               WORK( L ) = C
               WORK( N-1+L ) = S
               CALL <a name="SLASR.255"></a><a href="slasr.f.html#SLASR.1">SLASR</a>( <span class="string">'R'</span>, <span class="string">'V'</span>, <span class="string">'B'</span>, N, 2, WORK( L ),
     $                     WORK( N-1+L ), Z( 1, L ), LDZ )
            ELSE
               CALL <a name="SLAE2.258"></a><a href="slae2.f.html#SLAE2.1">SLAE2</a>( D( L ), E( L ), D( L+1 ), RT1, RT2 )
            END IF
            D( L ) = RT1
            D( L+1 ) = RT2
            E( L ) = ZERO
            L = L + 2
            IF( L.LE.LEND )
     $         GO TO 40
            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 ) )
         R = <a name="SLAPY2.276"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( G, ONE )
         G = D( M ) - P + ( E( L ) / ( 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>         MM1 = M - 1
         DO 70 I = MM1, L, -1
            F = S*E( I )
            B = C*E( I )
            CALL <a name="SLARTG.289"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( G, F, C, S, R )
            IF( I.NE.M-1 )
     $         E( I+1 ) = R
            G = D( I+1 ) - P
            R = ( D( I )-G )*S + TWO*C*B
            P = S*R
            D( I+1 ) = 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>   70    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 = M - L + 1
            CALL <a name="SLASR.311"></a><a href="slasr.f.html#SLASR.1">SLASR</a>( <span class="string">'R'</span>, <span class="string">'V'</span>, <span class="string">'B'</span>, N, MM, WORK( L ), WORK( N-1+L ),
     $                  Z( 1, L ), LDZ )
         END IF
<span class="comment">*</span><span class="comment">
</span>         D( L ) = D( L ) - P
         E( L ) = G
         GO TO 40
<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>   80    CONTINUE
         D( L ) = P
<span class="comment">*</span><span class="comment">
</span>         L = L + 1
         IF( L.LE.LEND )
     $      GO TO 40
         GO TO 140
<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