cgebd2.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 275 行 · 第 1/2 页
HTML
275 行
</span><span class="comment">*</span><span class="comment"> where d and e denote diagonal and off-diagonal elements of B, vi
</span><span class="comment">*</span><span class="comment"> denotes an element of the vector defining H(i), and ui an element of
</span><span class="comment">*</span><span class="comment"> the vector defining G(i).
</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> COMPLEX ZERO, ONE
PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ),
$ ONE = ( 1.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
COMPLEX ALPHA
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Subroutines ..
</span> EXTERNAL <a name="CLACGV.136"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>, <a name="CLARF.136"></a><a href="clarf.f.html#CLARF.1">CLARF</a>, <a name="CLARFG.136"></a><a href="clarfg.f.html#CLARFG.1">CLARFG</a>, <a name="XERBLA.136"></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"> .. Intrinsic Functions ..
</span> INTRINSIC CONJG, MAX, 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"> Test the input parameters
</span><span class="comment">*</span><span class="comment">
</span> INFO = 0
IF( M.LT.0 ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
INFO = -4
END IF
IF( INFO.LT.0 ) THEN
CALL <a name="XERBLA.154"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CGEBD2.154"></a><a href="cgebd2.f.html#CGEBD2.1">CGEBD2</a>'</span>, -INFO )
RETURN
END IF
<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, N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Generate elementary reflector H(i) to annihilate A(i+1:m,i)
</span><span class="comment">*</span><span class="comment">
</span> ALPHA = A( I, I )
CALL <a name="CLARFG.167"></a><a href="clarfg.f.html#CLARFG.1">CLARFG</a>( M-I+1, ALPHA, A( MIN( I+1, M ), I ), 1,
$ TAUQ( I ) )
D( I ) = ALPHA
A( I, I ) = ONE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Apply H(i)' to A(i:m,i+1:n) from the left
</span><span class="comment">*</span><span class="comment">
</span> IF( I.LT.N )
$ CALL <a name="CLARF.175"></a><a href="clarf.f.html#CLARF.1">CLARF</a>( <span class="string">'Left'</span>, M-I+1, N-I, A( I, I ), 1,
$ CONJG( TAUQ( I ) ), A( I, I+1 ), LDA, WORK )
A( I, I ) = D( I )
<span class="comment">*</span><span class="comment">
</span> IF( I.LT.N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Generate elementary reflector G(i) to annihilate
</span><span class="comment">*</span><span class="comment"> A(i,i+2:n)
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="CLACGV.184"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( N-I, A( I, I+1 ), LDA )
ALPHA = A( I, I+1 )
CALL <a name="CLARFG.186"></a><a href="clarfg.f.html#CLARFG.1">CLARFG</a>( N-I, ALPHA, A( I, MIN( I+2, N ) ),
$ LDA, TAUP( I ) )
E( I ) = ALPHA
A( I, I+1 ) = ONE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Apply G(i) to A(i+1:m,i+1:n) from the right
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="CLARF.193"></a><a href="clarf.f.html#CLARF.1">CLARF</a>( <span class="string">'Right'</span>, M-I, N-I, A( I, I+1 ), LDA,
$ TAUP( I ), A( I+1, I+1 ), LDA, WORK )
CALL <a name="CLACGV.195"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( N-I, A( I, I+1 ), LDA )
A( I, I+1 ) = E( I )
ELSE
TAUP( I ) = ZERO
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, M
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Generate elementary reflector G(i) to annihilate A(i,i+1:n)
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="CLACGV.209"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( N-I+1, A( I, I ), LDA )
ALPHA = A( I, I )
CALL <a name="CLARFG.211"></a><a href="clarfg.f.html#CLARFG.1">CLARFG</a>( N-I+1, ALPHA, A( I, MIN( I+1, N ) ), LDA,
$ TAUP( I ) )
D( I ) = ALPHA
A( I, I ) = ONE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Apply G(i) to A(i+1:m,i:n) from the right
</span><span class="comment">*</span><span class="comment">
</span> IF( I.LT.M )
$ CALL <a name="CLARF.219"></a><a href="clarf.f.html#CLARF.1">CLARF</a>( <span class="string">'Right'</span>, M-I, N-I+1, A( I, I ), LDA,
$ TAUP( I ), A( I+1, I ), LDA, WORK )
CALL <a name="CLACGV.221"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( N-I+1, A( I, I ), LDA )
A( I, I ) = D( I )
<span class="comment">*</span><span class="comment">
</span> IF( I.LT.M ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Generate elementary reflector H(i) to annihilate
</span><span class="comment">*</span><span class="comment"> A(i+2:m,i)
</span><span class="comment">*</span><span class="comment">
</span> ALPHA = A( I+1, I )
CALL <a name="CLARFG.230"></a><a href="clarfg.f.html#CLARFG.1">CLARFG</a>( M-I, ALPHA, A( MIN( I+2, M ), I ), 1,
$ TAUQ( I ) )
E( I ) = ALPHA
A( I+1, I ) = ONE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Apply H(i)' to A(i+1:m,i+1:n) from the left
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="CLARF.237"></a><a href="clarf.f.html#CLARF.1">CLARF</a>( <span class="string">'Left'</span>, M-I, N-I, A( I+1, I ), 1,
$ CONJG( TAUQ( I ) ), A( I+1, I+1 ), LDA,
$ WORK )
A( I+1, I ) = E( I )
ELSE
TAUQ( I ) = ZERO
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="CGEBD2.248"></a><a href="cgebd2.f.html#CGEBD2.1">CGEBD2</a>
</span><span class="comment">*</span><span class="comment">
</span> END
</pre>
</body>
</html>
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?