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 + -
显示快捷键?