claqr0.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 626 行 · 第 1/4 页

HTML
626
字号
<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.299"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 14, <span class="string">'<a name="CLAQR0.299"></a><a href="claqr0.f.html#CLAQR0.1">CLAQR0</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.305"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 16, <span class="string">'<a name="CLAQR0.305"></a><a href="claqr0.f.html#CLAQR0.1">CLAQR0</a>'</span>, JBCMPZ, N, ILO, IHI, LWORK )
         KACC22 = MAX( 0, KACC22 )
         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 70 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 80
<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( CABS1( H( KWTOP, KWTOP-1 ) ).GT.
     $                      CABS1( 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="CLAQR3.415"></a><a href="claqr3.f.html#CLAQR3.1">CLAQR3</a>( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
     $                   IHIZ, Z, LDZ, LS, LD, W, 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="CLAQR3.438"></a><a href="claqr3.f.html#CLAQR3.1">CLAQR3</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="CLAQR3.448"></a><a href="claqr3.f.html#CLAQR3.1">CLAQR3</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

⌨️ 快捷键说明

复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?