dlasd2.f.html

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

</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Based on contributions by
</span><span class="comment">*</span><span class="comment">     Ming Gu and Huan Ren, Computer Science Division, University of
</span><span class="comment">*</span><span class="comment">     California at Berkeley, USA
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      DOUBLE PRECISION   ZERO, ONE, TWO, EIGHT
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
     $                   EIGHT = 8.0D+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      INTEGER            CTOT( 4 ), PSM( 4 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      INTEGER            CT, I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M,
     $                   N, NLP1, NLP2
      DOUBLE PRECISION   C, EPS, HLFTOL, S, TAU, TOL, Z1
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      DOUBLE PRECISION   <a name="DLAMCH.183"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLAPY2.183"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>
      EXTERNAL           <a name="DLAMCH.184"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLAPY2.184"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           DCOPY, <a name="DLACPY.187"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>, <a name="DLAMRG.187"></a><a href="dlamrg.f.html#DLAMRG.1">DLAMRG</a>, <a name="DLASET.187"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>, DROT, <a name="XERBLA.187"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, MAX
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Test the input parameters.
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
<span class="comment">*</span><span class="comment">
</span>      IF( NL.LT.1 ) THEN
         INFO = -1
      ELSE IF( NR.LT.1 ) THEN
         INFO = -2
      ELSE IF( ( SQRE.NE.1 ) .AND. ( SQRE.NE.0 ) ) THEN
         INFO = -3
      END IF
<span class="comment">*</span><span class="comment">
</span>      N = NL + NR + 1
      M = N + SQRE
<span class="comment">*</span><span class="comment">
</span>      IF( LDU.LT.N ) THEN
         INFO = -10
      ELSE IF( LDVT.LT.M ) THEN
         INFO = -12
      ELSE IF( LDU2.LT.N ) THEN
         INFO = -15
      ELSE IF( LDVT2.LT.M ) THEN
         INFO = -17
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.219"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DLASD2.219"></a><a href="dlasd2.f.html#DLASD2.1">DLASD2</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      NLP1 = NL + 1
      NLP2 = NL + 2
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Generate the first part of the vector Z; and move the singular
</span><span class="comment">*</span><span class="comment">     values in the first part of D one position backward.
</span><span class="comment">*</span><span class="comment">
</span>      Z1 = ALPHA*VT( NLP1, NLP1 )
      Z( 1 ) = Z1
      DO 10 I = NL, 1, -1
         Z( I+1 ) = ALPHA*VT( I, NLP1 )
         D( I+1 ) = D( I )
         IDXQ( I+1 ) = IDXQ( I ) + 1
   10 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Generate the second part of the vector Z.
</span><span class="comment">*</span><span class="comment">
</span>      DO 20 I = NLP2, M
         Z( I ) = BETA*VT( I, NLP2 )
   20 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initialize some reference arrays.
</span><span class="comment">*</span><span class="comment">
</span>      DO 30 I = 2, NLP1
         COLTYP( I ) = 1
   30 CONTINUE
      DO 40 I = NLP2, N
         COLTYP( I ) = 2
   40 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Sort the singular values into increasing order
</span><span class="comment">*</span><span class="comment">
</span>      DO 50 I = NLP2, N
         IDXQ( I ) = IDXQ( I ) + NLP1
   50 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     DSIGMA, IDXC, IDXC, and the first column of U2
</span><span class="comment">*</span><span class="comment">     are used as storage space.
</span><span class="comment">*</span><span class="comment">
</span>      DO 60 I = 2, N
         DSIGMA( I ) = D( IDXQ( I ) )
         U2( I, 1 ) = Z( IDXQ( I ) )
         IDXC( I ) = COLTYP( IDXQ( I ) )
   60 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      CALL <a name="DLAMRG.267"></a><a href="dlamrg.f.html#DLAMRG.1">DLAMRG</a>( NL, NR, DSIGMA( 2 ), 1, 1, IDX( 2 ) )
<span class="comment">*</span><span class="comment">
</span>      DO 70 I = 2, N
         IDXI = 1 + IDX( I )
         D( I ) = DSIGMA( IDXI )
         Z( I ) = U2( IDXI, 1 )
         COLTYP( I ) = IDXC( IDXI )
   70 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Calculate the allowable deflation tolerance
</span><span class="comment">*</span><span class="comment">
</span>      EPS = <a name="DLAMCH.278"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Epsilon'</span> )
      TOL = MAX( ABS( ALPHA ), ABS( BETA ) )
      TOL = EIGHT*EPS*MAX( ABS( D( N ) ), TOL )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     There are 2 kinds of deflation -- first a value in the z-vector
</span><span class="comment">*</span><span class="comment">     is small, second two (or more) singular values are very close
</span><span class="comment">*</span><span class="comment">     together (their difference is small).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If the value in the z-vector is small, we simply permute the
</span><span class="comment">*</span><span class="comment">     array so that the corresponding singular value is moved to the
</span><span class="comment">*</span><span class="comment">     end.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If two values in the D-vector are close, we perform a two-sided
</span><span class="comment">*</span><span class="comment">     rotation designed to make one of the corresponding z-vector
</span><span class="comment">*</span><span class="comment">     entries zero, and then permute the array so that the deflated
</span><span class="comment">*</span><span class="comment">     singular value is moved to the end.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If there are multiple singular values then the problem deflates.
</span><span class="comment">*</span><span class="comment">     Here the number of equal singular values are found.  As each equal
</span><span class="comment">*</span><span class="comment">     singular value is found, an elementary reflector is computed to
</span><span class="comment">*</span><span class="comment">     rotate the corresponding singular subspace so that the
</span><span class="comment">*</span><span class="comment">     corresponding components of Z are zero in this new basis.
</span><span class="comment">*</span><span class="comment">
</span>      K = 1
      K2 = N + 1
      DO 80 J = 2, N
         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
            COLTYP( J ) = 4
            IF( J.EQ.N )
     $         GO TO 120
         ELSE
            JPREV = J
            GO TO 90
         END IF
   80 CONTINUE
   90 CONTINUE
      J = JPREV
  100 CONTINUE
      J = J + 1
      IF( J.GT.N )
     $   GO TO 110
      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
         COLTYP( J ) = 4
      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 )