dlasd2.f.html

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

            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="DLAPY2.345"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( C, S )
            C = C / TAU
            S = -S / TAU
            Z( J ) = TAU
            Z( JPREV ) = ZERO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Apply back the Givens rotation to the left and right
</span><span class="comment">*</span><span class="comment">           singular vector matrices.
</span><span class="comment">*</span><span class="comment">
</span>            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
            CALL DROT( N, U( 1, IDXJP ), 1, U( 1, IDXJ ), 1, C, S )
            CALL DROT( M, VT( IDXJP, 1 ), LDVT, VT( IDXJ, 1 ), LDVT, C,
     $                 S )
            IF( COLTYP( J ).NE.COLTYP( JPREV ) ) THEN
               COLTYP( J ) = 3
            END IF
            COLTYP( JPREV ) = 4
            K2 = K2 - 1
            IDXP( K2 ) = JPREV
            JPREV = J
         ELSE
            K = K + 1
            U2( K, 1 ) = Z( JPREV )
            DSIGMA( K ) = D( JPREV )
            IDXP( K ) = JPREV
            JPREV = J
         END IF
      END IF
      GO TO 100
  110 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
      U2( K, 1 ) = Z( JPREV )
      DSIGMA( K ) = D( JPREV )
      IDXP( K ) = JPREV
<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">     Count up the total number of the various types of columns, then
</span><span class="comment">*</span><span class="comment">     form a permutation which positions the four column types into
</span><span class="comment">*</span><span class="comment">     four groups of uniform structure (although one or more of these
</span><span class="comment">*</span><span class="comment">     groups may be empty).
</span><span class="comment">*</span><span class="comment">
</span>      DO 130 J = 1, 4
         CTOT( J ) = 0
  130 CONTINUE
      DO 140 J = 2, N
         CT = COLTYP( J )
         CTOT( CT ) = CTOT( CT ) + 1
  140 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     PSM(*) = Position in SubMatrix (of types 1 through 4)
</span><span class="comment">*</span><span class="comment">
</span>      PSM( 1 ) = 2
      PSM( 2 ) = 2 + CTOT( 1 )
      PSM( 3 ) = PSM( 2 ) + CTOT( 2 )
      PSM( 4 ) = PSM( 3 ) + CTOT( 3 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Fill out the IDXC array so that the permutation which it induces
</span><span class="comment">*</span><span class="comment">     will place all type-1 columns first, all type-2 columns next,
</span><span class="comment">*</span><span class="comment">     then all type-3's, and finally all type-4's, starting from the
</span><span class="comment">*</span><span class="comment">     second column. This applies similarly to the rows of VT.
</span><span class="comment">*</span><span class="comment">
</span>      DO 150 J = 2, N
         JP = IDXP( J )
         CT = COLTYP( JP )
         IDXC( PSM( CT ) ) = J
         PSM( CT ) = PSM( CT ) + 1
  150 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Sort the singular values and corresponding singular vectors into
</span><span class="comment">*</span><span class="comment">     DSIGMA, U2, and VT2 respectively.  The singular values/vectors
</span><span class="comment">*</span><span class="comment">     which were not deflated go into the first K slots of DSIGMA, U2,
</span><span class="comment">*</span><span class="comment">     and VT2 respectively, while those which were deflated go into the
</span><span class="comment">*</span><span class="comment">     last N - K slots, except that the first column/row will be treated
</span><span class="comment">*</span><span class="comment">     separately.
</span><span class="comment">*</span><span class="comment">
</span>      DO 160 J = 2, N
         JP = IDXP( J )
         DSIGMA( J ) = D( JP )
         IDXJ = IDXQ( IDX( IDXP( IDXC( J ) ) )+1 )
         IF( IDXJ.LE.NLP1 ) THEN
            IDXJ = IDXJ - 1
         END IF
         CALL DCOPY( N, U( 1, IDXJ ), 1, U2( 1, J ), 1 )
         CALL DCOPY( M, VT( IDXJ, 1 ), LDVT, VT2( J, 1 ), LDVT2 )
  160 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Determine DSIGMA(1), DSIGMA(2) and Z(1)
</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="DLAPY2.449"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</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
      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">     Move the rest of the updating row to Z.
</span><span class="comment">*</span><span class="comment">
</span>      CALL DCOPY( K-1, U2( 2, 1 ), 1, Z( 2 ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Determine the first column of U2, the first row of VT2 and the
</span><span class="comment">*</span><span class="comment">     last row of VT.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="DLASET.473"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'A'</span>, N, 1, ZERO, ZERO, U2, LDU2 )
      U2( NLP1, 1 ) = ONE
      IF( M.GT.N ) THEN
         DO 170 I = 1, NLP1
            VT( M, I ) = -S*VT( NLP1, I )
            VT2( 1, I ) = C*VT( NLP1, I )
  170    CONTINUE
         DO 180 I = NLP2, M
            VT2( 1, I ) = S*VT( M, I )
            VT( M, I ) = C*VT( M, I )
  180    CONTINUE
      ELSE
         CALL DCOPY( M, VT( NLP1, 1 ), LDVT, VT2( 1, 1 ), LDVT2 )
      END IF
      IF( M.GT.N ) THEN
         CALL DCOPY( M, VT( M, 1 ), LDVT, VT2( M, 1 ), LDVT2 )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     The deflated singular values and their corresponding vectors go
</span><span class="comment">*</span><span class="comment">     into the back of D, U, and V respectively.
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.GT.K ) THEN
         CALL DCOPY( N-K, DSIGMA( K+1 ), 1, D( K+1 ), 1 )
         CALL <a name="DLACPY.496"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'A'</span>, N, N-K, U2( 1, K+1 ), LDU2, U( 1, K+1 ),
     $                LDU )
         CALL <a name="DLACPY.498"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'A'</span>, N-K, M, VT2( K+1, 1 ), LDVT2, VT( K+1, 1 ),
     $                LDVT )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Copy CTOT into COLTYP for referencing in <a name="DLASD3.502"></a><a href="dlasd3.f.html#DLASD3.1">DLASD3</a>.
</span><span class="comment">*</span><span class="comment">
</span>      DO 190 J = 1, 4
         COLTYP( J ) = CTOT( J )
  190 CONTINUE
<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="DLASD2.510"></a><a href="dlasd2.f.html#DLASD2.1">DLASD2</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>