slaqr0.f.html

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

HTML
665 行
字号
1234
</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="SLAQR3.422"></a><a href="slaqr3.f.html#SLAQR3.1">SLAQR3</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="SLAQR3.445"></a><a href="slaqr3.f.html#SLAQR3.1">SLAQR3</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="SLAQR3.455"></a><a href="slaqr3.f.html#SLAQR3.1">SLAQR3</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="SLANV2.466"></a><a href="slanv2.f.html#SLANV2.1">SLANV2</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="SLAQR4.477"></a><a href="slaqr4.f.html#SLAQR4.1">SLAQR4</a> or
</span><span class="comment">*</span><span class="comment">                 .    <a name="SLAHQR.478"></a><a href="slahqr.f.html#SLAHQR.1">SLAHQR</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
</span><span class="comment">*</span><span class="comment">                 .    to fit an NS-by-NS scratch array.) ====
</span><span class="comment">*</span><span class="comment">
1234

slaqr0.f.html - 源码说明

本页面展示了「famous linear algebra library (LAPACK) ports to windows」中的 slaqr0.f.html 源码文件，采用 HTML 编程语言编写，共 665 行代码。您可以在线阅读完整代码内容，也可以返回资源详情页下载完整源码包进行本地学习和开发。

虫虫下载站收录了大量与algebra相关的技术资源，包括源代码、技术文档、电路图等，是电子工程师和嵌入式开发者的专业学习平台。

⌨️ 快捷键说明

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