slaqr5.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 837 行 · 第 1/5 页
HTML
837 行
$ JROW, JTOP, K, K1, KDU, KMS, KNZ, KRCOL, KZS,
$ M, M22, MBOT, MEND, MSTART, MTOP, NBMPS, NDCOL,
$ NS, NU
LOGICAL ACCUM, BLK22, BMP22
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Functions ..
</span> REAL <a name="SLAMCH.156"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
EXTERNAL <a name="SLAMCH.157"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Intrinsic Functions ..
</span><span class="comment">*</span><span class="comment">
</span> INTRINSIC ABS, MAX, MIN, MOD, REAL
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Local Arrays ..
</span> REAL VT( 3 )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Subroutines ..
</span> EXTERNAL SGEMM, <a name="SLABAD.167"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>, <a name="SLACPY.167"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>, <a name="SLAQR1.167"></a><a href="slaqr1.f.html#SLAQR1.1">SLAQR1</a>, <a name="SLARFG.167"></a><a href="slarfg.f.html#SLARFG.1">SLARFG</a>, <a name="SLASET.167"></a><a href="slaset.f.html#SLASET.1">SLASET</a>,
$ STRMM
<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"> ==== If there are no shifts, then there is nothing to do. ====
</span><span class="comment">*</span><span class="comment">
</span> IF( NSHFTS.LT.2 )
$ RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== If the active block is empty or 1-by-1, then there
</span><span class="comment">*</span><span class="comment"> . is nothing to do. ====
</span><span class="comment">*</span><span class="comment">
</span> IF( KTOP.GE.KBOT )
$ RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Shuffle shifts into pairs of real shifts and pairs
</span><span class="comment">*</span><span class="comment"> . of complex conjugate shifts assuming complex
</span><span class="comment">*</span><span class="comment"> . conjugate shifts are already adjacent to one
</span><span class="comment">*</span><span class="comment"> . another. ====
</span><span class="comment">*</span><span class="comment">
</span> DO 10 I = 1, NSHFTS - 2, 2
IF( SI( I ).NE.-SI( I+1 ) ) THEN
<span class="comment">*</span><span class="comment">
</span> SWAP = SR( I )
SR( I ) = SR( I+1 )
SR( I+1 ) = SR( I+2 )
SR( I+2 ) = SWAP
<span class="comment">*</span><span class="comment">
</span> SWAP = SI( I )
SI( I ) = SI( I+1 )
SI( I+1 ) = SI( I+2 )
SI( I+2 ) = SWAP
END IF
10 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== NSHFTS is supposed to be even, but if is odd,
</span><span class="comment">*</span><span class="comment"> . then simply reduce it by one. The shuffle above
</span><span class="comment">*</span><span class="comment"> . ensures that the dropped shift is real and that
</span><span class="comment">*</span><span class="comment"> . the remaining shifts are paired. ====
</span><span class="comment">*</span><span class="comment">
</span> NS = NSHFTS - MOD( NSHFTS, 2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Machine constants for deflation ====
</span><span class="comment">*</span><span class="comment">
</span> SAFMIN = <a name="SLAMCH.212"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'SAFE MINIMUM'</span> )
SAFMAX = ONE / SAFMIN
CALL <a name="SLABAD.214"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>( SAFMIN, SAFMAX )
ULP = <a name="SLAMCH.215"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'PRECISION'</span> )
SMLNUM = SAFMIN*( REAL( N ) / ULP )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Use accumulated reflections to update far-from-diagonal
</span><span class="comment">*</span><span class="comment"> . entries ? ====
</span><span class="comment">*</span><span class="comment">
</span> ACCUM = ( KACC22.EQ.1 ) .OR. ( KACC22.EQ.2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== If so, exploit the 2-by-2 block structure? ====
</span><span class="comment">*</span><span class="comment">
</span> BLK22 = ( NS.GT.2 ) .AND. ( KACC22.EQ.2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== clear trash ====
</span><span class="comment">*</span><span class="comment">
</span> IF( KTOP+2.LE.KBOT )
$ H( KTOP+2, KTOP ) = ZERO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== NBMPS = number of 2-shift bulges in the chain ====
</span><span class="comment">*</span><span class="comment">
</span> NBMPS = NS / 2
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== KDU = width of slab ====
</span><span class="comment">*</span><span class="comment">
</span> KDU = 6*NBMPS - 3
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Create and chase chains of NBMPS bulges ====
</span><span class="comment">*</span><span class="comment">
</span> DO 220 INCOL = 3*( 1-NBMPS ) + KTOP - 1, KBOT - 2, 3*NBMPS - 2
NDCOL = INCOL + KDU
IF( ACCUM )
$ CALL <a name="SLASET.245"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'ALL'</span>, KDU, KDU, ZERO, ONE, U, LDU )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Near-the-diagonal bulge chase. The following loop
</span><span class="comment">*</span><span class="comment"> . performs the near-the-diagonal part of a small bulge
</span><span class="comment">*</span><span class="comment"> . multi-shift QR sweep. Each 6*NBMPS-2 column diagonal
</span><span class="comment">*</span><span class="comment"> . chunk extends from column INCOL to column NDCOL
</span><span class="comment">*</span><span class="comment"> . (including both column INCOL and column NDCOL). The
</span><span class="comment">*</span><span class="comment"> . following loop chases a 3*NBMPS column long chain of
</span><span class="comment">*</span><span class="comment"> . NBMPS bulges 3*NBMPS-2 columns to the right. (INCOL
</span><span class="comment">*</span><span class="comment"> . may be less than KTOP and and NDCOL may be greater than
</span><span class="comment">*</span><span class="comment"> . KBOT indicating phantom columns from which to chase
</span><span class="comment">*</span><span class="comment"> . bulges before they are actually introduced or to which
</span><span class="comment">*</span><span class="comment"> . to chase bulges beyond column KBOT.) ====
</span><span class="comment">*</span><span class="comment">
</span> DO 150 KRCOL = INCOL, MIN( INCOL+3*NBMPS-3, KBOT-2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Bulges number MTOP to MBOT are active double implicit
</span><span class="comment">*</span><span class="comment"> . shift bulges. There may or may not also be small
</span><span class="comment">*</span><span class="comment"> . 2-by-2 bulge, if there is room. The inactive bulges
</span><span class="comment">*</span><span class="comment"> . (if any) must wait until the active bulges have moved
</span><span class="comment">*</span><span class="comment"> . down the diagonal to make room. The phantom matrix
</span><span class="comment">*</span><span class="comment"> . paradigm described above helps keep track. ====
</span><span class="comment">*</span><span class="comment">
</span> MTOP = MAX( 1, ( ( KTOP-1 )-KRCOL+2 ) / 3+1 )
MBOT = MIN( NBMPS, ( KBOT-KRCOL ) / 3 )
M22 = MBOT + 1
BMP22 = ( MBOT.LT.NBMPS ) .AND. ( KRCOL+3*( M22-1 ) ).EQ.
$ ( KBOT-2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Generate reflections to chase the chain right
</span><span class="comment">*</span><span class="comment"> . one column. (The minimum value of K is KTOP-1.) ====
</span><span class="comment">*</span><span class="comment">
</span> DO 20 M = MTOP, MBOT
K = KRCOL + 3*( M-1 )
IF( K.EQ.KTOP-1 ) THEN
CALL <a name="SLAQR1.280"></a><a href="slaqr1.f.html#SLAQR1.1">SLAQR1</a>( 3, H( KTOP, KTOP ), LDH, SR( 2*M-1 ),
$ SI( 2*M-1 ), SR( 2*M ), SI( 2*M ),
$ V( 1, M ) )
ALPHA = V( 1, M )
CALL <a name="SLARFG.284"></a><a href="slarfg.f.html#SLARFG.1">SLARFG</a>( 3, ALPHA, V( 2, M ), 1, V( 1, M ) )
ELSE
BETA = H( K+1, K )
V( 2, M ) = H( K+2, K )
V( 3, M ) = H( K+3, K )
CALL <a name="SLARFG.289"></a><a href="slarfg.f.html#SLARFG.1">SLARFG</a>( 3, BETA, V( 2, M ), 1, V( 1, M ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== A Bulge may collapse because of vigilant
</span><span class="comment">*</span><span class="comment"> . deflation or destructive underflow. (The
</span><span class="comment">*</span><span class="comment"> . initial bulge is always collapsed.) Use
</span><span class="comment">*</span><span class="comment"> . the two-small-subdiagonals trick to try
</span><span class="comment">*</span><span class="comment"> . to get it started again. If V(2,M).NE.0 and
</span><span class="comment">*</span><span class="comment"> . V(3,M) = H(K+3,K+1) = H(K+3,K+2) = 0, then
</span><span class="comment">*</span><span class="comment"> . this bulge is collapsing into a zero
</span><span class="comment">*</span><span class="comment"> . subdiagonal. It will be restarted next
</span><span class="comment">*</span><span class="comment"> . trip through the loop.)
</span><span class="comment">*</span><span class="comment">
</span> IF( V( 1, M ).NE.ZERO .AND.
$ ( V( 3, M ).NE.ZERO .OR. ( H( K+3,
$ K+1 ).EQ.ZERO .AND. H( K+3, K+2 ).EQ.ZERO ) ) )
$ THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Typical case: not collapsed (yet). ====
</span><span class="comment">*</span><span class="comment">
</span> H( K+1, K ) = BETA
H( K+2, K ) = ZERO
H( K+3, K ) = ZERO
ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Atypical case: collapsed. Attempt to
</span><span class="comment">*</span><span class="comment"> . reintroduce ignoring H(K+1,K). If the
</span><span class="comment">*</span><span class="comment"> . fill resulting from the new reflector
</span><span class="comment">*</span><span class="comment"> . is too large, then abandon it.
</span><span class="comment">*</span><span class="comment"> . Otherwise, use the new one. ====
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?