claqr0.f.html

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

                  DO 30 I = KBOT, KS + 1, -2
                     W( I ) = H( I, I ) + WILK1*CABS1( H( I, I-1 ) )
                     W( I-1 ) = W( I )
   30             CONTINUE
               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="CLAQR4.459"></a><a href="claqr4.f.html#CLAQR4.1">CLAQR4</a> or
</span><span class="comment">*</span><span class="comment">                 .    <a name="CLAHQR.460"></a><a href="clahqr.f.html#CLAHQR.1">CLAHQR</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">
</span>                  IF( KBOT-KS+1.LE.NS / 2 ) THEN
                     KS = KBOT - NS + 1
                     KT = N - NS + 1
                     CALL <a name="CLACPY.468"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'A'</span>, NS, NS, H( KS, KS ), LDH,
     $                            H( KT, 1 ), LDH )
                     IF( NS.GT.NMIN ) THEN
                        CALL <a name="CLAQR4.471"></a><a href="claqr4.f.html#CLAQR4.1">CLAQR4</a>( .false., .false., NS, 1, NS,
     $                               H( KT, 1 ), LDH, W( KS ), 1, 1,
     $                               ZDUM, 1, WORK, LWORK, INF )
                     ELSE
                        CALL <a name="CLAHQR.475"></a><a href="clahqr.f.html#CLAHQR.1">CLAHQR</a>( .false., .false., NS, 1, NS,
     $                               H( KT, 1 ), LDH, W( KS ), 1, 1,
     $                               ZDUM, 1, INF )
                     END IF
                     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.  Scale to avoid
</span><span class="comment">*</span><span class="comment">                    .    overflows, underflows and subnormals.
</span><span class="comment">*</span><span class="comment">                    .    (The scale factor S can not be zero,
</span><span class="comment">*</span><span class="comment">                    .    because H(KBOT,KBOT-1) is nonzero.) ====
</span><span class="comment">*</span><span class="comment">
</span>                     IF( KS.GE.KBOT ) THEN
                        S = CABS1( H( KBOT-1, KBOT-1 ) ) +
     $                      CABS1( H( KBOT, KBOT-1 ) ) +
     $                      CABS1( H( KBOT-1, KBOT ) ) +
     $                      CABS1( H( KBOT, KBOT ) )
                        AA = H( KBOT-1, KBOT-1 ) / S
                        CC = H( KBOT, KBOT-1 ) / S
                        BB = H( KBOT-1, KBOT ) / S
                        DD = H( KBOT, KBOT ) / S
                        TR2 = ( AA+DD ) / TWO
                        DET = ( AA-TR2 )*( DD-TR2 ) - BB*CC
                        RTDISC = SQRT( -DET )
                        W( KBOT-1 ) = ( TR2+RTDISC )*S
                        W( KBOT ) = ( TR2-RTDISC )*S
<span class="comment">*</span><span class="comment">
</span>                        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">
</span>                     SORTED = .false.
                     DO 50 K = KBOT, KS + 1, -1
                        IF( SORTED )
     $                     GO TO 60
                        SORTED = .true.
                        DO 40 I = KS, K - 1
                           IF( CABS1( W( I ) ).LT.CABS1( W( I+1 ) ) )
     $                          THEN
                              SORTED = .false.
                              SWAP = W( I )
                              W( I ) = W( I+1 )
                              W( I+1 ) = SWAP
                           END IF
   40                   CONTINUE
   50                CONTINUE
   60                CONTINUE
                  END IF
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== If there are only two shifts, then use
</span><span class="comment">*</span><span class="comment">              .    only one.  ====
</span><span class="comment">*</span><span class="comment">
</span>               IF( KBOT-KS+1.EQ.2 ) THEN
                  IF( CABS1( W( KBOT )-H( KBOT, KBOT ) ).LT.
     $                CABS1( W( KBOT-1 )-H( KBOT, KBOT ) ) ) THEN
                     W( KBOT-1 ) = W( KBOT )
                  ELSE
                     W( KBOT ) = W( KBOT-1 )
                  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="CLAQR5.571"></a><a href="claqr5.f.html#CLAQR5.1">CLAQR5</a>( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NS,
     $                      W( 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>   70    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
   80    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 ) = CMPLX( LWKOPT, 0 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== End of <a name="CLAQR0.599"></a><a href="claqr0.f.html#CLAQR0.1">CLAQR0</a> ====
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>