claqr3.f.html

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

</span><span class="comment">*</span><span class="comment">          If LWORK = -1, then a workspace query is assumed; <a name="CLAQR3.140"></a><a href="claqr3.f.html#CLAQR3.1">CLAQR3</a>
</span><span class="comment">*</span><span class="comment">          only estimates the optimal workspace size for the given
</span><span class="comment">*</span><span class="comment">          values of N, NW, KTOP and KBOT.  The estimate is returned
</span><span class="comment">*</span><span class="comment">          in WORK(1).  No error message related to LWORK is issued
</span><span class="comment">*</span><span class="comment">          by <a name="XERBLA.144"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>.  Neither H nor Z are accessed.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ================================================================
</span><span class="comment">*</span><span class="comment">     Based on contributions by
</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">     .. Parameters ..
</span>      COMPLEX            ZERO, ONE
      PARAMETER          ( ZERO = ( 0.0e0, 0.0e0 ),
     $                   ONE = ( 1.0e0, 0.0e0 ) )
      REAL               RZERO, RONE
      PARAMETER          ( RZERO = 0.0e0, RONE = 1.0e0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      COMPLEX            BETA, CDUM, S, TAU
      REAL               FOO, SAFMAX, SAFMIN, SMLNUM, ULP
      INTEGER            I, IFST, ILST, INFO, INFQR, J, JW, KCOL, KLN,
     $                   KNT, KROW, KWTOP, LTOP, LWK1, LWK2, LWK3,
     $                   LWKOPT, NMIN
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      REAL               <a name="SLAMCH.167"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
      INTEGER            <a name="ILAENV.168"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
      EXTERNAL           <a name="SLAMCH.169"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="ILAENV.169"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           CCOPY, <a name="CGEHRD.172"></a><a href="cgehrd.f.html#CGEHRD.1">CGEHRD</a>, CGEMM, <a name="CLACPY.172"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>, <a name="CLAHQR.172"></a><a href="clahqr.f.html#CLAHQR.1">CLAHQR</a>, <a name="CLAQR4.172"></a><a href="claqr4.f.html#CLAQR4.1">CLAQR4</a>,
     $                   <a name="CLARF.173"></a><a href="clarf.f.html#CLARF.1">CLARF</a>, <a name="CLARFG.173"></a><a href="clarfg.f.html#CLARFG.1">CLARFG</a>, <a name="CLASET.173"></a><a href="claset.f.html#CLASET.1">CLASET</a>, <a name="CTREXC.173"></a><a href="ctrexc.f.html#CTREXC.1">CTREXC</a>, <a name="CUNGHR.173"></a><a href="cunghr.f.html#CUNGHR.1">CUNGHR</a>, <a name="SLABAD.173"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, AIMAG, CMPLX, CONJG, INT, MAX, MIN, REAL
<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><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>      JW = MIN( NW, KBOT-KTOP+1 )
      IF( JW.LE.2 ) THEN
         LWKOPT = 1
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Workspace query call to <a name="CGEHRD.193"></a><a href="cgehrd.f.html#CGEHRD.1">CGEHRD</a> ====
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CGEHRD.195"></a><a href="cgehrd.f.html#CGEHRD.1">CGEHRD</a>( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO )
         LWK1 = INT( WORK( 1 ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Workspace query call to <a name="CUNGHR.198"></a><a href="cunghr.f.html#CUNGHR.1">CUNGHR</a> ====
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CUNGHR.200"></a><a href="cunghr.f.html#CUNGHR.1">CUNGHR</a>( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO )
         LWK2 = INT( WORK( 1 ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Workspace query call to <a name="CLAQR4.203"></a><a href="claqr4.f.html#CLAQR4.1">CLAQR4</a> ====
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CLAQR4.205"></a><a href="claqr4.f.html#CLAQR4.1">CLAQR4</a>( .true., .true., JW, 1, JW, T, LDT, SH, 1, JW, V,
     $                LDV, WORK, -1, INFQR )
         LWK3 = INT( WORK( 1 ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Optimal workspace ====
</span><span class="comment">*</span><span class="comment">
</span>         LWKOPT = MAX( JW+MAX( LWK1, LWK2 ), LWK3 )
      END IF
<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">     ==== Nothing to do ...
</span><span class="comment">*</span><span class="comment">     ... for an empty active block ... ====
</span>      NS = 0
      ND = 0
      IF( KTOP.GT.KBOT )
     $   RETURN
<span class="comment">*</span><span class="comment">     ... nor for an empty deflation window. ====
</span>      IF( NW.LT.1 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== Machine constants ====
</span><span class="comment">*</span><span class="comment">
</span>      SAFMIN = <a name="SLAMCH.233"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'SAFE MINIMUM'</span> )
      SAFMAX = RONE / SAFMIN
      CALL <a name="SLABAD.235"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>( SAFMIN, SAFMAX )
      ULP = <a name="SLAMCH.236"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'PRECISION'</span> )
      SMLNUM = SAFMIN*( REAL( N ) / ULP )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== Setup deflation window ====
</span><span class="comment">*</span><span class="comment">
</span>      JW = MIN( NW, KBOT-KTOP+1 )
      KWTOP = KBOT - JW + 1
      IF( KWTOP.EQ.KTOP ) THEN
         S = ZERO
      ELSE
         S = H( KWTOP, KWTOP-1 )
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( KBOT.EQ.KWTOP ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== 1-by-1 deflation window: not much to do ====
</span><span class="comment">*</span><span class="comment">
</span>         SH( KWTOP ) = H( KWTOP, KWTOP )
         NS = 1
         ND = 0
         IF( CABS1( S ).LE.MAX( SMLNUM, ULP*CABS1( H( KWTOP,
     $       KWTOP ) ) ) ) THEN

            NS = 0
            ND = 1
            IF( KWTOP.GT.KTOP )
     $         H( KWTOP, KWTOP-1 ) = ZERO
         END IF
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== Convert to spike-triangular form.  (In case of a
</span><span class="comment">*</span><span class="comment">     .    rare QR failure, this routine continues to do
</span><span class="comment">*</span><span class="comment">     .    aggressive early deflation using that part of
</span><span class="comment">*</span><span class="comment">     .    the deflation window that converged using INFQR
</span><span class="comment">*</span><span class="comment">     .    here and there to keep track.) ====
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="CLACPY.273"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'U'</span>, JW, JW, H( KWTOP, KWTOP ), LDH, T, LDT )
      CALL CCOPY( JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ), LDT+1 )
<span class="comment">*</span><span class="comment">
</span>      CALL <a name="CLASET.276"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'A'</span>, JW, JW, ZERO, ONE, V, LDV )
      NMIN = <a name="ILAENV.277"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 12, <span class="string">'<a name="CLAQR3.277"></a><a href="claqr3.f.html#CLAQR3.1">CLAQR3</a>'</span>, <span class="string">'SV'</span>, JW, 1, JW, LWORK )
      IF( JW.GT.NMIN ) THEN
         CALL <a name="CLAQR4.279"></a><a href="claqr4.f.html#CLAQR4.1">CLAQR4</a>( .true., .true., JW, 1, JW, T, LDT, SH( KWTOP ), 1,
     $                JW, V, LDV, WORK, LWORK, INFQR )
      ELSE
         CALL <a name="CLAHQR.282"></a><a href="clahqr.f.html#CLAHQR.1">CLAHQR</a>( .true., .true., JW, 1, JW, T, LDT, SH( KWTOP ), 1,
     $                JW, V, LDV, INFQR )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== Deflation detection loop ====
</span><span class="comment">*</span><span class="comment">
</span>      NS = JW
      ILST = INFQR + 1
      DO 10 KNT = INFQR + 1, JW
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Small spike tip deflation test ====
</span><span class="comment">*</span><span class="comment">
</span>         FOO = CABS1( T( NS, NS ) )
         IF( FOO.EQ.RZERO )
     $      FOO = CABS1( S )
         IF( CABS1( S )*CABS1( V( 1, NS ) ).LE.MAX( SMLNUM, ULP*FOO ) )