slaqr0.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 665 行 · 第 1/4 页
HTML
665 行
</span><span class="comment">*</span><span class="comment"> where U is the orthogonal matrix in (*) (regard-
</span><span class="comment">*</span><span class="comment"> less of the value of WANTT.)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
</span><span class="comment">*</span><span class="comment"> accessed.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ================================================================
</span><span class="comment">*</span><span class="comment"> Based on contributions by
</span><span class="comment">*</span><span class="comment"> Karen Braman and Ralph Byers, Department of Mathematics,
</span><span class="comment">*</span><span class="comment"> University of Kansas, USA
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ================================================================
</span><span class="comment">*</span><span class="comment"> References:
</span><span class="comment">*</span><span class="comment"> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
</span><span class="comment">*</span><span class="comment"> Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
</span><span class="comment">*</span><span class="comment"> Performance, SIAM Journal of Matrix Analysis, volume 23, pages
</span><span class="comment">*</span><span class="comment"> 929--947, 2002.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
</span><span class="comment">*</span><span class="comment"> Algorithm Part II: Aggressive Early Deflation, SIAM Journal
</span><span class="comment">*</span><span class="comment"> of Matrix Analysis, volume 23, pages 948--973, 2002.
</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><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Matrices of order NTINY or smaller must be processed by
</span><span class="comment">*</span><span class="comment"> . <a name="SLAHQR.175"></a><a href="slahqr.f.html#SLAHQR.1">SLAHQR</a> because of insufficient subdiagonal scratch space.
</span><span class="comment">*</span><span class="comment"> . (This is a hard limit.) ====
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Exceptional deflation windows: try to cure rare
</span><span class="comment">*</span><span class="comment"> . slow convergence by increasing the size of the
</span><span class="comment">*</span><span class="comment"> . deflation window after KEXNW iterations. =====
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Exceptional shifts: try to cure rare slow convergence
</span><span class="comment">*</span><span class="comment"> . with ad-hoc exceptional shifts every KEXSH iterations.
</span><span class="comment">*</span><span class="comment"> . The constants WILK1 and WILK2 are used to form the
</span><span class="comment">*</span><span class="comment"> . exceptional shifts. ====
</span><span class="comment">*</span><span class="comment">
</span> INTEGER NTINY
PARAMETER ( NTINY = 11 )
INTEGER KEXNW, KEXSH
PARAMETER ( KEXNW = 5, KEXSH = 6 )
REAL WILK1, WILK2
PARAMETER ( WILK1 = 0.75e0, WILK2 = -0.4375e0 )
REAL ZERO, ONE
PARAMETER ( ZERO = 0.0e0, ONE = 1.0e0 )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Local Scalars ..
</span> REAL AA, BB, CC, CS, DD, SN, SS, SWAP
INTEGER I, INF, IT, ITMAX, K, KACC22, KBOT, KDU, KS,
$ KT, KTOP, KU, KV, KWH, KWTOP, KWV, LD, LS,
$ LWKOPT, NDFL, NH, NHO, NIBBLE, NMIN, NS, NSMAX,
$ NSR, NVE, NW, NWMAX, NWR
LOGICAL NWINC, SORTED
CHARACTER JBCMPZ*2
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Functions ..
</span> INTEGER <a name="ILAENV.206"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
EXTERNAL <a name="ILAENV.207"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Local Arrays ..
</span> REAL ZDUM( 1, 1 )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Subroutines ..
</span> EXTERNAL <a name="SLACPY.213"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>, <a name="SLAHQR.213"></a><a href="slahqr.f.html#SLAHQR.1">SLAHQR</a>, <a name="SLANV2.213"></a><a href="slanv2.f.html#SLANV2.1">SLANV2</a>, <a name="SLAQR3.213"></a><a href="slaqr3.f.html#SLAQR3.1">SLAQR3</a>, <a name="SLAQR4.213"></a><a href="slaqr4.f.html#SLAQR4.1">SLAQR4</a>, <a name="SLAQR5.213"></a><a href="slaqr5.f.html#SLAQR5.1">SLAQR5</a>
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Intrinsic Functions ..
</span> INTRINSIC ABS, INT, MAX, MIN, MOD, REAL
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Executable Statements ..
</span> INFO = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Quick return for N = 0: nothing to do. ====
</span><span class="comment">*</span><span class="comment">
</span> IF( N.EQ.0 ) THEN
WORK( 1 ) = ONE
RETURN
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Set up job flags for <a name="ILAENV.228"></a><a href="hfy-index.html#ILAENV">ILAENV</a>. ====
</span><span class="comment">*</span><span class="comment">
</span> IF( WANTT ) THEN
JBCMPZ( 1: 1 ) = <span class="string">'S'</span>
ELSE
JBCMPZ( 1: 1 ) = <span class="string">'E'</span>
END IF
IF( WANTZ ) THEN
JBCMPZ( 2: 2 ) = <span class="string">'V'</span>
ELSE
JBCMPZ( 2: 2 ) = <span class="string">'N'</span>
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Tiny matrices must use <a name="SLAHQR.241"></a><a href="slahqr.f.html#SLAHQR.1">SLAHQR</a>. ====
</span><span class="comment">*</span><span class="comment">
</span> IF( N.LE.NTINY ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Estimate optimal workspace. ====
</span><span class="comment">*</span><span class="comment">
</span> LWKOPT = 1
IF( LWORK.NE.-1 )
$ CALL <a name="SLAHQR.249"></a><a href="slahqr.f.html#SLAHQR.1">SLAHQR</a>( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
$ ILOZ, IHIZ, Z, LDZ, INFO )
ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Use small bulge multi-shift QR with aggressive early
</span><span class="comment">*</span><span class="comment"> . deflation on larger-than-tiny matrices. ====
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Hope for the best. ====
</span><span class="comment">*</span><span class="comment">
</span> INFO = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== NWR = recommended deflation window size. At this
</span><span class="comment">*</span><span class="comment"> . point, N .GT. NTINY = 11, so there is enough
</span><span class="comment">*</span><span class="comment"> . subdiagonal workspace for NWR.GE.2 as required.
</span><span class="comment">*</span><span class="comment"> . (In fact, there is enough subdiagonal space for
</span><span class="comment">*</span><span class="comment"> . NWR.GE.3.) ====
</span><span class="comment">*</span><span class="comment">
</span> NWR = <a name="ILAENV.266"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 13, <span class="string">'<a name="SLAQR0.266"></a><a href="slaqr0.f.html#SLAQR0.1">SLAQR0</a>'</span>, JBCMPZ, N, ILO, IHI, LWORK )
NWR = MAX( 2, NWR )
NWR = MIN( IHI-ILO+1, ( N-1 ) / 3, NWR )
NW = NWR
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== NSR = recommended number of simultaneous shifts.
</span><span class="comment">*</span><span class="comment"> . At this point N .GT. NTINY = 11, so there is at
</span><span class="comment">*</span><span class="comment"> . enough subdiagonal workspace for NSR to be even
</span><span class="comment">*</span><span class="comment"> . and greater than or equal to two as required. ====
</span><span class="comment">*</span><span class="comment">
</span> NSR = <a name="ILAENV.276"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 15, <span class="string">'<a name="SLAQR0.276"></a><a href="slaqr0.f.html#SLAQR0.1">SLAQR0</a>'</span>, JBCMPZ, N, ILO, IHI, LWORK )
NSR = MIN( NSR, ( N+6 ) / 9, IHI-ILO )
NSR = MAX( 2, NSR-MOD( NSR, 2 ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Estimate optimal workspace ====
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Workspace query call to <a name="SLAQR3.282"></a><a href="slaqr3.f.html#SLAQR3.1">SLAQR3</a> ====
</span><span class="comment">*</span><span class="comment">
</span> CALL <a name="SLAQR3.284"></a><a href="slaqr3.f.html#SLAQR3.1">SLAQR3</a>( WANTT, WANTZ, N, ILO, IHI, NWR+1, H, LDH, ILOZ,
$ IHIZ, Z, LDZ, LS, LD, WR, WI, H, LDH, N, H, LDH,
$ N, H, LDH, WORK, -1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Optimal workspace = MAX(<a name="SLAQR5.288"></a><a href="slaqr5.f.html#SLAQR5.1">SLAQR5</a>, <a name="SLAQR3.288"></a><a href="slaqr3.f.html#SLAQR3.1">SLAQR3</a>) ====
</span><span class="comment">*</span><span class="comment">
</span> LWKOPT = MAX( 3*NSR / 2, INT( WORK( 1 ) ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Quick return in case of workspace query. ====
</span><span class="comment">*</span><span class="comment">
</span> IF( LWORK.EQ.-1 ) THEN
WORK( 1 ) = REAL( LWKOPT )
RETURN
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== <a name="SLAHQR.299"></a><a href="slahqr.f.html#SLAHQR.1">SLAHQR</a>/<a name="SLAQR0.299"></a><a href="slaqr0.f.html#SLAQR0.1">SLAQR0</a> crossover point ====
</span><span class="comment">*</span><span class="comment">
</span> NMIN = <a name="ILAENV.301"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 12, <span class="string">'<a name="SLAQR0.301"></a><a href="slaqr0.f.html#SLAQR0.1">SLAQR0</a>'</span>, JBCMPZ, N, ILO, IHI, LWORK )
NMIN = MAX( NTINY, NMIN )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Nibble crossover point ====
</span><span class="comment">*</span><span class="comment">
</span> NIBBLE = <a name="ILAENV.306"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 14, <span class="string">'<a name="SLAQR0.306"></a><a href="slaqr0.f.html#SLAQR0.1">SLAQR0</a>'</span>, JBCMPZ, N, ILO, IHI, LWORK )
NIBBLE = MAX( 0, NIBBLE )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ==== Accumulate reflections during ttswp? Use block
</span><span class="comment">*</span><span class="comment"> . 2-by-2 structure during matrix-matrix multiply? ====
</span><span class="comment">*</span><span class="comment">
</span> KACC22 = <a name="ILAENV.312"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 16, <span class="string">'<a name="SLAQR0.312"></a><a href="slaqr0.f.html#SLAQR0.1">SLAQR0</a>'</span>, JBCMPZ, N, ILO, IHI, LWORK )
KACC22 = MAX( 0, KACC22 )
KACC22 = MIN( 2, KACC22 )
<span class="comment">*</span><span class="comment">
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?