sbdsdc.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 453 行 · 第 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">     .. Parameters ..
</span>      REAL               ZERO, ONE, TWO
      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      INTEGER            DIFL, DIFR, GIVCOL, GIVNUM, GIVPTR, I, IC,
     $                   ICOMPQ, IERR, II, IS, IU, IUPLO, IVT, J, K, KK,
     $                   MLVL, NM1, NSIZE, PERM, POLES, QSTART, SMLSIZ,
     $                   SMLSZP, SQRE, START, WSTART, Z
      REAL               CS, EPS, ORGNRM, P, R, SN
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.147"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      INTEGER            <a name="ILAENV.148"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
      REAL               <a name="SLAMCH.149"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLANST.149"></a><a href="slanst.f.html#SLANST.1">SLANST</a>
      EXTERNAL           <a name="SLAMCH.150"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLANST.150"></a><a href="slanst.f.html#SLANST.1">SLANST</a>, <a name="ILAENV.150"></a><a href="hfy-index.html#ILAENV">ILAENV</a>, <a name="LSAME.150"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           SCOPY, <a name="SLARTG.153"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>, <a name="SLASCL.153"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>, <a name="SLASD0.153"></a><a href="slasd0.f.html#SLASD0.1">SLASD0</a>, <a name="SLASDA.153"></a><a href="slasda.f.html#SLASDA.1">SLASDA</a>, <a name="SLASDQ.153"></a><a href="slasdq.f.html#SLASDQ.1">SLASDQ</a>,
     $                   <a name="SLASET.154"></a><a href="slaset.f.html#SLASET.1">SLASET</a>, <a name="SLASR.154"></a><a href="slasr.f.html#SLASR.1">SLASR</a>, SSWAP, <a name="XERBLA.154"></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          REAL, ABS, INT, LOG, SIGN
<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>      IUPLO = 0
      IF( <a name="LSAME.166"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'U'</span> ) )
     $   IUPLO = 1
      IF( <a name="LSAME.168"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'L'</span> ) )
     $   IUPLO = 2
      IF( <a name="LSAME.170"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPQ, <span class="string">'N'</span> ) ) THEN
         ICOMPQ = 0
      ELSE IF( <a name="LSAME.172"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPQ, <span class="string">'P'</span> ) ) THEN
         ICOMPQ = 1
      ELSE IF( <a name="LSAME.174"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPQ, <span class="string">'I'</span> ) ) THEN
         ICOMPQ = 2
      ELSE
         ICOMPQ = -1
      END IF
      IF( IUPLO.EQ.0 ) THEN
         INFO = -1
      ELSE IF( ICOMPQ.LT.0 ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( ( LDU.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDU.LT.
     $         N ) ) ) THEN
         INFO = -7
      ELSE IF( ( LDVT.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDVT.LT.
     $         N ) ) ) THEN
         INFO = -9
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.193"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SBDSDC.193"></a><a href="sbdsdc.f.html#SBDSDC.1">SBDSDC</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.0 )
     $   RETURN
      SMLSIZ = <a name="ILAENV.201"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 9, <span class="string">'<a name="SBDSDC.201"></a><a href="sbdsdc.f.html#SBDSDC.1">SBDSDC</a>'</span>, <span class="string">' '</span>, 0, 0, 0, 0 )
      IF( N.EQ.1 ) THEN
         IF( ICOMPQ.EQ.1 ) THEN
            Q( 1 ) = SIGN( ONE, D( 1 ) )
            Q( 1+SMLSIZ*N ) = ONE
         ELSE IF( ICOMPQ.EQ.2 ) THEN
            U( 1, 1 ) = SIGN( ONE, D( 1 ) )
            VT( 1, 1 ) = ONE
         END IF
         D( 1 ) = ABS( D( 1 ) )
         RETURN
      END IF
      NM1 = N - 1
<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>      WSTART = 1
      QSTART = 3
      IF( ICOMPQ.EQ.1 ) THEN
         CALL SCOPY( N, D, 1, Q( 1 ), 1 )
         CALL SCOPY( N-1, E, 1, Q( N+1 ), 1 )
      END IF
      IF( IUPLO.EQ.2 ) THEN
         QSTART = 5
         WSTART = 2*N - 1
         DO 10 I = 1, N - 1
            CALL <a name="SLARTG.228"></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 )
            IF( ICOMPQ.EQ.1 ) THEN
               Q( I+2*N ) = CS
               Q( I+3*N ) = SN
            ELSE IF( ICOMPQ.EQ.2 ) THEN
               WORK( I ) = CS
               WORK( NM1+I ) = -SN
            END IF
   10    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If ICOMPQ = 0, use <a name="SLASDQ.242"></a><a href="slasdq.f.html#SLASDQ.1">SLASDQ</a> to compute the singular values.
</span><span class="comment">*</span><span class="comment">
</span>      IF( ICOMPQ.EQ.0 ) THEN
         CALL <a name="SLASDQ.245"></a><a href="slasdq.f.html#SLASDQ.1">SLASDQ</a>( <span class="string">'U'</span>, 0, N, 0, 0, 0, D, E, VT, LDVT, U, LDU, U,
     $                LDU, WORK( WSTART ), INFO )
         GO TO 40
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If N is smaller than the minimum divide size SMLSIZ, then solve
</span><span class="comment">*</span><span class="comment">     the problem with another solver.
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.LE.SMLSIZ ) THEN
         IF( ICOMPQ.EQ.2 ) THEN
            CALL <a name="SLASET.255"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, U, LDU )
            CALL <a name="SLASET.256"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, VT, LDVT )
            CALL <a name="SLASDQ.257"></a><a href="slasdq.f.html#SLASDQ.1">SLASDQ</a>( <span class="string">'U'</span>, 0, N, N, N, 0, D, E, VT, LDVT, U, LDU, U,
     $                   LDU, WORK( WSTART ), INFO )
         ELSE IF( ICOMPQ.EQ.1 ) THEN
            IU = 1
            IVT = IU + N
            CALL <a name="SLASET.262"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, Q( IU+( QSTART-1 )*N ),
     $                   N )
            CALL <a name="SLASET.264"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, Q( IVT+( QSTART-1 )*N ),
     $                   N )
            CALL <a name="SLASDQ.266"></a><a href="slasdq.f.html#SLASDQ.1">SLASDQ</a>( <span class="string">'U'</span>, 0, N, N, N, 0, D, E,
     $                   Q( IVT+( QSTART-1 )*N ), N,
     $                   Q( IU+( QSTART-1 )*N ), N,
     $                   Q( IU+( QSTART-1 )*N ), N, WORK( WSTART ),
     $                   INFO )
         END IF
         GO TO 40
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( ICOMPQ.EQ.2 ) THEN
         CALL <a name="SLASET.276"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, U, LDU )
         CALL <a name="SLASET.277"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, VT, LDVT )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Scale.
</span><span class="comment">*</span><span class="comment">
</span>      ORGNRM = <a name="SLANST.282"></a><a href="slanst.f.html#SLANST.1">SLANST</a>( <span class="string">'M'</span>, N, D, E )
      IF( ORGNRM.EQ.ZERO )