sbdsqr.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 767 行 · 第 1/4 页
HTML
767 行
E( LL ) = H*OLDSN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Update singular vectors
</span><span class="comment">*</span><span class="comment">
</span> IF( NCVT.GT.0 )
$ CALL <a name="SLASR.563"></a><a href="slasr.f.html#SLASR.1">SLASR</a>( <span class="string">'L'</span>, <span class="string">'V'</span>, <span class="string">'B'</span>, M-LL+1, NCVT, WORK( NM12+1 ),
$ WORK( NM13+1 ), VT( LL, 1 ), LDVT )
IF( NRU.GT.0 )
$ CALL <a name="SLASR.566"></a><a href="slasr.f.html#SLASR.1">SLASR</a>( <span class="string">'R'</span>, <span class="string">'V'</span>, <span class="string">'B'</span>, NRU, M-LL+1, WORK( 1 ),
$ WORK( N ), U( 1, LL ), LDU )
IF( NCC.GT.0 )
$ CALL <a name="SLASR.569"></a><a href="slasr.f.html#SLASR.1">SLASR</a>( <span class="string">'L'</span>, <span class="string">'V'</span>, <span class="string">'B'</span>, M-LL+1, NCC, WORK( 1 ),
$ WORK( N ), C( LL, 1 ), LDC )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Test convergence
</span><span class="comment">*</span><span class="comment">
</span> IF( ABS( E( LL ) ).LE.THRESH )
$ E( LL ) = ZERO
END IF
ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Use nonzero shift
</span><span class="comment">*</span><span class="comment">
</span> IF( IDIR.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Chase bulge from top to bottom
</span><span class="comment">*</span><span class="comment"> Save cosines and sines for later singular vector updates
</span><span class="comment">*</span><span class="comment">
</span> F = ( ABS( D( LL ) )-SHIFT )*
$ ( SIGN( ONE, D( LL ) )+SHIFT / D( LL ) )
G = E( LL )
DO 140 I = LL, M - 1
CALL <a name="SLARTG.590"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( F, G, COSR, SINR, R )
IF( I.GT.LL )
$ E( I-1 ) = R
F = COSR*D( I ) + SINR*E( I )
E( I ) = COSR*E( I ) - SINR*D( I )
G = SINR*D( I+1 )
D( I+1 ) = COSR*D( I+1 )
CALL <a name="SLARTG.597"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( F, G, COSL, SINL, R )
D( I ) = R
F = COSL*E( I ) + SINL*D( I+1 )
D( I+1 ) = COSL*D( I+1 ) - SINL*E( I )
IF( I.LT.M-1 ) THEN
G = SINL*E( I+1 )
E( I+1 ) = COSL*E( I+1 )
END IF
WORK( I-LL+1 ) = COSR
WORK( I-LL+1+NM1 ) = SINR
WORK( I-LL+1+NM12 ) = COSL
WORK( I-LL+1+NM13 ) = SINL
140 CONTINUE
E( M-1 ) = F
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Update singular vectors
</span><span class="comment">*</span><span class="comment">
</span> IF( NCVT.GT.0 )
$ CALL <a name="SLASR.615"></a><a href="slasr.f.html#SLASR.1">SLASR</a>( <span class="string">'L'</span>, <span class="string">'V'</span>, <span class="string">'F'</span>, M-LL+1, NCVT, WORK( 1 ),
$ WORK( N ), VT( LL, 1 ), LDVT )
IF( NRU.GT.0 )
$ CALL <a name="SLASR.618"></a><a href="slasr.f.html#SLASR.1">SLASR</a>( <span class="string">'R'</span>, <span class="string">'V'</span>, <span class="string">'F'</span>, NRU, M-LL+1, WORK( NM12+1 ),
$ WORK( NM13+1 ), U( 1, LL ), LDU )
IF( NCC.GT.0 )
$ CALL <a name="SLASR.621"></a><a href="slasr.f.html#SLASR.1">SLASR</a>( <span class="string">'L'</span>, <span class="string">'V'</span>, <span class="string">'F'</span>, M-LL+1, NCC, WORK( NM12+1 ),
$ WORK( NM13+1 ), C( LL, 1 ), LDC )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Test convergence
</span><span class="comment">*</span><span class="comment">
</span> IF( ABS( E( M-1 ) ).LE.THRESH )
$ E( M-1 ) = ZERO
<span class="comment">*</span><span class="comment">
</span> ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Chase bulge from bottom to top
</span><span class="comment">*</span><span class="comment"> Save cosines and sines for later singular vector updates
</span><span class="comment">*</span><span class="comment">
</span> F = ( ABS( D( M ) )-SHIFT )*( SIGN( ONE, D( M ) )+SHIFT /
$ D( M ) )
G = E( M-1 )
DO 150 I = M, LL + 1, -1
CALL <a name="SLARTG.638"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( F, G, COSR, SINR, R )
IF( I.LT.M )
$ E( I ) = R
F = COSR*D( I ) + SINR*E( I-1 )
E( I-1 ) = COSR*E( I-1 ) - SINR*D( I )
G = SINR*D( I-1 )
D( I-1 ) = COSR*D( I-1 )
CALL <a name="SLARTG.645"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( F, G, COSL, SINL, R )
D( I ) = R
F = COSL*E( I-1 ) + SINL*D( I-1 )
D( I-1 ) = COSL*D( I-1 ) - SINL*E( I-1 )
IF( I.GT.LL+1 ) THEN
G = SINL*E( I-2 )
E( I-2 ) = COSL*E( I-2 )
END IF
WORK( I-LL ) = COSR
WORK( I-LL+NM1 ) = -SINR
WORK( I-LL+NM12 ) = COSL
WORK( I-LL+NM13 ) = -SINL
150 CONTINUE
E( LL ) = F
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Test convergence
</span><span class="comment">*</span><span class="comment">
</span> IF( ABS( E( LL ) ).LE.THRESH )
$ E( LL ) = ZERO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Update singular vectors if desired
</span><span class="comment">*</span><span class="comment">
</span> IF( NCVT.GT.0 )
$ CALL <a name="SLASR.668"></a><a href="slasr.f.html#SLASR.1">SLASR</a>( <span class="string">'L'</span>, <span class="string">'V'</span>, <span class="string">'B'</span>, M-LL+1, NCVT, WORK( NM12+1 ),
$ WORK( NM13+1 ), VT( LL, 1 ), LDVT )
IF( NRU.GT.0 )
$ CALL <a name="SLASR.671"></a><a href="slasr.f.html#SLASR.1">SLASR</a>( <span class="string">'R'</span>, <span class="string">'V'</span>, <span class="string">'B'</span>, NRU, M-LL+1, WORK( 1 ),
$ WORK( N ), U( 1, LL ), LDU )
IF( NCC.GT.0 )
$ CALL <a name="SLASR.674"></a><a href="slasr.f.html#SLASR.1">SLASR</a>( <span class="string">'L'</span>, <span class="string">'V'</span>, <span class="string">'B'</span>, M-LL+1, NCC, WORK( 1 ),
$ WORK( N ), C( LL, 1 ), LDC )
END IF
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> QR iteration finished, go back and check convergence
</span><span class="comment">*</span><span class="comment">
</span> GO TO 60
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> All singular values converged, so make them positive
</span><span class="comment">*</span><span class="comment">
</span> 160 CONTINUE
DO 170 I = 1, N
IF( D( I ).LT.ZERO ) THEN
D( I ) = -D( I )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Change sign of singular vectors, if desired
</span><span class="comment">*</span><span class="comment">
</span> IF( NCVT.GT.0 )
$ CALL SSCAL( NCVT, NEGONE, VT( I, 1 ), LDVT )
END IF
170 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Sort the singular values into decreasing order (insertion sort on
</span><span class="comment">*</span><span class="comment"> singular values, but only one transposition per singular vector)
</span><span class="comment">*</span><span class="comment">
</span> DO 190 I = 1, N - 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Scan for smallest D(I)
</span><span class="comment">*</span><span class="comment">
</span> ISUB = 1
SMIN = D( 1 )
DO 180 J = 2, N + 1 - I
IF( D( J ).LE.SMIN ) THEN
ISUB = J
SMIN = D( J )
END IF
180 CONTINUE
IF( ISUB.NE.N+1-I ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Swap singular values and vectors
</span><span class="comment">*</span><span class="comment">
</span> D( ISUB ) = D( N+1-I )
D( N+1-I ) = SMIN
IF( NCVT.GT.0 )
$ CALL SSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( N+1-I, 1 ),
$ LDVT )
IF( NRU.GT.0 )
$ CALL SSWAP( NRU, U( 1, ISUB ), 1, U( 1, N+1-I ), 1 )
IF( NCC.GT.0 )
$ CALL SSWAP( NCC, C( ISUB, 1 ), LDC, C( N+1-I, 1 ), LDC )
END IF
190 CONTINUE
GO TO 220
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Maximum number of iterations exceeded, failure to converge
</span><span class="comment">*</span><span class="comment">
</span> 200 CONTINUE
INFO = 0
DO 210 I = 1, N - 1
IF( E( I ).NE.ZERO )
$ INFO = INFO + 1
210 CONTINUE
220 CONTINUE
RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> End of <a name="SBDSQR.740"></a><a href="sbdsqr.f.html#SBDSQR.1">SBDSQR</a>
</span><span class="comment">*</span><span class="comment">
</span> END
</pre>
</body>
</html>
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?