dlabrd.f.html

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

<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      INTEGER            I
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           DGEMV, <a name="DLARFG.145"></a><a href="dlarfg.f.html#DLARFG.1">DLARFG</a>, DSCAL
<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( M.LE.0 .OR. N.LE.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span>      IF( M.GE.N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Reduce to upper bidiagonal form
</span><span class="comment">*</span><span class="comment">
</span>         DO 10 I = 1, NB
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Update A(i:m,i)
</span><span class="comment">*</span><span class="comment">
</span>            CALL DGEMV( <span class="string">'No transpose'</span>, M-I+1, I-1, -ONE, A( I, 1 ),
     $                  LDA, Y( I, 1 ), LDY, ONE, A( I, I ), 1 )
            CALL DGEMV( <span class="string">'No transpose'</span>, M-I+1, I-1, -ONE, X( I, 1 ),
     $                  LDX, A( 1, I ), 1, ONE, A( I, I ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Generate reflection Q(i) to annihilate A(i+1:m,i)
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DLARFG.172"></a><a href="dlarfg.f.html#DLARFG.1">DLARFG</a>( M-I+1, A( I, I ), A( MIN( I+1, M ), I ), 1,
     $                   TAUQ( I ) )
            D( I ) = A( I, I )
            IF( I.LT.N ) THEN
               A( I, I ) = ONE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute Y(i+1:n,i)
</span><span class="comment">*</span><span class="comment">
</span>               CALL DGEMV( <span class="string">'Transpose'</span>, M-I+1, N-I, ONE, A( I, I+1 ),
     $                     LDA, A( I, I ), 1, ZERO, Y( I+1, I ), 1 )
               CALL DGEMV( <span class="string">'Transpose'</span>, M-I+1, I-1, ONE, A( I, 1 ), LDA,
     $                     A( I, I ), 1, ZERO, Y( 1, I ), 1 )
               CALL DGEMV( <span class="string">'No transpose'</span>, N-I, I-1, -ONE, Y( I+1, 1 ),
     $                     LDY, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
               CALL DGEMV( <span class="string">'Transpose'</span>, M-I+1, I-1, ONE, X( I, 1 ), LDX,
     $                     A( I, I ), 1, ZERO, Y( 1, I ), 1 )
               CALL DGEMV( <span class="string">'Transpose'</span>, I-1, N-I, -ONE, A( 1, I+1 ),
     $                     LDA, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
               CALL DSCAL( N-I, TAUQ( I ), Y( I+1, I ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Update A(i,i+1:n)
</span><span class="comment">*</span><span class="comment">
</span>               CALL DGEMV( <span class="string">'No transpose'</span>, N-I, I, -ONE, Y( I+1, 1 ),
     $                     LDY, A( I, 1 ), LDA, ONE, A( I, I+1 ), LDA )
               CALL DGEMV( <span class="string">'Transpose'</span>, I-1, N-I, -ONE, A( 1, I+1 ),
     $                     LDA, X( I, 1 ), LDX, ONE, A( I, I+1 ), LDA )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Generate reflection P(i) to annihilate A(i,i+2:n)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DLARFG.201"></a><a href="dlarfg.f.html#DLARFG.1">DLARFG</a>( N-I, A( I, I+1 ), A( I, MIN( I+2, N ) ),
     $                      LDA, TAUP( I ) )
               E( I ) = A( I, I+1 )
               A( I, I+1 ) = ONE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute X(i+1:m,i)
</span><span class="comment">*</span><span class="comment">
</span>               CALL DGEMV( <span class="string">'No transpose'</span>, M-I, N-I, ONE, A( I+1, I+1 ),
     $                     LDA, A( I, I+1 ), LDA, ZERO, X( I+1, I ), 1 )
               CALL DGEMV( <span class="string">'Transpose'</span>, N-I, I, ONE, Y( I+1, 1 ), LDY,
     $                     A( I, I+1 ), LDA, ZERO, X( 1, I ), 1 )
               CALL DGEMV( <span class="string">'No transpose'</span>, M-I, I, -ONE, A( I+1, 1 ),
     $                     LDA, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
               CALL DGEMV( <span class="string">'No transpose'</span>, I-1, N-I, ONE, A( 1, I+1 ),
     $                     LDA, A( I, I+1 ), LDA, ZERO, X( 1, I ), 1 )
               CALL DGEMV( <span class="string">'No transpose'</span>, M-I, I-1, -ONE, X( I+1, 1 ),
     $                     LDX, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
               CALL DSCAL( M-I, TAUP( I ), X( I+1, I ), 1 )
            END IF
   10    CONTINUE
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Reduce to lower bidiagonal form
</span><span class="comment">*</span><span class="comment">
</span>         DO 20 I = 1, NB
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Update A(i,i:n)
</span><span class="comment">*</span><span class="comment">
</span>            CALL DGEMV( <span class="string">'No transpose'</span>, N-I+1, I-1, -ONE, Y( I, 1 ),
     $                  LDY, A( I, 1 ), LDA, ONE, A( I, I ), LDA )
            CALL DGEMV( <span class="string">'Transpose'</span>, I-1, N-I+1, -ONE, A( 1, I ), LDA,
     $                  X( I, 1 ), LDX, ONE, A( I, I ), LDA )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Generate reflection P(i) to annihilate A(i,i+1:n)
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DLARFG.236"></a><a href="dlarfg.f.html#DLARFG.1">DLARFG</a>( N-I+1, A( I, I ), A( I, MIN( I+1, N ) ), LDA,
     $                   TAUP( I ) )
            D( I ) = A( I, I )
            IF( I.LT.M ) THEN
               A( I, I ) = ONE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute X(i+1:m,i)
</span><span class="comment">*</span><span class="comment">
</span>               CALL DGEMV( <span class="string">'No transpose'</span>, M-I, N-I+1, ONE, A( I+1, I ),
     $                     LDA, A( I, I ), LDA, ZERO, X( I+1, I ), 1 )
               CALL DGEMV( <span class="string">'Transpose'</span>, N-I+1, I-1, ONE, Y( I, 1 ), LDY,
     $                     A( I, I ), LDA, ZERO, X( 1, I ), 1 )
               CALL DGEMV( <span class="string">'No transpose'</span>, M-I, I-1, -ONE, A( I+1, 1 ),
     $                     LDA, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
               CALL DGEMV( <span class="string">'No transpose'</span>, I-1, N-I+1, ONE, A( 1, I ),
     $                     LDA, A( I, I ), LDA, ZERO, X( 1, I ), 1 )
               CALL DGEMV( <span class="string">'No transpose'</span>, M-I, I-1, -ONE, X( I+1, 1 ),
     $                     LDX, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
               CALL DSCAL( M-I, TAUP( I ), X( I+1, I ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Update A(i+1:m,i)
</span><span class="comment">*</span><span class="comment">
</span>               CALL DGEMV( <span class="string">'No transpose'</span>, M-I, I-1, -ONE, A( I+1, 1 ),
     $                     LDA, Y( I, 1 ), LDY, ONE, A( I+1, I ), 1 )
               CALL DGEMV( <span class="string">'No transpose'</span>, M-I, I, -ONE, X( I+1, 1 ),
     $                     LDX, A( 1, I ), 1, ONE, A( I+1, I ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Generate reflection Q(i) to annihilate A(i+2:m,i)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DLARFG.265"></a><a href="dlarfg.f.html#DLARFG.1">DLARFG</a>( M-I, A( I+1, I ), A( MIN( I+2, M ), I ), 1,
     $                      TAUQ( I ) )
               E( I ) = A( I+1, I )
               A( I+1, I ) = ONE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute Y(i+1:n,i)
</span><span class="comment">*</span><span class="comment">
</span>               CALL DGEMV( <span class="string">'Transpose'</span>, M-I, N-I, ONE, A( I+1, I+1 ),
     $                     LDA, A( I+1, I ), 1, ZERO, Y( I+1, I ), 1 )
               CALL DGEMV( <span class="string">'Transpose'</span>, M-I, I-1, ONE, A( I+1, 1 ), LDA,
     $                     A( I+1, I ), 1, ZERO, Y( 1, I ), 1 )
               CALL DGEMV( <span class="string">'No transpose'</span>, N-I, I-1, -ONE, Y( I+1, 1 ),
     $                     LDY, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
               CALL DGEMV( <span class="string">'Transpose'</span>, M-I, I, ONE, X( I+1, 1 ), LDX,
     $                     A( I+1, I ), 1, ZERO, Y( 1, I ), 1 )
               CALL DGEMV( <span class="string">'Transpose'</span>, I, N-I, -ONE, A( 1, I+1 ), LDA,
     $                     Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
               CALL DSCAL( N-I, TAUQ( I ), Y( I+1, I ), 1 )
            END IF
   20    CONTINUE
      END IF
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="DLABRD.288"></a><a href="dlabrd.f.html#DLABRD.1">DLABRD</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>