dhgeqz.f.html

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

            ELSE
               H( ILAST, ILAST ) = -H( ILAST, ILAST )
               T( ILAST, ILAST ) = -T( ILAST, ILAST )
            END IF
            IF( ILZ ) THEN
               DO 100 J = 1, N
                  Z( J, ILAST ) = -Z( J, ILAST )
  100          CONTINUE
            END IF
         END IF
         ALPHAR( ILAST ) = H( ILAST, ILAST )
         ALPHAI( ILAST ) = ZERO
         BETA( ILAST ) = T( ILAST, ILAST )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Go to next block -- exit if finished.
</span><span class="comment">*</span><span class="comment">
</span>         ILAST = ILAST - 1
         IF( ILAST.LT.ILO )
     $      GO TO 380
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Reset counters
</span><span class="comment">*</span><span class="comment">
</span>         IITER = 0
         ESHIFT = ZERO
         IF( .NOT.ILSCHR ) THEN
            ILASTM = ILAST
            IF( IFRSTM.GT.ILAST )
     $         IFRSTM = ILO
         END IF
         GO TO 350
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        QZ step
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        This iteration only involves rows/columns IFIRST:ILAST. We
</span><span class="comment">*</span><span class="comment">        assume IFIRST &lt; ILAST, and that the diagonal of B is non-zero.
</span><span class="comment">*</span><span class="comment">
</span>  110    CONTINUE
         IITER = IITER + 1
         IF( .NOT.ILSCHR ) THEN
            IFRSTM = IFIRST
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute single shifts.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        At this point, IFIRST &lt; ILAST, and the diagonal elements of
</span><span class="comment">*</span><span class="comment">        T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
</span><span class="comment">*</span><span class="comment">        magnitude)
</span><span class="comment">*</span><span class="comment">
</span>         IF( ( IITER / 10 )*10.EQ.IITER ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Exceptional shift.  Chosen for no particularly good reason.
</span><span class="comment">*</span><span class="comment">           (Single shift only.)
</span><span class="comment">*</span><span class="comment">
</span>            IF( ( DBLE( MAXIT )*SAFMIN )*ABS( H( ILAST-1, ILAST ) ).LT.
     $          ABS( T( ILAST-1, ILAST-1 ) ) ) THEN
               ESHIFT = ESHIFT + H( ILAST-1, ILAST ) /
     $                  T( ILAST-1, ILAST-1 )
            ELSE
               ESHIFT = ESHIFT + ONE / ( SAFMIN*DBLE( MAXIT ) )
            END IF
            S1 = ONE
            WR = ESHIFT
<span class="comment">*</span><span class="comment">
</span>         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Shifts based on the generalized eigenvalues of the
</span><span class="comment">*</span><span class="comment">           bottom-right 2x2 block of A and B. The first eigenvalue
</span><span class="comment">*</span><span class="comment">           returned by <a name="DLAG2.643"></a><a href="dlag2.f.html#DLAG2.1">DLAG2</a> is the Wilkinson shift (AEP p.512),
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DLAG2.645"></a><a href="dlag2.f.html#DLAG2.1">DLAG2</a>( H( ILAST-1, ILAST-1 ), LDH,
     $                  T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1,
     $                  S2, WR, WR2, WI )
<span class="comment">*</span><span class="comment">
</span>            TEMP = MAX( S1, SAFMIN*MAX( ONE, ABS( WR ), ABS( WI ) ) )
            IF( WI.NE.ZERO )
     $         GO TO 200
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Fiddle with shift to avoid overflow
</span><span class="comment">*</span><span class="comment">
</span>         TEMP = MIN( ASCALE, ONE )*( HALF*SAFMAX )
         IF( S1.GT.TEMP ) THEN
            SCALE = TEMP / S1
         ELSE
            SCALE = ONE
         END IF
<span class="comment">*</span><span class="comment">
</span>         TEMP = MIN( BSCALE, ONE )*( HALF*SAFMAX )
         IF( ABS( WR ).GT.TEMP )
     $      SCALE = MIN( SCALE, TEMP / ABS( WR ) )
         S1 = SCALE*S1
         WR = SCALE*WR
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Now check for two consecutive small subdiagonals.
</span><span class="comment">*</span><span class="comment">
</span>         DO 120 J = ILAST - 1, IFIRST + 1, -1
            ISTART = J
            TEMP = ABS( S1*H( J, J-1 ) )
            TEMP2 = ABS( S1*H( J, J )-WR*T( J, J ) )
            TEMPR = MAX( TEMP, TEMP2 )
            IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
               TEMP = TEMP / TEMPR
               TEMP2 = TEMP2 / TEMPR
            END IF
            IF( ABS( ( ASCALE*H( J+1, J ) )*TEMP ).LE.( ASCALE*ATOL )*
     $          TEMP2 )GO TO 130
  120    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         ISTART = IFIRST
  130    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Do an implicit single-shift QZ sweep.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Initial Q
