dlals0.f.html

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

<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>      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
      END IF
<span class="comment">*</span><span class="comment">
</span>      N = NL + NR + 1
<span class="comment">*</span><span class="comment">
</span>      IF( NRHS.LT.1 ) THEN
         INFO = -5
      ELSE IF( LDB.LT.N ) THEN
         INFO = -7
      ELSE IF( LDBX.LT.N ) THEN
         INFO = -9
      ELSE IF( GIVPTR.LT.0 ) THEN
         INFO = -11
      ELSE IF( LDGCOL.LT.N ) THEN
         INFO = -13
      ELSE IF( LDGNUM.LT.N ) THEN
         INFO = -15
      ELSE IF( K.LT.1 ) THEN
         INFO = -20
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.218"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DLALS0.218"></a><a href="dlals0.f.html#DLALS0.1">DLALS0</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      M = N + SQRE
      NLP1 = NL + 1
<span class="comment">*</span><span class="comment">
</span>      IF( ICOMPQ.EQ.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Apply back orthogonal transformations from the left.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Step (1L): apply back the Givens rotations performed.
</span><span class="comment">*</span><span class="comment">
</span>         DO 10 I = 1, GIVPTR
            CALL DROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
     $                 B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
     $                 GIVNUM( I, 1 ) )
   10    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Step (2L): permute rows of B.
</span><span class="comment">*</span><span class="comment">
</span>         CALL DCOPY( NRHS, B( NLP1, 1 ), LDB, BX( 1, 1 ), LDBX )
         DO 20 I = 2, N
            CALL DCOPY( NRHS, B( PERM( I ), 1 ), LDB, BX( I, 1 ), LDBX )
   20    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Step (3L): apply the inverse of the left singular vector
