clarft.f.html

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

      PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ),
     $                   ZERO = ( 0.0E+0, 0.0E+0 ) )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      INTEGER            I, J
      COMPLEX            VII
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           CGEMV, <a name="CLACGV.115"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>, CTRMV
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.118"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      EXTERNAL           <a name="LSAME.119"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
<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.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span>      IF( <a name="LSAME.128"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( DIRECT, <span class="string">'F'</span> ) ) THEN
         DO 20 I = 1, K
            IF( TAU( I ).EQ.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              H(i)  =  I
</span><span class="comment">*</span><span class="comment">
</span>               DO 10 J = 1, I
                  T( J, I ) = ZERO
   10          CONTINUE
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              general case
</span><span class="comment">*</span><span class="comment">
</span>               VII = V( I, I )
               V( I, I ) = ONE
               IF( <a name="LSAME.143"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( STOREV, <span class="string">'C'</span> ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 T(1:i-1,i) := - tau(i) * V(i:n,1:i-1)' * V(i:n,i)
</span><span class="comment">*</span><span class="comment">
</span>                  CALL CGEMV( <span class="string">'Conjugate transpose'</span>, N-I+1, I-1,
     $                        -TAU( I ), V( I, 1 ), LDV, V( I, I ), 1,
     $                        ZERO, T( 1, I ), 1 )
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 T(1:i-1,i) := - tau(i) * V(1:i-1,i:n) * V(i,i:n)'
</span><span class="comment">*</span><span class="comment">
</span>                  IF( I.LT.N )
     $               CALL <a name="CLACGV.155"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( N-I, V( I, I+1 ), LDV )
                  CALL CGEMV( <span class="string">'No transpose'</span>, I-1, N-I+1, -TAU( I ),
     $                        V( 1, I ), LDV, V( I, I ), LDV, ZERO,
     $                        T( 1, I ), 1 )
                  IF( I.LT.N )
     $               CALL <a name="CLACGV.160"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( N-I, V( I, I+1 ), LDV )
               END IF
               V( I, I ) = VII
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
</span><span class="comment">*</span><span class="comment">
</span>               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 )
            END IF
   20    CONTINUE
      ELSE
         DO 40 I = K, 1, -1
            IF( TAU( I ).EQ.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              H(i)  =  I
</span><span class="comment">*</span><span class="comment">
</span>               DO 30 J = I, K
                  T( J, I ) = ZERO
   30          CONTINUE
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              general case
</span><span class="comment">*</span><span class="comment">
</span>               IF( I.LT.K ) THEN
                  IF( <a name="LSAME.185"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( STOREV, <span class="string">'C'</span> ) ) THEN
                     VII = V( N-K+I, I )
                     V( N-K+I, I ) = ONE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    T(i+1:k,i) :=
</span><span class="comment">*</span><span class="comment">                            - tau(i) * V(1:n-k+i,i+1:k)' * V(1:n-k+i,i)
</span><span class="comment">*</span><span class="comment">
</span>                     CALL CGEMV( <span class="string">'Conjugate transpose'</span>, N-K+I, K-I,
     $                           -TAU( I ), V( 1, I+1 ), LDV, V( 1, I ),
     $                           1, ZERO, T( I+1, I ), 1 )
                     V( N-K+I, I ) = VII
                  ELSE
                     VII = V( I, N-K+I )
                     V( I, N-K+I ) = ONE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    T(i+1:k,i) :=
</span><span class="comment">*</span><span class="comment">                            - tau(i) * V(i+1:k,1:n-k+i) * V(i,1:n-k+i)'
</span><span class="comment">*</span><span class="comment">
</span>                     CALL <a name="CLACGV.203"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( N-K+I-1, V( I, 1 ), LDV )
                     CALL CGEMV( <span class="string">'No transpose'</span>, K-I, N-K+I, -TAU( I ),
     $                           V( I+1, 1 ), LDV, V( I, 1 ), LDV, ZERO,
     $                           T( I+1, I ), 1 )
                     CALL <a name="CLACGV.207"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( N-K+I-1, V( I, 1 ), LDV )
                     V( I, N-K+I ) = VII
                  END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i)
</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">'Non-unit'</span>, K-I,
     $                        T( I+1, I+1 ), LDT, T( I+1, I ), 1 )
               END IF
               T( I, I ) = TAU( I )
            END IF
   40    CONTINUE
      END IF
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="CLARFT.222"></a><a href="clarft.f.html#CLARFT.1">CLARFT</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>