slasd7.f.html

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

</span><span class="comment">*</span><span class="comment">  GIVNUM (output) REAL array, dimension ( LDGNUM, 2 )
</span><span class="comment">*</span><span class="comment">         Each number indicates the C or S value to be used in the
</span><span class="comment">*</span><span class="comment">         corresponding Givens rotation. Not referenced if ICOMPQ = 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDGNUM (input) INTEGER
</span><span class="comment">*</span><span class="comment">         The leading dimension of GIVNUM, must be at least N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  C      (output) REAL
</span><span class="comment">*</span><span class="comment">         C contains garbage if SQRE =0 and the C-value of a Givens
</span><span class="comment">*</span><span class="comment">         rotation related to the right null space if SQRE = 1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  S      (output) REAL
</span><span class="comment">*</span><span class="comment">         S contains garbage if SQRE =0 and the S-value of a Givens
</span><span class="comment">*</span><span class="comment">         rotation related to the right null space if SQRE = 1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO   (output) INTEGER
</span><span class="comment">*</span><span class="comment">         = 0:  successful exit.
</span><span class="comment">*</span><span class="comment">         &lt; 0:  if INFO = -i, the i-th argument had an illegal value.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</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>      REAL               ZERO, ONE, TWO, EIGHT
      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0,
     $                   EIGHT = 8.0E+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span><span class="comment">*</span><span class="comment">
</span>      INTEGER            I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M, N,
     $                   NLP1, NLP2
      REAL               EPS, HLFTOL, TAU, TOL, Z1
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           SCOPY, <a name="SLAMRG.179"></a><a href="slamrg.f.html#SLAMRG.1">SLAMRG</a>, SROT, <a name="XERBLA.179"></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">     .. External Functions ..
</span>      REAL               <a name="SLAMCH.182"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLAPY2.182"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>
      EXTERNAL           <a name="SLAMCH.183"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLAPY2.183"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</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
      N = NL + NR + 1
      M = N + SQRE
<span class="comment">*</span><span class="comment">
</span>      IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
         INFO = -1
      ELSE IF( NL.LT.1 ) THEN
         INFO = -2
      ELSE IF( NR.LT.1 ) THEN
         INFO = -3
      ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
         INFO = -4
      ELSE IF( LDGCOL.LT.N ) THEN
         INFO = -22
      ELSE IF( LDGNUM.LT.N ) THEN
         INFO = -24
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.210"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SLASD7.210"></a><a href="slasd7.f.html#SLASD7.1">SLASD7</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      NLP1 = NL + 1
      NLP2 = NL + 2
      IF( ICOMPQ.EQ.1 ) THEN
         GIVPTR = 0
      END IF
<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*VL( NLP1 )
      VL( NLP1 ) = ZERO
      TAU = VF( NLP1 )
      DO 10 I = NL, 1, -1
         Z( I+1 ) = ALPHA*VL( I )
         VL( I ) = ZERO
         VF( I+1 ) = VF( I )
         D( I+1 ) = D( I )
         IDXQ( I+1 ) = IDXQ( I ) + 1
   10 CONTINUE
      VF( 1 ) = TAU
<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*VF( I )
         VF( I ) = ZERO
   20 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 30 I = NLP2, N
         IDXQ( I ) = IDXQ( I ) + NLP1
   30 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     DSIGMA, IDXC, IDXC, and ZW are used as storage space.
</span><span class="comment">*</span><span class="comment">
</span>      DO 40 I = 2, N
         DSIGMA( I ) = D( IDXQ( I ) )
         ZW( I ) = Z( IDXQ( I ) )
         VFW( I ) = VF( IDXQ( I ) )
         VLW( I ) = VL( IDXQ( I ) )
   40 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      CALL <a name="SLAMRG.257"></a><a href="slamrg.f.html#SLAMRG.1">SLAMRG</a>( NL, NR, DSIGMA( 2 ), 1, 1, IDX( 2 ) )
<span class="comment">*</span><span class="comment">
</span>      DO 50 I = 2, N
         IDXI = 1 + IDX( I )
         D( I ) = DSIGMA( IDXI )
         Z( I ) = ZW( IDXI )
         VF( I ) = VFW( IDXI )
         VL( I ) = VLW( IDXI )
   50 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Calculate the allowable deflation tolerence
</span><span class="comment">*</span><span class="comment">
</span>      EPS = <a name="SLAMCH.269"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'Epsilon'</span> )
      TOL = MAX( ABS( ALPHA ), ABS( BETA ) )
      TOL = EIGHT*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 60 J = 2, N
         IF( ABS( Z( J ) ).LE.TOL ) THEN