clals0.f.html

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

                  END IF
   40          CONTINUE
               RWORK( 1 ) = NEGONE
               TEMP = SNRM2( K, RWORK, 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Since B and BX are complex, the following call to SGEMV
</span><span class="comment">*</span><span class="comment">              is performed in two steps (real and imaginary parts).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              CALL SGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
</span><span class="comment">*</span><span class="comment">    $                     B( J, 1 ), LDB )
</span><span class="comment">*</span><span class="comment">
</span>               I = K + NRHS*2
               DO 60 JCOL = 1, NRHS
                  DO 50 JROW = 1, K
                     I = I + 1
                     RWORK( I ) = REAL( BX( JROW, JCOL ) )
   50             CONTINUE
   60          CONTINUE
               CALL SGEMV( <span class="string">'T'</span>, K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
     $                     RWORK( 1 ), 1, ZERO, RWORK( 1+K ), 1 )
               I = K + NRHS*2
               DO 80 JCOL = 1, NRHS
                  DO 70 JROW = 1, K
                     I = I + 1
                     RWORK( I ) = AIMAG( BX( JROW, JCOL ) )
   70             CONTINUE
   80          CONTINUE
               CALL SGEMV( <span class="string">'T'</span>, K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
     $                     RWORK( 1 ), 1, ZERO, RWORK( 1+K+NRHS ), 1 )
               DO 90 JCOL = 1, NRHS
                  B( J, JCOL ) = CMPLX( RWORK( JCOL+K ),
     $                           RWORK( JCOL+K+NRHS ) )
   90          CONTINUE
               CALL <a name="CLASCL.321"></a><a href="clascl.f.html#CLASCL.1">CLASCL</a>( <span class="string">'G'</span>, 0, 0, TEMP, ONE, 1, NRHS, B( J, 1 ),
     $                      LDB, INFO )
  100       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="CLACPY.329"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</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 CCOPY( NRHS, B, LDB, BX, LDBX )
         ELSE
            DO 180 J = 1, K
               DSIGJ = POLES( J, 2 )
               IF( Z( J ).EQ.ZERO ) THEN
                  RWORK( J ) = ZERO
               ELSE
                  RWORK( J ) = -Z( J ) / DIFL( J ) /
     $                         ( DSIGJ+POLES( J, 1 ) ) / DIFR( J, 2 )
               END IF
               DO 110 I = 1, J - 1
                  IF( Z( J ).EQ.ZERO ) THEN
                     RWORK( I ) = ZERO
                  ELSE
                     RWORK( I ) = Z( J ) / ( <a name="SLAMC3.353"></a><a href="slamch.f.html#SLAMC3.574">SLAMC3</a>( DSIGJ, -POLES( I+1,
     $                            2 ) )-DIFR( I, 1 ) ) /
     $                            ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
                  END IF
  110          CONTINUE
               DO 120 I = J + 1, K
                  IF( Z( J ).EQ.ZERO ) THEN
                     RWORK( I ) = ZERO
                  ELSE
                     RWORK( I ) = Z( J ) / ( <a name="SLAMC3.362"></a><a href="slamch.f.html#SLAMC3.574">SLAMC3</a>( DSIGJ, -POLES( I,
     $                            2 ) )-DIFL( I ) ) /
     $                            ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
                  END IF
  120          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Since B and BX are complex, the following call to SGEMV
</span><span class="comment">*</span><span class="comment">              is performed in two steps (real and imaginary parts).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              CALL SGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
</span><span class="comment">*</span><span class="comment">    $                     BX( J, 1 ), LDBX )
</span><span class="comment">*</span><span class="comment">
</span>               I = K + NRHS*2
               DO 140 JCOL = 1, NRHS
                  DO 130 JROW = 1, K
                     I = I + 1
                     RWORK( I ) = REAL( B( JROW, JCOL ) )
  130             CONTINUE
  140          CONTINUE
               CALL SGEMV( <span class="string">'T'</span>, K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
     $                     RWORK( 1 ), 1, ZERO, RWORK( 1+K ), 1 )
               I = K + NRHS*2
               DO 160 JCOL = 1, NRHS
                  DO 150 JROW = 1, K
                     I = I + 1
                     RWORK( I ) = AIMAG( B( JROW, JCOL ) )
  150             CONTINUE
  160          CONTINUE
               CALL SGEMV( <span class="string">'T'</span>, K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
     $                     RWORK( 1 ), 1, ZERO, RWORK( 1+K+NRHS ), 1 )
               DO 170 JCOL = 1, NRHS
                  BX( J, JCOL ) = CMPLX( RWORK( JCOL+K ),
     $                            RWORK( JCOL+K+NRHS ) )
  170          CONTINUE
  180       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 CCOPY( NRHS, B( M, 1 ), LDB, BX( M, 1 ), LDBX )
            CALL CSROT( NRHS, BX( 1, 1 ), LDBX, BX( M, 1 ), LDBX, C, S )
         END IF
         IF( K.LT.MAX( M, N ) )
     $      CALL <a name="CLACPY.407"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</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 CCOPY( NRHS, BX( 1, 1 ), LDBX, B( NLP1, 1 ), LDB )
         IF( SQRE.EQ.1 ) THEN
            CALL CCOPY( NRHS, BX( M, 1 ), LDBX, B( M, 1 ), LDB )
         END IF
         DO 190 I = 2, N
            CALL CCOPY( NRHS, BX( I, 1 ), LDBX, B( PERM( I ), 1 ), LDB )
  190    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 200 I = GIVPTR, 1, -1
            CALL CSROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
     $                  B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
     $                  -GIVNUM( I, 1 ) )
  200    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="CLALS0.431"></a><a href="clals0.f.html#CLALS0.1">CLALS0</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>