sbdsqr.f.html

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

      REAL               <a name="SLAMCH.174"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
      EXTERNAL           <a name="LSAME.175"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="SLAMCH.175"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="SLARTG.178"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>, <a name="SLAS2.178"></a><a href="slas2.f.html#SLAS2.1">SLAS2</a>, <a name="SLASQ1.178"></a><a href="slasq1.f.html#SLASQ1.1">SLASQ1</a>, <a name="SLASR.178"></a><a href="slasr.f.html#SLASR.1">SLASR</a>, <a name="SLASV2.178"></a><a href="slasv2.f.html#SLASV2.1">SLASV2</a>, SROT,
     $                   SSCAL, SSWAP, <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">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, MAX, MIN, REAL, SIGN, SQRT
<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
      LOWER = <a name="LSAME.189"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'L'</span> )
      IF( .NOT.<a name="LSAME.190"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'U'</span> ) .AND. .NOT.LOWER ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( NCVT.LT.0 ) THEN
         INFO = -3
      ELSE IF( NRU.LT.0 ) THEN
         INFO = -4
      ELSE IF( NCC.LT.0 ) THEN
         INFO = -5
      ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR.
     $         ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN
         INFO = -9
      ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN
         INFO = -11
      ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR.
     $         ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN
         INFO = -13
      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="SBDSQR.210"></a><a href="sbdsqr.f.html#SBDSQR.1">SBDSQR</a>'</span>, -INFO )
         RETURN
      END IF
      IF( N.EQ.0 )
     $   RETURN
      IF( N.EQ.1 )
     $   GO TO 160
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ROTATE is true if any singular vectors desired, false otherwise
</span><span class="comment">*</span><span class="comment">
</span>      ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If no singular vectors desired, use qd algorithm
</span><span class="comment">*</span><span class="comment">
</span>      IF( .NOT.ROTATE ) THEN
         CALL <a name="SLASQ1.225"></a><a href="slasq1.f.html#SLASQ1.1">SLASQ1</a>( N, D, E, WORK, INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      NM1 = N - 1
      NM12 = NM1 + NM1
      NM13 = NM12 + NM1
      IDIR = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Get machine constants
</span><span class="comment">*</span><span class="comment">
</span>      EPS = <a name="SLAMCH.236"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'Epsilon'</span> )
      UNFL = <a name="SLAMCH.237"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'Safe minimum'</span> )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If matrix lower bidiagonal, rotate to be upper bidiagonal
</span><span class="comment">*</span><span class="comment">     by applying Givens rotations on the left
</span><span class="comment">*</span><span class="comment">
</span>      IF( LOWER ) THEN
         DO 10 I = 1, N - 1
            CALL <a name="SLARTG.244"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( D( I ), E( I ), CS, SN, R )
            D( I ) = R
            E( I ) = SN*D( I+1 )
            D( I+1 ) = CS*D( I+1 )
            WORK( I ) = CS
            WORK( NM1+I ) = SN
   10    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Update singular vectors if desired
</span><span class="comment">*</span><span class="comment">
</span>         IF( NRU.GT.0 )
     $      CALL <a name="SLASR.255"></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, N, WORK( 1 ), WORK( N ), U,
     $                  LDU )
         IF( NCC.GT.0 )
     $      CALL <a name="SLASR.258"></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>, N, NCC, WORK( 1 ), WORK( N ), C,
     $                  LDC )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute singular values to relative accuracy TOL
</span><span class="comment">*</span><span class="comment">     (By setting TOL to be negative, algorithm will compute
</span><span class="comment">*</span><span class="comment">     singular values to absolute accuracy ABS(TOL)*norm(input matrix))
</span><span class="comment">*</span><span class="comment">
</span>      TOLMUL = MAX( TEN, MIN( HNDRD, EPS**MEIGTH ) )
      TOL = TOLMUL*EPS
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute approximate maximum, minimum singular values
</span><span class="comment">*</span><span class="comment">
</span>      SMAX = ZERO
      DO 20 I = 1, N
         SMAX = MAX( SMAX, ABS( D( I ) ) )
   20 CONTINUE
      DO 30 I = 1, N - 1
         SMAX = MAX( SMAX, ABS( E( I ) ) )
   30 CONTINUE
      SMINL = ZERO
      IF( TOL.GE.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Relative accuracy desired
</span><span class="comment">*</span><span class="comment">
</span>         SMINOA = ABS( D( 1 ) )
         IF( SMINOA.EQ.ZERO )
     $      GO TO 50
         MU = SMINOA
         DO 40 I = 2, N
            MU = ABS( D( I ) )*( MU / ( MU+ABS( E( I-1 ) ) ) )
            SMINOA = MIN( SMINOA, MU )
            IF( SMINOA.EQ.ZERO )
     $         GO TO 50
   40    CONTINUE
   50    CONTINUE
         SMINOA = SMINOA / SQRT( REAL( N ) )
         THRESH = MAX( TOL*SMINOA, MAXITR*N*N*UNFL )
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Absolute accuracy desired
</span><span class="comment">*</span><span class="comment">
</span>         THRESH = MAX( ABS( TOL )*SMAX, MAXITR*N*N*UNFL )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Prepare for main iteration loop for the singular values
</span><span class="comment">*</span><span class="comment">     (MAXIT is the maximum number of passes through the inner
</span><span class="comment">*</span><span class="comment">     loop permitted before nonconvergence signalled.)
</span><span class="comment">*</span><span class="comment">
</span>      MAXIT = MAXITR*N*N
      ITER = 0
      OLDLL = -1
      OLDM = -1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     M points to last element of unconverged part of matrix
</span><span class="comment">*</span><span class="comment">
</span>      M = N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Begin main iteration loop
</span><span class="comment">*</span><span class="comment">
</span>   60 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Check for convergence or exceeding iteration count
</span><span class="comment">*</span><span class="comment">
</span>      IF( M.LE.1 )
     $   GO TO 160
      IF( ITER.GT.MAXIT )
     $   GO TO 200
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Find diagonal block of matrix to work on
</span><span class="comment">*</span><span class="comment">
</span>      IF( TOL.LT.ZERO .AND. ABS( D( M ) ).LE.THRESH )
     $   D( M ) = ZERO
      SMAX = ABS( D( M ) )
      SMIN = SMAX
      DO 70 LLL = 1, M - 1
         LL = M - LLL
         ABSS = ABS( D( LL ) )
         ABSE = ABS( E( LL ) )
         IF( TOL.LT.ZERO .AND. ABSS.LE.THRESH )
     $      D( LL ) = ZERO
         IF( ABSE.LE.THRESH )
     $      GO TO 80
         SMIN = MIN( SMIN, ABSS )
         SMAX = MAX( SMAX, ABSS, ABSE )
   70 CONTINUE
      LL = 0
      GO TO 90
   80 CONTINUE
      E( LL ) = ZERO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Matrix splits since E(LL) = 0
</span><span class="comment">*</span><span class="comment">
</span>      IF( LL.EQ.M-1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Convergence of bottom singular value, return to top of loop
</span><span class="comment">*</span><span class="comment">
</span>         M = M - 1
         GO TO 60
      END IF
   90 CONTINUE
      LL = LL + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     E(LL) through E(M-1) are nonzero, E(LL-1) is zero
</span><span class="comment">*</span><span class="comment">
</span>      IF( LL.EQ.M-1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        2 by 2 block, handle separately