zlaqr2.f.html

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

</span><span class="comment">*</span><span class="comment">          On exit, WORK(1) is set to an estimate of the optimal value
</span><span class="comment">*</span><span class="comment">          of LWORK for the given values of N, NW, KTOP and KBOT.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     LWORK   (input) integer
</span><span class="comment">*</span><span class="comment">          The dimension of the work array WORK.  LWORK = 2*NW
</span><span class="comment">*</span><span class="comment">          suffices, but greater efficiency may result from larger
</span><span class="comment">*</span><span class="comment">          values of LWORK.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          If LWORK = -1, then a workspace query is assumed; <a name="ZLAQR2.144"></a><a href="zlaqr2.f.html#ZLAQR2.1">ZLAQR2</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.148"></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*16         ZERO, ONE
      PARAMETER          ( ZERO = ( 0.0d0, 0.0d0 ),
     $                   ONE = ( 1.0d0, 0.0d0 ) )
      DOUBLE PRECISION   RZERO, RONE
      PARAMETER          ( RZERO = 0.0d0, RONE = 1.0d0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      COMPLEX*16         BETA, CDUM, S, TAU
      DOUBLE PRECISION   FOO, SAFMAX, SAFMIN, SMLNUM, ULP
      INTEGER            I, IFST, ILST, INFO, INFQR, J, JW, KCOL, KLN,
     $                   KNT, KROW, KWTOP, LTOP, LWK1, LWK2, LWKOPT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      DOUBLE PRECISION   <a name="DLAMCH.170"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
      EXTERNAL           <a name="DLAMCH.171"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="DLABAD.174"></a><a href="dlabad.f.html#DLABAD.1">DLABAD</a>, ZCOPY, <a name="ZGEHRD.174"></a><a href="zgehrd.f.html#ZGEHRD.1">ZGEHRD</a>, ZGEMM, <a name="ZLACPY.174"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>, <a name="ZLAHQR.174"></a><a href="zlahqr.f.html#ZLAHQR.1">ZLAHQR</a>,
     $                   <a name="ZLARF.175"></a><a href="zlarf.f.html#ZLARF.1">ZLARF</a>, <a name="ZLARFG.175"></a><a href="zlarfg.f.html#ZLARFG.1">ZLARFG</a>, <a name="ZLASET.175"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>, <a name="ZTREXC.175"></a><a href="ztrexc.f.html#ZTREXC.1">ZTREXC</a>, <a name="ZUNGHR.175"></a><a href="zunghr.f.html#ZUNGHR.1">ZUNGHR</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, MAX, MIN
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Functions ..
</span>      DOUBLE PRECISION   CABS1
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Function definitions ..
</span>      CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( 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="ZGEHRD.195"></a><a href="zgehrd.f.html#ZGEHRD.1">ZGEHRD</a> ====
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZGEHRD.197"></a><a href="zgehrd.f.html#ZGEHRD.1">ZGEHRD</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="ZUNGHR.200"></a><a href="zunghr.f.html#ZUNGHR.1">ZUNGHR</a> ====
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZUNGHR.202"></a><a href="zunghr.f.html#ZUNGHR.1">ZUNGHR</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">        ==== Optimal workspace ====
</span><span class="comment">*</span><span class="comment">
</span>         LWKOPT = JW + MAX( LWK1, LWK2 )
      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 ) = DCMPLX( 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="DLAMCH.229"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'SAFE MINIMUM'</span> )
      SAFMAX = RONE / SAFMIN
      CALL <a name="DLABAD.231"></a><a href="dlabad.f.html#DLABAD.1">DLABAD</a>( SAFMIN, SAFMAX )
      ULP = <a name="DLAMCH.232"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'PRECISION'</span> )
      SMLNUM = SAFMIN*( DBLE( 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="ZLACPY.268"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'U'</span>, JW, JW, H( KWTOP, KWTOP ), LDH, T, LDT )
      CALL ZCOPY( JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ), LDT+1 )
<span class="comment">*</span><span class="comment">
</span>      CALL <a name="ZLASET.271"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>( <span class="string">'A'</span>, JW, JW, ZERO, ONE, V, LDV )
      CALL <a name="ZLAHQR.272"></a><a href="zlahqr.f.html#ZLAHQR.1">ZLAHQR</a>( .true., .true., JW, 1, JW, T, LDT, SH( KWTOP ), 1,
     $             JW, V, LDV, INFQR )
<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 ) )
     $        THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== One more converged eigenvalue ====