dlaqr4.f.html

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

<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 ==== Got NS/2 or fewer shifts? Use <a name="DLAHQR.484"></a><a href="dlahqr.f.html#DLAHQR.1">DLAHQR</a>
</span><span class="comment">*</span><span class="comment">                 .    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">
</span>                  IF( KBOT-KS+1.LE.NS / 2 ) THEN
                     KS = KBOT - NS + 1
                     KT = N - NS + 1
                     CALL <a name="DLACPY.493"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'A'</span>, NS, NS, H( KS, KS ), LDH,
     $                            H( KT, 1 ), LDH )
                     CALL <a name="DLAHQR.495"></a><a href="dlahqr.f.html#DLAHQR.1">DLAHQR</a>( .false., .false., NS, 1, NS,
     $                            H( KT, 1 ), LDH, WR( KS ), WI( KS ),
     $                            1, 1, ZDUM, 1, INF )
                     KS = KS + INF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    ==== In case of a rare QR failure use
</span><span class="comment">*</span><span class="comment">                    .    eigenvalues of the trailing 2-by-2
</span><span class="comment">*</span><span class="comment">                    .    principal submatrix.  ====
</span><span class="comment">*</span><span class="comment">
</span>                     IF( KS.GE.KBOT ) THEN
                        AA = H( KBOT-1, KBOT-1 )
                        CC = H( KBOT, KBOT-1 )
                        BB = H( KBOT-1, KBOT )
                        DD = H( KBOT, KBOT )
                        CALL <a name="DLANV2.509"></a><a href="dlanv2.f.html#DLANV2.1">DLANV2</a>( AA, BB, CC, DD, WR( KBOT-1 ),
     $                               WI( KBOT-1 ), WR( KBOT ),
     $                               WI( KBOT ), CS, SN )
                        KS = KBOT - 1
                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  IF( KBOT-KS+1.GT.NS ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    ==== Sort the shifts (Helps a little)
</span><span class="comment">*</span><span class="comment">                    .    Bubble sort keeps complex conjugate
</span><span class="comment">*</span><span class="comment">                    .    pairs together. ====
</span><span class="comment">*</span><span class="comment">
</span>                     SORTED = .false.
                     DO 50 K = KBOT, KS + 1, -1
                        IF( SORTED )
     $                     GO TO 60
                        SORTED = .true.
                        DO 40 I = KS, K - 1
                           IF( ABS( WR( I ) )+ABS( WI( I ) ).LT.
     $                         ABS( WR( I+1 ) )+ABS( WI( I+1 ) ) ) THEN
                              SORTED = .false.
<span class="comment">*</span><span class="comment">
</span>                              SWAP = WR( I )
                              WR( I ) = WR( I+1 )
                              WR( I+1 ) = SWAP
<span class="comment">*</span><span class="comment">
</span>                              SWAP = WI( I )
                              WI( I ) = WI( I+1 )
                              WI( I+1 ) = SWAP
                           END IF
   40                   CONTINUE
   50                CONTINUE
   60                CONTINUE
                  END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 ==== Shuffle shifts into pairs of real shifts
</span><span class="comment">*</span><span class="comment">                 .    and pairs of complex conjugate shifts
</span><span class="comment">*</span><span class="comment">                 .    assuming complex conjugate shifts are
</span><span class="comment">*</span><span class="comment">                 .    already adjacent to one another. (Yes,
</span><span class="comment">*</span><span class="comment">                 .    they are.)  ====
</span><span class="comment">*</span><span class="comment">
</span>                  DO 70 I = KBOT, KS + 2, -2
                     IF( WI( I ).NE.-WI( I-1 ) ) THEN
<span class="comment">*</span><span class="comment">
</span>                        SWAP = WR( I )
                        WR( I ) = WR( I-1 )
                        WR( I-1 ) = WR( I-2 )
                        WR( I-2 ) = SWAP
<span class="comment">*</span><span class="comment">
</span>                        SWAP = WI( I )
                        WI( I ) = WI( I-1 )
                        WI( I-1 ) = WI( I-2 )
                        WI( I-2 ) = SWAP
                     END IF
   70             CONTINUE
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== If there are only two shifts and both are
</span><span class="comment">*</span><span class="comment">              .    real, then use only one.  ====
</span><span class="comment">*</span><span class="comment">
</span>               IF( KBOT-KS+1.EQ.2 ) THEN
                  IF( WI( KBOT ).EQ.ZERO ) THEN
                     IF( ABS( WR( KBOT )-H( KBOT, KBOT ) ).LT.
     $                   ABS( WR( KBOT-1 )-H( KBOT, KBOT ) ) ) THEN
                        WR( KBOT-1 ) = WR( KBOT )
                     ELSE
                        WR( KBOT ) = WR( KBOT-1 )
                     END IF
                  END IF
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Use up to NS of the the smallest magnatiude
</span><span class="comment">*</span><span class="comment">              .    shifts.  If there aren't NS shifts available,
</span><span class="comment">*</span><span class="comment">              .    then use them all, possibly dropping one to
</span><span class="comment">*</span><span class="comment">              .    make the number of shifts even. ====
</span><span class="comment">*</span><span class="comment">
</span>               NS = MIN( NS, KBOT-KS+1 )
               NS = NS - MOD( NS, 2 )
               KS = KBOT - NS + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Small-bulge multi-shift QR sweep:
</span><span class="comment">*</span><span class="comment">              .    split workspace under the subdiagonal into
</span><span class="comment">*</span><span class="comment">              .    - a KDU-by-KDU work array U in the lower
</span><span class="comment">*</span><span class="comment">              .      left-hand-corner,
</span><span class="comment">*</span><span class="comment">              .    - a KDU-by-at-least-KDU-but-more-is-better
</span><span class="comment">*</span><span class="comment">              .      (KDU-by-NHo) horizontal work array WH along
</span><span class="comment">*</span><span class="comment">              .      the bottom edge,
</span><span class="comment">*</span><span class="comment">              .    - and an at-least-KDU-but-more-is-better-by-KDU
</span><span class="comment">*</span><span class="comment">              .      (NVE-by-KDU) vertical work WV arrow along
</span><span class="comment">*</span><span class="comment">              .      the left-hand-edge. ====
</span><span class="comment">*</span><span class="comment">
</span>               KDU = 3*NS - 3
               KU = N - KDU + 1
               KWH = KDU + 1
               NHO = ( N-KDU+1-4 ) - ( KDU+1 ) + 1
               KWV = KDU + 4
               NVE = N - KDU - KWV + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Small-bulge multi-shift QR sweep ====
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DLAQR5.610"></a><a href="dlaqr5.f.html#DLAQR5.1">DLAQR5</a>( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NS,
     $                      WR( KS ), WI( KS ), H, LDH, ILOZ, IHIZ, Z,
     $                      LDZ, WORK, 3, H( KU, 1 ), LDH, NVE,
     $                      H( KWV, 1 ), LDH, NHO, H( KU, KWH ), LDH )
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== Note progress (or the lack of it). ====
</span><span class="comment">*</span><span class="comment">
</span>            IF( LD.GT.0 ) THEN
               NDFL = 1
            ELSE
               NDFL = NDFL + 1
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== End of main loop ====
</span>   80    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Iteration limit exceeded.  Set INFO to show where
</span><span class="comment">*</span><span class="comment">        .    the problem occurred and exit. ====
</span><span class="comment">*</span><span class="comment">
</span>         INFO = KBOT
   90    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== Return the optimal value of LWORK. ====
</span><span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = DBLE( LWKOPT )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== End of <a name="DLAQR4.638"></a><a href="dlaqr4.f.html#DLAQR4.1">DLAQR4</a> ====
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>