claqr0.f.html

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

</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="CLAHQR.157"></a><a href="clahqr.f.html#CLAHQR.1">CLAHQR</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
      PARAMETER          ( WILK1 = 0.75e0 )
      COMPLEX            ZERO, ONE
      PARAMETER          ( ZERO = ( 0.0e0, 0.0e0 ),
     $                   ONE = ( 1.0e0, 0.0e0 ) )
      REAL               TWO
      PARAMETER          ( TWO = 2.0e0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      COMPLEX            AA, BB, CC, CDUM, DD, DET, RTDISC, SWAP, TR2
      REAL               S
      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.192"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
      EXTERNAL           <a name="ILAENV.193"></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>      COMPLEX            ZDUM( 1, 1 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="CLACPY.199"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>, <a name="CLAHQR.199"></a><a href="clahqr.f.html#CLAHQR.1">CLAHQR</a>, <a name="CLAQR3.199"></a><a href="claqr3.f.html#CLAQR3.1">CLAQR3</a>, <a name="CLAQR4.199"></a><a href="claqr4.f.html#CLAQR4.1">CLAQR4</a>, <a name="CLAQR5.199"></a><a href="claqr5.f.html#CLAQR5.1">CLAQR5</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, AIMAG, CMPLX, INT, MAX, MIN, MOD, REAL,
     $                   SQRT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Functions ..
</span>      REAL               CABS1
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Function definitions ..
</span>      CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
<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.221"></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="CLAHQR.234"></a><a href="clahqr.f.html#CLAHQR.1">CLAHQR</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="CLAHQR.242"></a><a href="clahqr.f.html#CLAHQR.1">CLAHQR</a>( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, 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.259"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 13, <span class="string">'<a name="CLAQR0.259"></a><a href="claqr0.f.html#CLAQR0.1">CLAQR0</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.269"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 15, <span class="string">'<a name="CLAQR0.269"></a><a href="claqr0.f.html#CLAQR0.1">CLAQR0</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="CLAQR3.275"></a><a href="claqr3.f.html#CLAQR3.1">CLAQR3</a> ====
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CLAQR3.277"></a><a href="claqr3.f.html#CLAQR3.1">CLAQR3</a>( WANTT, WANTZ, N, ILO, IHI, NWR+1, H, LDH, ILOZ,
     $                IHIZ, Z, LDZ, LS, LD, W, 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="CLAQR5.281"></a><a href="claqr5.f.html#CLAQR5.1">CLAQR5</a>, <a name="CLAQR3.281"></a><a href="claqr3.f.html#CLAQR3.1">CLAQR3</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 ) = CMPLX( LWKOPT, 0 )
            RETURN
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== <a name="CLAHQR.292"></a><a href="clahqr.f.html#CLAHQR.1">CLAHQR</a>/<a name="CLAQR0.292"></a><a href="claqr0.f.html#CLAQR0.1">CLAQR0</a> crossover point ====
</span><span class="comment">*</span><span class="comment">
</span>         NMIN = <a name="ILAENV.294"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 12, <span class="string">'<a name="CLAQR0.294"></a><a href="claqr0.f.html#CLAQR0.1">CLAQR0</a>'</span>, JBCMPZ, N, ILO, IHI, LWORK )
         NMIN = MAX( NTINY, NMIN )