slasd7.f.html

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

<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Deflate due to small z component.
</span><span class="comment">*</span><span class="comment">
</span>            K2 = K2 - 1
            IDXP( K2 ) = J
            IF( J.EQ.N )
     $         GO TO 100
         ELSE
            JPREV = J
            GO TO 70
         END IF
   60 CONTINUE
   70 CONTINUE
      J = JPREV
   80 CONTINUE
      J = J + 1
      IF( J.GT.N )
     $   GO TO 90
      IF( ABS( Z( J ) ).LE.TOL ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Deflate due to small z component.
</span><span class="comment">*</span><span class="comment">
</span>         K2 = K2 - 1
         IDXP( K2 ) = J
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Check if singular values are close enough to allow deflation.
</span><span class="comment">*</span><span class="comment">
</span>         IF( ABS( D( J )-D( JPREV ) ).LE.TOL ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Deflation is possible.
</span><span class="comment">*</span><span class="comment">
</span>            S = Z( JPREV )
            C = Z( J )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Find sqrt(a**2+b**2) without overflow or
</span><span class="comment">*</span><span class="comment">           destructive underflow.
</span><span class="comment">*</span><span class="comment">
</span>            TAU = <a name="SLAPY2.334"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( C, S )
            Z( J ) = TAU
            Z( JPREV ) = ZERO
            C = C / TAU
            S = -S / TAU
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Record the appropriate Givens rotation
</span><span class="comment">*</span><span class="comment">
</span>            IF( ICOMPQ.EQ.1 ) THEN
               GIVPTR = GIVPTR + 1
               IDXJP = IDXQ( IDX( JPREV )+1 )
               IDXJ = IDXQ( IDX( J )+1 )
               IF( IDXJP.LE.NLP1 ) THEN
                  IDXJP = IDXJP - 1
               END IF
               IF( IDXJ.LE.NLP1 ) THEN
                  IDXJ = IDXJ - 1
               END IF
               GIVCOL( GIVPTR, 2 ) = IDXJP
               GIVCOL( GIVPTR, 1 ) = IDXJ
               GIVNUM( GIVPTR, 2 ) = C
               GIVNUM( GIVPTR, 1 ) = S
            END IF
            CALL SROT( 1, VF( JPREV ), 1, VF( J ), 1, C, S )
            CALL SROT( 1, VL( JPREV ), 1, VL( J ), 1, C, S )
            K2 = K2 - 1
            IDXP( K2 ) = JPREV
            JPREV = J
         ELSE
            K = K + 1
            ZW( K ) = Z( JPREV )
            DSIGMA( K ) = D( JPREV )
            IDXP( K ) = JPREV
            JPREV = J
         END IF
      END IF
      GO TO 80
   90 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Record the last singular value.
</span><span class="comment">*</span><span class="comment">
</span>      K = K + 1
      ZW( K ) = Z( JPREV )
      DSIGMA( K ) = D( JPREV )
      IDXP( K ) = JPREV
<span class="comment">*</span><span class="comment">
</span>  100 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Sort the singular values into DSIGMA. The singular values which
</span><span class="comment">*</span><span class="comment">     were not deflated go into the first K slots of DSIGMA, except
</span><span class="comment">*</span><span class="comment">     that DSIGMA(1) is treated separately.
</span><span class="comment">*</span><span class="comment">
</span>      DO 110 J = 2, N
         JP = IDXP( J )
         DSIGMA( J ) = D( JP )
         VFW( J ) = VF( JP )
         VLW( J ) = VL( JP )
  110 CONTINUE
      IF( ICOMPQ.EQ.1 ) THEN
         DO 120 J = 2, N
            JP = IDXP( J )
            PERM( J ) = IDXQ( IDX( JP )+1 )
            IF( PERM( J ).LE.NLP1 ) THEN
               PERM( J ) = PERM( J ) - 1
            END IF
  120    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     The deflated singular values go back into the last N - K slots of
</span><span class="comment">*</span><span class="comment">     D.
</span><span class="comment">*</span><span class="comment">
</span>      CALL SCOPY( N-K, DSIGMA( K+1 ), 1, D( K+1 ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and
</span><span class="comment">*</span><span class="comment">     VL(M).
</span><span class="comment">*</span><span class="comment">
</span>      DSIGMA( 1 ) = ZERO
      HLFTOL = TOL / TWO
      IF( ABS( DSIGMA( 2 ) ).LE.HLFTOL )
     $   DSIGMA( 2 ) = HLFTOL
      IF( M.GT.N ) THEN
         Z( 1 ) = <a name="SLAPY2.415"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( Z1, Z( M ) )
         IF( Z( 1 ).LE.TOL ) THEN
            C = ONE
            S = ZERO
            Z( 1 ) = TOL
         ELSE
            C = Z1 / Z( 1 )
            S = -Z( M ) / Z( 1 )
         END IF
         CALL SROT( 1, VF( M ), 1, VF( 1 ), 1, C, S )
         CALL SROT( 1, VL( M ), 1, VL( 1 ), 1, C, S )
      ELSE
         IF( ABS( Z1 ).LE.TOL ) THEN
            Z( 1 ) = TOL
         ELSE
            Z( 1 ) = Z1
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Restore Z, VF, and VL.
</span><span class="comment">*</span><span class="comment">
</span>      CALL SCOPY( K-1, ZW( 2 ), 1, Z( 2 ), 1 )
      CALL SCOPY( N-1, VFW( 2 ), 1, VF( 2 ), 1 )
      CALL SCOPY( N-1, VLW( 2 ), 1, VL( 2 ), 1 )
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="SLASD7.442"></a><a href="slasd7.f.html#SLASD7.1">SLASD7</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>