</span><span class="comment">*</span><span class="comment">        matrix to BX.
</span><span class="comment">*</span><span class="comment">
</span>         IF( K.EQ.1 ) THEN
            CALL DCOPY( NRHS, BX, LDBX, B, LDB )
            IF( Z( 1 ).LT.ZERO ) THEN
               CALL DSCAL( NRHS, NEGONE, B, LDB )
            END IF
         ELSE
            DO 50 J = 1, K
               DIFLJ = DIFL( J )
               DJ = POLES( J, 1 )
               DSIGJ = -POLES( J, 2 )
               IF( J.LT.K ) THEN
                  DIFRJ = -DIFR( J, 1 )
                  DSIGJP = -POLES( J+1, 2 )
               END IF
               IF( ( Z( J ).EQ.ZERO ) .OR. ( POLES( J, 2 ).EQ.ZERO ) )
     $              THEN
                  WORK( J ) = ZERO
               ELSE
                  WORK( J ) = -POLES( J, 2 )*Z( J ) / DIFLJ /
     $                        ( POLES( J, 2 )+DJ )
               END IF
               DO 30 I = 1, J - 1
                  IF( ( Z( I ).EQ.ZERO ) .OR.
     $                ( POLES( I, 2 ).EQ.ZERO ) ) THEN
                     WORK( I ) = ZERO
                  ELSE
                     WORK( I ) = POLES( I, 2 )*Z( I ) /
     $                           ( <a name="DLAMC3.274"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( POLES( I, 2 ), DSIGJ )-
     $                           DIFLJ ) / ( POLES( I, 2 )+DJ )
                  END IF
   30          CONTINUE
               DO 40 I = J + 1, K
                  IF( ( Z( I ).EQ.ZERO ) .OR.
     $                ( POLES( I, 2 ).EQ.ZERO ) ) THEN
                     WORK( I ) = ZERO
                  ELSE
                     WORK( I ) = POLES( I, 2 )*Z( I ) /
     $                           ( <a name="DLAMC3.284"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( POLES( I, 2 ), DSIGJP )+
     $                           DIFRJ ) / ( POLES( I, 2 )+DJ )
                  END IF
   40          CONTINUE
               WORK( 1 ) = NEGONE
               TEMP = DNRM2( K, WORK, 1 )
               CALL DGEMV( <span class="string">'T'</span>, K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
     $                     B( J, 1 ), LDB )
               CALL <a name="DLASCL.292"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, TEMP, ONE, 1, NRHS, B( J, 1 ),
     $                      LDB, INFO )
   50       CONTINUE
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Move the deflated rows of BX to B also.
</span><span class="comment">*</span><span class="comment">
</span>         IF( K.LT.MAX( M, N ) )
     $      CALL <a name="DLACPY.300"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'A'</span>, N-K, NRHS, BX( K+1, 1 ), LDBX,
     $                   B( K+1, 1 ), LDB )
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Apply back the right orthogonal transformations.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Step (1R): apply back the new right singular vector matrix
</span><span class="comment">*</span><span class="comment">        to B.
</span><span class="comment">*</span><span class="comment">
</span>         IF( K.EQ.1 ) THEN
            CALL DCOPY( NRHS, B, LDB, BX, LDBX )
         ELSE
            DO 80 J = 1, K
               DSIGJ = POLES( J, 2 )
               IF( Z( J ).EQ.ZERO ) THEN
                  WORK( J ) = ZERO
               ELSE
                  WORK( J ) = -Z( J ) / DIFL( J ) /
     $                        ( DSIGJ+POLES( J, 1 ) ) / DIFR( J, 2 )
               END IF
               DO 60 I = 1, J - 1
                  IF( Z( J ).EQ.ZERO ) THEN
                     WORK( I ) = ZERO
                  ELSE
                     WORK( I ) = Z( J ) / ( <a name="DLAMC3.324"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( DSIGJ, -POLES( I+1,
     $                           2 ) )-DIFR( I, 1 ) ) /
     $                           ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
                  END IF
   60          CONTINUE
               DO 70 I = J + 1, K
                  IF( Z( J ).EQ.ZERO ) THEN
                     WORK( I ) = ZERO
                  ELSE
                     WORK( I ) = Z( J ) / ( <a name="DLAMC3.333"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( DSIGJ, -POLES( I,
     $                           2 ) )-DIFL( I ) ) /
     $                           ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
                  END IF
   70          CONTINUE
               CALL DGEMV( <span class="string">'T'</span>, K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
     $                     BX( J, 1 ), LDBX )
   80       CONTINUE
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Step (2R): if SQRE = 1, apply back the rotation that is
</span><span class="comment">*</span><span class="comment">        related to the right null space of the subproblem.
</span><span class="comment">*</span><span class="comment">
</span>         IF( SQRE.EQ.1 ) THEN
            CALL DCOPY( NRHS, B( M, 1 ), LDB, BX( M, 1 ), LDBX )
            CALL DROT( NRHS, BX( 1, 1 ), LDBX, BX( M, 1 ), LDBX, C, S )
         END IF
         IF( K.LT.MAX( M, N ) )
     $      CALL <a name="DLACPY.351"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'A'</span>, N-K, NRHS, B( K+1, 1 ), LDB, BX( K+1, 1 ),
     $                   LDBX )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Step (3R): permute rows of B.
</span><span class="comment">*</span><span class="comment">
</span>         CALL DCOPY( NRHS, BX( 1, 1 ), LDBX, B( NLP1, 1 ), LDB )
         IF( SQRE.EQ.1 ) THEN
            CALL DCOPY( NRHS, BX( M, 1 ), LDBX, B( M, 1 ), LDB )
         END IF
         DO 90 I = 2, N
            CALL DCOPY( NRHS, BX( I, 1 ), LDBX, B( PERM( I ), 1 ), LDB )
   90    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Step (4R): apply back the Givens rotations performed.
</span><span class="comment">*</span><span class="comment">
</span>         DO 100 I = GIVPTR, 1, -1
            CALL DROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
     $                 B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
     $                 -GIVNUM( I, 1 ) )
  100    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="DLALS0.375"></a><a href="dlals0.f.html#DLALS0.1">DLALS0</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>