clahr2.f.html

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

<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      INTEGER            I
      COMPLEX            EI
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           CAXPY, CCOPY, CGEMM, CGEMV, <a name="CLACPY.121"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>,
     $                   <a name="CLARFG.122"></a><a href="clarfg.f.html#CLARFG.1">CLARFG</a>, CSCAL, CTRMM, CTRMV, <a name="CLACGV.122"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          MIN
<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">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.LE.1 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span>      DO 10 I = 1, NB
         IF( I.GT.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Update A(K+1:N,I)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Update I-th column of A - Y * V'
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="CLACGV.141"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( I-1, A( K+I-1, 1 ), LDA ) 
            CALL CGEMV( <span class="string">'NO TRANSPOSE'</span>, N-K, I-1, -ONE, Y(K+1,1), LDY,
     $                  A( K+I-1, 1 ), LDA, ONE, A( K+1, I ), 1 )
            CALL <a name="CLACGV.144"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( I-1, A( K+I-1, 1 ), LDA ) 
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Apply I - V * T' * V' to this column (call it b) from the
</span><span class="comment">*</span><span class="comment">           left, using the last column of T as workspace
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Let  V = ( V1 )   and   b = ( b1 )   (first I-1 rows)
</span><span class="comment">*</span><span class="comment">                    ( V2 )             ( b2 )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           where V1 is unit lower triangular
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           w := V1' * b1
</span><span class="comment">*</span><span class="comment">
</span>            CALL CCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 )
            CALL CTRMV( <span class="string">'Lower'</span>, <span class="string">'Conjugate transpose'</span>, <span class="string">'UNIT'</span>, 
     $                  I-1, A( K+1, 1 ),
     $                  LDA, T( 1, NB ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           w := w + V2'*b2
</span><span class="comment">*</span><span class="comment">
</span>            CALL CGEMV( <span class="string">'Conjugate transpose'</span>, N-K-I+1, I-1, 
     $                  ONE, A( K+I, 1 ),
     $                  LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           w := T'*w
</span><span class="comment">*</span><span class="comment">
</span>            CALL CTRMV( <span class="string">'Upper'</span>, <span class="string">'Conjugate transpose'</span>, <span class="string">'NON-UNIT'</span>, 
     $                  I-1, T, LDT,
     $                  T( 1, NB ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           b2 := b2 - V2*w
</span><span class="comment">*</span><span class="comment">
</span>            CALL CGEMV( <span class="string">'NO TRANSPOSE'</span>, N-K-I+1, I-1, -ONE, 
     $                  A( K+I, 1 ),
     $                  LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           b1 := b1 - V1*w
</span><span class="comment">*</span><span class="comment">
</span>            CALL CTRMV( <span class="string">'Lower'</span>, <span class="string">'NO TRANSPOSE'</span>, 
     $                  <span class="string">'UNIT'</span>, I-1,
     $                  A( K+1, 1 ), LDA, T( 1, NB ), 1 )
            CALL CAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 )
<span class="comment">*</span><span class="comment">
</span>            A( K+I-1, I-1 ) = EI
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Generate the elementary reflector H(I) to annihilate
</span><span class="comment">*</span><span class="comment">        A(K+I+1:N,I)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CLARFG.192"></a><a href="clarfg.f.html#CLARFG.1">CLARFG</a>( N-K-I+1, A( K+I, I ), A( MIN( K+I+1, N ), I ), 1,
     $                TAU( I ) )
         EI = A( K+I, I )
         A( K+I, I ) = ONE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute  Y(K+1:N,I)
</span><span class="comment">*</span><span class="comment">
</span>         CALL CGEMV( <span class="string">'NO TRANSPOSE'</span>, N-K, N-K-I+1, 
     $               ONE, A( K+1, I+1 ),
     $               LDA, A( K+I, I ), 1, ZERO, Y( K+1, I ), 1 )
         CALL CGEMV( <span class="string">'Conjugate transpose'</span>, N-K-I+1, I-1, 
     $               ONE, A( K+I, 1 ), LDA,
     $               A( K+I, I ), 1, ZERO, T( 1, I ), 1 )
         CALL CGEMV( <span class="string">'NO TRANSPOSE'</span>, N-K, I-1, -ONE, 
     $               Y( K+1, 1 ), LDY,
     $               T( 1, I ), 1, ONE, Y( K+1, I ), 1 )
         CALL CSCAL( N-K, TAU( I ), Y( K+1, I ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute T(1:I,I)
</span><span class="comment">*</span><span class="comment">
</span>         CALL CSCAL( I-1, -TAU( I ), T( 1, I ), 1 )
         CALL CTRMV( <span class="string">'Upper'</span>, <span class="string">'No Transpose'</span>, <span class="string">'NON-UNIT'</span>, 
     $               I-1, T, LDT,
     $               T( 1, I ), 1 )
         T( I, I ) = TAU( I )
<span class="comment">*</span><span class="comment">
</span>   10 CONTINUE
      A( K+NB, NB ) = EI
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute Y(1:K,1:NB)
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="CLACPY.223"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'ALL'</span>, K, NB, A( 1, 2 ), LDA, Y, LDY )
      CALL CTRMM( <span class="string">'RIGHT'</span>, <span class="string">'Lower'</span>, <span class="string">'NO TRANSPOSE'</span>, 
     $            <span class="string">'UNIT'</span>, K, NB,
     $            ONE, A( K+1, 1 ), LDA, Y, LDY )
      IF( N.GT.K+NB )
     $   CALL CGEMM( <span class="string">'NO TRANSPOSE'</span>, <span class="string">'NO TRANSPOSE'</span>, K, 
     $               NB, N-K-NB, ONE,
     $               A( 1, 2+NB ), LDA, A( K+1+NB, 1 ), LDA, ONE, Y,
     $               LDY )
      CALL CTRMM( <span class="string">'RIGHT'</span>, <span class="string">'Upper'</span>, <span class="string">'NO TRANSPOSE'</span>, 
     $            <span class="string">'NON-UNIT'</span>, K, NB,
     $            ONE, T, LDT, Y, LDY )
<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="CLAHR2.238"></a><a href="clahr2.f.html#CLAHR2.1">CLAHR2</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>