sbdsqr.f.html

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

</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SLASV2.367"></a><a href="slasv2.f.html#SLASV2.1">SLASV2</a>( D( M-1 ), E( M-1 ), D( M ), SIGMN, SIGMX, SINR,
     $                COSR, SINL, COSL )
         D( M-1 ) = SIGMX
         E( M-1 ) = ZERO
         D( M ) = SIGMN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute singular vectors, if desired
</span><span class="comment">*</span><span class="comment">
</span>         IF( NCVT.GT.0 )
     $      CALL SROT( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT, COSR,
     $                 SINR )
         IF( NRU.GT.0 )
     $      CALL SROT( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, SINL )
         IF( NCC.GT.0 )
     $      CALL SROT( NCC, C( M-1, 1 ), LDC, C( M, 1 ), LDC, COSL,
     $                 SINL )
         M = M - 2
         GO TO 60
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If working on new submatrix, choose shift direction
</span><span class="comment">*</span><span class="comment">     (from larger end diagonal element towards smaller)
</span><span class="comment">*</span><span class="comment">
</span>      IF( LL.GT.OLDM .OR. M.LT.OLDLL ) THEN
         IF( ABS( D( LL ) ).GE.ABS( D( M ) ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Chase bulge from top (big end) to bottom (small end)
</span><span class="comment">*</span><span class="comment">
</span>            IDIR = 1
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Chase bulge from bottom (big end) to top (small end)
</span><span class="comment">*</span><span class="comment">
</span>            IDIR = 2
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Apply convergence tests
</span><span class="comment">*</span><span class="comment">
</span>      IF( IDIR.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Run convergence test in forward direction
</span><span class="comment">*</span><span class="comment">        First apply standard test to bottom of matrix
</span><span class="comment">*</span><span class="comment">
</span>         IF( ABS( E( M-1 ) ).LE.ABS( TOL )*ABS( D( M ) ) .OR.
     $       ( TOL.LT.ZERO .AND. ABS( E( M-1 ) ).LE.THRESH ) ) THEN
            E( M-1 ) = ZERO
            GO TO 60
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( TOL.GE.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           If relative accuracy desired,
</span><span class="comment">*</span><span class="comment">           apply convergence criterion forward
</span><span class="comment">*</span><span class="comment">
</span>            MU = ABS( D( LL ) )
            SMINL = MU
            DO 100 LLL = LL, M - 1
               IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
                  E( LLL ) = ZERO
                  GO TO 60
               END IF
               MU = ABS( D( LLL+1 ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
               SMINL = MIN( SMINL, MU )
  100       CONTINUE
         END IF
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Run convergence test in backward direction
</span><span class="comment">*</span><span class="comment">        First apply standard test to top of matrix
</span><span class="comment">*</span><span class="comment">
</span>         IF( ABS( E( LL ) ).LE.ABS( TOL )*ABS( D( LL ) ) .OR.
     $       ( TOL.LT.ZERO .AND. ABS( E( LL ) ).LE.THRESH ) ) THEN
            E( LL ) = ZERO
            GO TO 60
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( TOL.GE.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           If relative accuracy desired,
</span><span class="comment">*</span><span class="comment">           apply convergence criterion backward
</span><span class="comment">*</span><span class="comment">
</span>            MU = ABS( D( M ) )
            SMINL = MU
            DO 110 LLL = M - 1, LL, -1
               IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
                  E( LLL ) = ZERO
                  GO TO 60
               END IF
               MU = ABS( D( LLL ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
               SMINL = MIN( SMINL, MU )
  110       CONTINUE
         END IF
      END IF
      OLDLL = LL
      OLDM = M
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute shift.  First, test if shifting would ruin relative
</span><span class="comment">*</span><span class="comment">     accuracy, and if so set the shift to zero.
</span><span class="comment">*</span><span class="comment">
</span>      IF( TOL.GE.ZERO .AND. N*TOL*( SMINL / SMAX ).LE.
     $    MAX( EPS, HNDRTH*TOL ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Use a zero shift to avoid loss of relative accuracy
</span><span class="comment">*</span><span class="comment">
</span>         SHIFT = ZERO
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute the shift from 2-by-2 block at end of matrix
</span><span class="comment">*</span><span class="comment">
</span>         IF( IDIR.EQ.1 ) THEN
            SLL = ABS( D( LL ) )
            CALL <a name="SLAS2.480"></a><a href="slas2.f.html#SLAS2.1">SLAS2</a>( D( M-1 ), E( M-1 ), D( M ), SHIFT, R )
         ELSE
            SLL = ABS( D( M ) )
            CALL <a name="SLAS2.483"></a><a href="slas2.f.html#SLAS2.1">SLAS2</a>( D( LL ), E( LL ), D( LL+1 ), SHIFT, R )
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Test if shift negligible, and if so set to zero
</span><span class="comment">*</span><span class="comment">
</span>         IF( SLL.GT.ZERO ) THEN
            IF( ( SHIFT / SLL )**2.LT.EPS )
     $         SHIFT = ZERO
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Increment iteration count
</span><span class="comment">*</span><span class="comment">
</span>      ITER = ITER + M - LL
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If SHIFT = 0, do simplified QR iteration
</span><span class="comment">*</span><span class="comment">
</span>      IF( SHIFT.EQ.ZERO ) THEN
         IF( IDIR.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Chase bulge from top to bottom
</span><span class="comment">*</span><span class="comment">           Save cosines and sines for later singular vector updates
</span><span class="comment">*</span><span class="comment">
</span>            CS = ONE
            OLDCS = ONE
            DO 120 I = LL, M - 1
               CALL <a name="SLARTG.509"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( D( I )*CS, E( I ), CS, SN, R )
               IF( I.GT.LL )
     $            E( I-1 ) = OLDSN*R
               CALL <a name="SLARTG.512"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( OLDCS*R, D( I+1 )*SN, OLDCS, OLDSN, D( I ) )
               WORK( I-LL+1 ) = CS
               WORK( I-LL+1+NM1 ) = SN
               WORK( I-LL+1+NM12 ) = OLDCS
               WORK( I-LL+1+NM13 ) = OLDSN
  120       CONTINUE
            H = D( M )*CS
            D( M ) = H*OLDCS
            E( M-1 ) = H*OLDSN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Update singular vectors
</span><span class="comment">*</span><span class="comment">
</span>            IF( NCVT.GT.0 )
     $         CALL <a name="SLASR.525"></a><a href="slasr.f.html#SLASR.1">SLASR</a>( <span class="string">'L'</span>, <span class="string">'V'</span>, <span class="string">'F'</span>, M-LL+1, NCVT, WORK( 1 ),
     $                     WORK( N ), VT( LL, 1 ), LDVT )
            IF( NRU.GT.0 )
     $         CALL <a name="SLASR.528"></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>, NRU, M-LL+1, WORK( NM12+1 ),
     $                     WORK( NM13+1 ), U( 1, LL ), LDU )
            IF( NCC.GT.0 )
     $         CALL <a name="SLASR.531"></a><a href="slasr.f.html#SLASR.1">SLASR</a>( <span class="string">'L'</span>, <span class="string">'V'</span>, <span class="string">'F'</span>, M-LL+1, NCC, WORK( NM12+1 ),
     $                     WORK( NM13+1 ), C( LL, 1 ), LDC )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Test convergence
</span><span class="comment">*</span><span class="comment">
</span>            IF( ABS( E( M-1 ) ).LE.THRESH )
     $         E( M-1 ) = ZERO
<span class="comment">*</span><span class="comment">
</span>         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Chase bulge from bottom to top
</span><span class="comment">*</span><span class="comment">           Save cosines and sines for later singular vector updates
</span><span class="comment">*</span><span class="comment">
</span>            CS = ONE
            OLDCS = ONE
            DO 130 I = M, LL + 1, -1
               CALL <a name="SLARTG.547"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( D( I )*CS, E( I-1 ), CS, SN, R )
               IF( I.LT.M )
     $            E( I ) = OLDSN*R
               CALL <a name="SLARTG.550"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( OLDCS*R, D( I-1 )*SN, OLDCS, OLDSN, D( I ) )
               WORK( I-LL ) = CS
               WORK( I-LL+NM1 ) = -SN
               WORK( I-LL+NM12 ) = OLDCS
               WORK( I-LL+NM13 ) = -OLDSN
  130       CONTINUE
            H = D( LL )*CS
            D( LL ) = H*OLDCS