dbdsdc.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 454 行 · 第 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>      DOUBLE PRECISION   ZERO, ONE, TWO
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+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
      DOUBLE PRECISION   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.148"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      INTEGER            <a name="ILAENV.149"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
      DOUBLE PRECISION   <a name="DLAMCH.150"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANST.150"></a><a href="dlanst.f.html#DLANST.1">DLANST</a>
      EXTERNAL           <a name="LSAME.151"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="ILAENV.151"></a><a href="hfy-index.html#ILAENV">ILAENV</a>, <a name="DLAMCH.151"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANST.151"></a><a href="dlanst.f.html#DLANST.1">DLANST</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           DCOPY, <a name="DLARTG.154"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>, <a name="DLASCL.154"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>, <a name="DLASD0.154"></a><a href="dlasd0.f.html#DLASD0.1">DLASD0</a>, <a name="DLASDA.154"></a><a href="dlasda.f.html#DLASDA.1">DLASDA</a>, <a name="DLASDQ.154"></a><a href="dlasdq.f.html#DLASDQ.1">DLASDQ</a>,
     $                   <a name="DLASET.155"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>, <a name="DLASR.155"></a><a href="dlasr.f.html#DLASR.1">DLASR</a>, DSWAP, <a name="XERBLA.155"></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, DBLE, 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.167"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'U'</span> ) )
     $   IUPLO = 1
      IF( <a name="LSAME.169"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'L'</span> ) )
     $   IUPLO = 2
      IF( <a name="LSAME.171"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPQ, <span class="string">'N'</span> ) ) THEN
         ICOMPQ = 0
      ELSE IF( <a name="LSAME.173"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPQ, <span class="string">'P'</span> ) ) THEN
         ICOMPQ = 1
      ELSE IF( <a name="LSAME.175"></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.194"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DBDSDC.194"></a><a href="dbdsdc.f.html#DBDSDC.1">DBDSDC</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.202"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 9, <span class="string">'<a name="DBDSDC.202"></a><a href="dbdsdc.f.html#DBDSDC.1">DBDSDC</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 DCOPY( N, D, 1, Q( 1 ), 1 )
         CALL DCOPY( 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="DLARTG.229"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</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="DLASDQ.243"></a><a href="dlasdq.f.html#DLASDQ.1">DLASDQ</a> to compute the singular values.
</span><span class="comment">*</span><span class="comment">
</span>      IF( ICOMPQ.EQ.0 ) THEN
         CALL <a name="DLASDQ.246"></a><a href="dlasdq.f.html#DLASDQ.1">DLASDQ</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="DLASET.256"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, U, LDU )
            CALL <a name="DLASET.257"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, VT, LDVT )
            CALL <a name="DLASDQ.258"></a><a href="dlasdq.f.html#DLASDQ.1">DLASDQ</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="DLASET.263"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, Q( IU+( QSTART-1 )*N ),
     $                   N )
            CALL <a name="DLASET.265"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, Q( IVT+( QSTART-1 )*N ),
     $                   N )
            CALL <a name="DLASDQ.267"></a><a href="dlasdq.f.html#DLASDQ.1">DLASDQ</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="DLASET.277"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, U, LDU )
         CALL <a name="DLASET.278"></a><a href="dlaset.f.html#DLASET.1">DLASET</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="DLANST.283"></a><a href="dlanst.f.html#DLANST.1">DLANST</a>( <span class="string">'M'</span>, N, D, E )
      IF( ORGNRM.EQ.ZERO )
     $   RETURN