</span><span class="comment">*</span><span class="comment">
</span>         TEMP = S1*H( ISTART, ISTART ) - WR*T( ISTART, ISTART )
         TEMP2 = S1*H( ISTART+1, ISTART )
         CALL <a name="DLARTG.693"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>( TEMP, TEMP2, C, S, TEMPR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Sweep
</span><span class="comment">*</span><span class="comment">
</span>         DO 190 J = ISTART, ILAST - 1
            IF( J.GT.ISTART ) THEN
               TEMP = H( J, J-1 )
               CALL <a name="DLARTG.700"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>( TEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) )
               H( J+1, J-1 ) = ZERO
            END IF
<span class="comment">*</span><span class="comment">
</span>            DO 140 JC = J, ILASTM
               TEMP = C*H( J, JC ) + S*H( J+1, JC )
               H( J+1, JC ) = -S*H( J, JC ) + C*H( J+1, JC )
               H( J, JC ) = TEMP
               TEMP2 = C*T( J, JC ) + S*T( J+1, JC )
               T( J+1, JC ) = -S*T( J, JC ) + C*T( J+1, JC )
               T( J, JC ) = TEMP2
  140       CONTINUE
            IF( ILQ ) THEN
               DO 150 JR = 1, N
                  TEMP = C*Q( JR, J ) + S*Q( JR, J+1 )
                  Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
                  Q( JR, J ) = TEMP
  150          CONTINUE
            END IF
<span class="comment">*</span><span class="comment">
</span>            TEMP = T( J+1, J+1 )
            CALL <a name="DLARTG.721"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>( TEMP, T( J+1, J ), C, S, T( J+1, J+1 ) )
            T( J+1, J ) = ZERO
<span class="comment">*</span><span class="comment">
</span>            DO 160 JR = IFRSTM, MIN( J+2, ILAST )
               TEMP = C*H( JR, J+1 ) + S*H( JR, J )
               H( JR, J ) = -S*H( JR, J+1 ) + C*H( JR, J )
               H( JR, J+1 ) = TEMP
  160       CONTINUE
            DO 170 JR = IFRSTM, J
               TEMP = C*T( JR, J+1 ) + S*T( JR, J )
               T( JR, J ) = -S*T( JR, J+1 ) + C*T( JR, J )
               T( JR, J+1 ) = TEMP
  170       CONTINUE
            IF( ILZ ) THEN
               DO 180 JR = 1, N
                  TEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
                  Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J )
                  Z( JR, J+1 ) = TEMP
  180          CONTINUE
            END IF
  190    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         GO TO 350
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Use Francis double-shift
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Note: the Francis double-shift should work with real shifts,
</span><span class="comment">*</span><span class="comment">              but only if the block is at least 3x3.
</span><span class="comment">*</span><span class="comment">              This code may break if this point is reached with
</span><span class="comment">*</span><span class="comment">              a 2x2 block with real eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span>  200    CONTINUE
         IF( IFIRST+1.EQ.ILAST ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Special case -- 2x2 block with complex eigenvectors
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Step 1: Standardize, that is, rotate so that
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                       ( B11  0  )
</span><span class="comment">*</span><span class="comment">                   B = (         )  with B11 non-negative.
</span><span class="comment">*</span><span class="comment">                       (  0  B22 )
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DLASV2.763"></a><a href="dlasv2.f.html#DLASV2.1">DLASV2</a>( T( ILAST-1, ILAST-1 ), T( ILAST-1, ILAST ),
     $                   T( ILAST, ILAST ), B22, B11, SR, CR, SL, CL )
<span class="comment">*</span><span class="comment">
</span>            IF( B11.LT.ZERO ) THEN
               CR = -CR
               SR = -SR
               B11 = -B11
               B22 = -B22
            END IF
<span class="comment">*</span><span class="comment">
</span>            CALL DROT( ILASTM+1-IFIRST, H( ILAST-1, ILAST-1 ), LDH,