clals0.f.html

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

</span><span class="comment">*</span><span class="comment">         Contain the components of the deflation-adjusted updating row
</span><span class="comment">*</span><span class="comment">         vector.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  K      (input) INTEGER
</span><span class="comment">*</span><span class="comment">         Contains the dimension of the non-deflated matrix,
</span><span class="comment">*</span><span class="comment">         This is the order of the related secular equation. 1 &lt;= K &lt;=N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  C      (input) REAL
</span><span class="comment">*</span><span class="comment">         C contains garbage if SQRE =0 and the C-value of a Givens
</span><span class="comment">*</span><span class="comment">         rotation related to the right null space if SQRE = 1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  S      (input) REAL
</span><span class="comment">*</span><span class="comment">         S contains garbage if SQRE =0 and the S-value of a Givens
</span><span class="comment">*</span><span class="comment">         rotation related to the right null space if SQRE = 1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RWORK  (workspace) REAL array, dimension
</span><span class="comment">*</span><span class="comment">         ( K*(1+NRHS) + 2*NRHS )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO   (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0:  successful exit.
</span><span class="comment">*</span><span class="comment">          &lt; 0:  if INFO = -i, the i-th argument had an illegal value.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Based on contributions by
</span><span class="comment">*</span><span class="comment">     Ming Gu and Ren-Cang Li, Computer Science Division, University of
</span><span class="comment">*</span><span class="comment">       California at Berkeley, USA
</span><span class="comment">*</span><span class="comment">     Osni Marques, LBNL/NERSC, USA
</span><span class="comment">*</span><span class="comment">
</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               ONE, ZERO, NEGONE
      PARAMETER          ( ONE = 1.0E0, ZERO = 0.0E0, NEGONE = -1.0E0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      INTEGER            I, J, JCOL, JROW, M, N, NLP1
      REAL               DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           CCOPY, <a name="CLACPY.176"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>, <a name="CLASCL.176"></a><a href="clascl.f.html#CLASCL.1">CLASCL</a>, CSROT, CSSCAL, SGEMV,
     $                   <a name="XERBLA.177"></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">     .. External Functions ..
</span>      REAL               <a name="SLAMC3.180"></a><a href="slamch.f.html#SLAMC3.574">SLAMC3</a>, SNRM2
      EXTERNAL           <a name="SLAMC3.181"></a><a href="slamch.f.html#SLAMC3.574">SLAMC3</a>, SNRM2
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          AIMAG, CMPLX, MAX, REAL
<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.220"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CLALS0.220"></a><a href="clals0.f.html#CLALS0.1">CLALS0</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 CSROT( 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 CCOPY( NRHS, B( NLP1, 1 ), LDB, BX( 1, 1 ), LDBX )
         DO 20 I = 2, N
            CALL CCOPY( 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 CCOPY( NRHS, BX, LDBX, B, LDB )
            IF( Z( 1 ).LT.ZERO ) THEN
               CALL CSSCAL( NRHS, NEGONE, B, LDB )
            END IF
         ELSE
            DO 100 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
                  RWORK( J ) = ZERO
               ELSE
                  RWORK( 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
                     RWORK( I ) = ZERO
                  ELSE
                     RWORK( I ) = POLES( I, 2 )*Z( I ) /
     $                            ( <a name="SLAMC3.276"></a><a href="slamch.f.html#SLAMC3.574">SLAMC3</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
                     RWORK( I ) = ZERO
                  ELSE
                     RWORK( I ) = POLES( I, 2 )*Z( I ) /
     $                            ( <a name="SLAMC3.286"></a><a href="slamch.f.html#SLAMC3.574">SLAMC3</a>( POLES( I, 2 ), DSIGJP )+
     $                            DIFRJ ) / ( POLES( I, 2 )+DJ )