dlaqr0.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 667 行 · 第 1/4 页
HTML
667 行
KACC22 = MIN( 2, KACC22 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== NWMAX = the largest possible deflation window for
</span><span class="comment">*</span><span class="comment"> . which there is sufficient workspace. ====
</span><span class="comment">*</span><span class="comment">
</span> NWMAX = MIN( ( N-1 ) / 3, LWORK / 2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== NSMAX = the Largest number of simultaneous shifts
</span><span class="comment">*</span><span class="comment"> . for which there is sufficient workspace. ====
</span><span class="comment">*</span><span class="comment">
</span> NSMAX = MIN( ( N+6 ) / 9, 2*LWORK / 3 )
NSMAX = NSMAX - MOD( NSMAX, 2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== NDFL: an iteration count restarted at deflation. ====
</span><span class="comment">*</span><span class="comment">
</span> NDFL = 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== ITMAX = iteration limit ====
</span><span class="comment">*</span><span class="comment">
</span> ITMAX = MAX( 30, 2*KEXSH )*MAX( 10, ( IHI-ILO+1 ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Last row and column in the active block ====
</span><span class="comment">*</span><span class="comment">
</span> KBOT = IHI
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Main Loop ====
</span><span class="comment">*</span><span class="comment">
</span> DO 80 IT = 1, ITMAX
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Done when KBOT falls below ILO ====
</span><span class="comment">*</span><span class="comment">
</span> IF( KBOT.LT.ILO )
$ GO TO 90
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Locate active block ====
</span><span class="comment">*</span><span class="comment">
</span> DO 10 K = KBOT, ILO + 1, -1
IF( H( K, K-1 ).EQ.ZERO )
$ GO TO 20
10 CONTINUE
K = ILO
20 CONTINUE
KTOP = K
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Select deflation window size ====
</span><span class="comment">*</span><span class="comment">
</span> NH = KBOT - KTOP + 1
IF( NDFL.LT.KEXNW .OR. NH.LT.NW ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Typical deflation window. If possible and
</span><span class="comment">*</span><span class="comment"> . advisable, nibble the entire active block.
</span><span class="comment">*</span><span class="comment"> . If not, use size NWR or NWR+1 depending upon
</span><span class="comment">*</span><span class="comment"> . which has the smaller corresponding subdiagonal
</span><span class="comment">*</span><span class="comment"> . entry (a heuristic). ====
</span><span class="comment">*</span><span class="comment">
</span> NWINC = .TRUE.
IF( NH.LE.MIN( NMIN, NWMAX ) ) THEN
NW = NH
ELSE
NW = MIN( NWR, NH, NWMAX )
IF( NW.LT.NWMAX ) THEN
IF( NW.GE.NH-1 ) THEN
NW = NH
ELSE
KWTOP = KBOT - NW + 1
IF( ABS( H( KWTOP, KWTOP-1 ) ).GT.
$ ABS( H( KWTOP-1, KWTOP-2 ) ) )NW = NW + 1
END IF
END IF
END IF
ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Exceptional deflation window. If there have
</span><span class="comment">*</span><span class="comment"> . been no deflations in KEXNW or more iterations,
</span><span class="comment">*</span><span class="comment"> . then vary the deflation window size. At first,
</span><span class="comment">*</span><span class="comment"> . because, larger windows are, in general, more
</span><span class="comment">*</span><span class="comment"> . powerful than smaller ones, rapidly increase the
</span><span class="comment">*</span><span class="comment"> . window up to the maximum reasonable and possible.
</span><span class="comment">*</span><span class="comment"> . Then maybe try a slightly smaller window. ====
</span><span class="comment">*</span><span class="comment">
</span> IF( NWINC .AND. NW.LT.MIN( NWMAX, NH ) ) THEN
NW = MIN( NWMAX, NH, 2*NW )
ELSE
NWINC = .FALSE.
IF( NW.EQ.NH .AND. NH.GT.2 )
$ NW = NH - 1
END IF
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Aggressive early deflation:
</span><span class="comment">*</span><span class="comment"> . split workspace under the subdiagonal into
</span><span class="comment">*</span><span class="comment"> . - an nw-by-nw work array V in the lower
</span><span class="comment">*</span><span class="comment"> . left-hand-corner,
</span><span class="comment">*</span><span class="comment"> . - an NW-by-at-least-NW-but-more-is-better
</span><span class="comment">*</span><span class="comment"> . (NW-by-NHO) horizontal work array along
</span><span class="comment">*</span><span class="comment"> . the bottom edge,
</span><span class="comment">*</span><span class="comment"> . - an at-least-NW-but-more-is-better (NHV-by-NW)
</span><span class="comment">*</span><span class="comment"> . vertical work array along the left-hand-edge.
</span><span class="comment">*</span><span class="comment"> . ====
</span><span class="comment">*</span><span class="comment">
</span> KV = N - NW + 1
KT = NW + 1
NHO = ( N-NW-1 ) - KT + 1
KWV = NW + 2
NVE = ( N-NW ) - KWV + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Aggressive early deflation ====
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="DLAQR3.424"></a><a href="dlaqr3.f.html#DLAQR3.1">DLAQR3</a>( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
$ IHIZ, Z, LDZ, LS, LD, WR, WI, H( KV, 1 ), LDH,
$ NHO, H( KV, KT ), LDH, NVE, H( KWV, 1 ), LDH,
$ WORK, LWORK )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Adjust KBOT accounting for new deflations. ====
</span><span class="comment">*</span><span class="comment">
</span> KBOT = KBOT - LD
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== KS points to the shifts. ====
</span><span class="comment">*</span><span class="comment">
</span> KS = KBOT - LS + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Skip an expensive QR sweep if there is a (partly
</span><span class="comment">*</span><span class="comment"> . heuristic) reason to expect that many eigenvalues
</span><span class="comment">*</span><span class="comment"> . will deflate without it. Here, the QR sweep is
</span><span class="comment">*</span><span class="comment"> . skipped if many eigenvalues have just been deflated
</span><span class="comment">*</span><span class="comment"> . or if the remaining active block is small.
</span><span class="comment">*</span><span class="comment">
</span> IF( ( LD.EQ.0 ) .OR. ( ( 100*LD.LE.NW*NIBBLE ) .AND. ( KBOT-
$ KTOP+1.GT.MIN( NMIN, NWMAX ) ) ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== NS = nominal number of simultaneous shifts.
</span><span class="comment">*</span><span class="comment"> . This may be lowered (slightly) if <a name="DLAQR3.447"></a><a href="dlaqr3.f.html#DLAQR3.1">DLAQR3</a>
</span><span class="comment">*</span><span class="comment"> . did not provide that many shifts. ====
</span><span class="comment">*</span><span class="comment">
</span> NS = MIN( NSMAX, NSR, MAX( 2, KBOT-KTOP ) )
NS = NS - MOD( NS, 2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== If there have been no deflations
</span><span class="comment">*</span><span class="comment"> . in a multiple of KEXSH iterations,
</span><span class="comment">*</span><span class="comment"> . then try exceptional shifts.
</span><span class="comment">*</span><span class="comment"> . Otherwise use shifts provided by
</span><span class="comment">*</span><span class="comment"> . <a name="DLAQR3.457"></a><a href="dlaqr3.f.html#DLAQR3.1">DLAQR3</a> above or from the eigenvalues
</span><span class="comment">*</span><span class="comment"> . of a trailing principal submatrix. ====
</span><span class="comment">*</span><span class="comment">
</span> IF( MOD( NDFL, KEXSH ).EQ.0 ) THEN
KS = KBOT - NS + 1
DO 30 I = KBOT, MAX( KS+1, KTOP+2 ), -2
SS = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) )
AA = WILK1*SS + H( I, I )
BB = SS
CC = WILK2*SS
DD = AA
CALL <a name="DLANV2.468"></a><a href="dlanv2.f.html#DLANV2.1">DLANV2</a>( AA, BB, CC, DD, WR( I-1 ), WI( I-1 ),
$ WR( I ), WI( I ), CS, SN )
30 CONTINUE
IF( KS.EQ.KTOP ) THEN
WR( KS+1 ) = H( KS+1, KS+1 )
WI( KS+1 ) = ZERO
WR( KS ) = WR( KS+1 )
WI( KS ) = WI( KS+1 )
END IF
ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Got NS/2 or fewer shifts? Use <a name="DLAQR4.479"></a><a href="dlaqr4.f.html#DLAQR4.1">DLAQR4</a> or
</span><span class="comment">*</span><span class="comment"> . <a name="DLAHQR.480"></a><a href="dlahqr.f.html#DLAHQR.1">DLAHQR</a> on a trailing principal submatrix to
</span><span class="comment">*</span><span class="comment"> . get more. (Since NS.LE.NSMAX.LE.(N+6)/9,
</span><span class="comment">*</span><span class="comment"> . there is enough space below the subdiagonal
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?