slalsd.f.html

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

</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( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( NRHS.LT.1 ) THEN
         INFO = -4
      ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN
         INFO = -8
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.147"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SLALSD.147"></a><a href="slalsd.f.html#SLALSD.1">SLALSD</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      EPS = <a name="SLAMCH.151"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'Epsilon'</span> )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set up the tolerance.
</span><span class="comment">*</span><span class="comment">
</span>      IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
         RCND = EPS
      ELSE
         RCND = RCOND
      END IF
<span class="comment">*</span><span class="comment">
</span>      RANK = 0
<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 ) THEN
         RETURN
      ELSE IF( N.EQ.1 ) THEN
         IF( D( 1 ).EQ.ZERO ) THEN
            CALL <a name="SLASET.169"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'A'</span>, 1, NRHS, ZERO, ZERO, B, LDB )
         ELSE
            RANK = 1
            CALL <a name="SLASCL.172"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO )
            D( 1 ) = ABS( D( 1 ) )
         END IF
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Rotate the matrix if it is lower bidiagonal.
</span><span class="comment">*</span><span class="comment">
</span>      IF( UPLO.EQ.<span class="string">'L'</span> ) THEN
         DO 10 I = 1, N - 1
            CALL <a name="SLARTG.182"></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( NRHS.EQ.1 ) THEN
               CALL SROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN )
            ELSE
               WORK( I*2-1 ) = CS
               WORK( I*2 ) = SN
            END IF
   10    CONTINUE
         IF( NRHS.GT.1 ) THEN
            DO 30 I = 1, NRHS
               DO 20 J = 1, N - 1
                  CS = WORK( J*2-1 )
                  SN = WORK( J*2 )
                  CALL SROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN )
   20          CONTINUE
   30       CONTINUE
         END IF
      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>      NM1 = N - 1
      ORGNRM = <a name="SLANST.207"></a><a href="slanst.f.html#SLANST.1">SLANST</a>( <span class="string">'M'</span>, N, D, E )
      IF( ORGNRM.EQ.ZERO ) THEN
         CALL <a name="SLASET.209"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'A'</span>, N, NRHS, ZERO, ZERO, B, LDB )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      CALL <a name="SLASCL.213"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
      CALL <a name="SLASCL.214"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO )
<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
         NWORK = 1 + N*N
         CALL <a name="SLASET.221"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, WORK, N )
         CALL <a name="SLASDQ.222"></a><a href="slasdq.f.html#SLASDQ.1">SLASDQ</a>( <span class="string">'U'</span>, 0, N, N, 0, NRHS, D, E, WORK, N, WORK, N, B,
     $                LDB, WORK( NWORK ), INFO )
         IF( INFO.NE.0 ) THEN
            RETURN
         END IF
         TOL = RCND*ABS( D( ISAMAX( N, D, 1 ) ) )
         DO 40 I = 1, N
            IF( D( I ).LE.TOL ) THEN
               CALL <a name="SLASET.230"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'A'</span>, 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
            ELSE
               CALL <a name="SLASCL.232"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ),
     $                      LDB, INFO )
               RANK = RANK + 1
            END IF
   40    CONTINUE
         CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
     $               WORK( NWORK ), N )
         CALL <a name="SLACPY.239"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'A'</span>, N, NRHS, WORK( NWORK ), N, B, LDB )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Unscale.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SLASCL.243"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
         CALL <a name="SLASRT.244"></a><a href="slasrt.f.html#SLASRT.1">SLASRT</a>( <span class="string">'D'</span>, N, D, INFO )
         CALL <a name="SLASCL.245"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
<span class="comment">*</span><span class="comment">
</span>         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Book-keeping and setting up some constants.
</span><span class="comment">*</span><span class="comment">
</span>      NLVL = INT( LOG( REAL( N ) / REAL( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
<span class="comment">*</span><span class="comment">
</span>      SMLSZP = SMLSIZ + 1
<span class="comment">*</span><span class="comment">
</span>      U = 1
      VT = 1 + SMLSIZ*N
      DIFL = VT + SMLSZP*N
      DIFR = DIFL + NLVL*N
      Z = DIFR + NLVL*N*2
      C = Z + NLVL*N
      S = C + N
      POLES = S + N
      GIVNUM = POLES + 2*NLVL*N
      BX = GIVNUM + 2*NLVL*N
      NWORK = BX + N*NRHS
<span class="comment">*</span><span class="comment">
</span>      SIZEI = 1 + N
      K = SIZEI + N
      GIVPTR = K + N
      PERM = GIVPTR + N
      GIVCOL = PERM + NLVL*N
      IWK = GIVCOL + NLVL*N*2
<span class="comment">*</span><span class="comment">
</span>      ST = 1
      SQRE = 0
      ICMPQ1 = 1
      ICMPQ2 = 0
      NSUB = 0
<span class="comment">*</span><span class="comment">
</span>      DO 50 I = 1, N
         IF( ABS( D( I ) ).LT.EPS ) THEN
            D( I ) = SIGN( EPS, D( I ) )
         END IF
   50 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      DO 60 I = 1, NM1