clahqr.f.html

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

</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           CCOPY, <a name="CLARFG.148"></a><a href="clarfg.f.html#CLARFG.1">CLARFG</a>, CSCAL, <a name="SLABAD.148"></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">     .. Statement Functions ..
</span>      REAL               CABS1
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, AIMAG, CONJG, MAX, MIN, REAL, SQRT
<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>      INFO = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.0 )
     $   RETURN
      IF( ILO.EQ.IHI ) THEN
         W( ILO ) = H( ILO, ILO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== clear out the trash ====
</span>      DO 10 J = ILO, IHI - 3
         H( J+2, J ) = ZERO
         H( J+3, J ) = ZERO
   10 CONTINUE
      IF( ILO.LE.IHI-2 )
     $   H( IHI, IHI-2 ) = ZERO
<span class="comment">*</span><span class="comment">     ==== ensure that subdiagonal entries are real ====
</span>      DO 20 I = ILO + 1, IHI
         IF( AIMAG( H( I, I-1 ) ).NE.RZERO ) THEN
<span class="comment">*</span><span class="comment">           ==== The following redundant normalization
</span><span class="comment">*</span><span class="comment">           .    avoids problems with both gradual and
</span><span class="comment">*</span><span class="comment">           .    sudden underflow in ABS(H(I,I-1)) ====
</span>            SC = H( I, I-1 ) / CABS1( H( I, I-1 ) )
            SC = CONJG( SC ) / ABS( SC )
            H( I, I-1 ) = ABS( H( I, I-1 ) )
            IF( WANTT ) THEN
               JLO = 1
               JHI = N
            ELSE
               JLO = ILO
               JHI = IHI
            END IF
            CALL CSCAL( JHI-I+1, SC, H( I, I ), LDH )
            CALL CSCAL( MIN( JHI, I+1 )-JLO+1, CONJG( SC ), H( JLO, I ),
     $                  1 )
            IF( WANTZ )
     $         CALL CSCAL( IHIZ-ILOZ+1, CONJG( SC ), Z( ILOZ, I ), 1 )
         END IF
   20 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      NH = IHI - ILO + 1
      NZ = IHIZ - ILOZ + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set machine-dependent constants for the stopping criterion.
</span><span class="comment">*</span><span class="comment">
</span>      SAFMIN = <a name="SLAMCH.208"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'SAFE MINIMUM'</span> )
      SAFMAX = RONE / SAFMIN
      CALL <a name="SLABAD.210"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>( SAFMIN, SAFMAX )
      ULP = <a name="SLAMCH.211"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'PRECISION'</span> )
      SMLNUM = SAFMIN*( REAL( NH ) / ULP )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     I1 and I2 are the indices of the first row and last column of H
</span><span class="comment">*</span><span class="comment">     to which transformations must be applied. If eigenvalues only are
</span><span class="comment">*</span><span class="comment">     being computed, I1 and I2 are set inside the main loop.
</span><span class="comment">*</span><span class="comment">
</span>      IF( WANTT ) THEN
         I1 = 1
         I2 = N
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     The main loop begins here. I is the loop index and decreases from
</span><span class="comment">*</span><span class="comment">     IHI to ILO in steps of 1. Each iteration of the loop works
</span><span class="comment">*</span><span class="comment">     with the active submatrix in rows and columns L to I.
</span><span class="comment">*</span><span class="comment">     Eigenvalues I+1 to IHI have already converged. Either L = ILO, or
</span><span class="comment">*</span><span class="comment">     H(L,L-1) is negligible so that the matrix splits.
</span><span class="comment">*</span><span class="comment">
</span>      I = IHI
   30 CONTINUE
      IF( I.LT.ILO )
     $   GO TO 150
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Perform QR iterations on rows and columns ILO to I until a
</span><span class="comment">*</span><span class="comment">     submatrix of order 1 splits off at the bottom because a
</span><span class="comment">*</span><span class="comment">     subdiagonal element has become negligible.
</span><span class="comment">*</span><span class="comment">
</span>      L = ILO
      DO 130 ITS = 0, ITMAX
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Look for a single small subdiagonal element.
</span><span class="comment">*</span><span class="comment">
</span>         DO 40 K = I, L + 1, -1
            IF( CABS1( H( K, K-1 ) ).LE.SMLNUM )
     $         GO TO 50
            TST = CABS1( H( K-1, K-1 ) ) + CABS1( H( K, K ) )
            IF( TST.EQ.ZERO ) THEN
               IF( K-2.GE.ILO )
     $            TST = TST + ABS( REAL( H( K-1, K-2 ) ) )
               IF( K+1.LE.IHI )
     $            TST = TST + ABS( REAL( H( K+1, K ) ) )
            END IF
<span class="comment">*</span><span class="comment">           ==== The following is a conservative small subdiagonal
</span><span class="comment">*</span><span class="comment">           .    deflation criterion due to Ahues &amp; Tisseur (LAWN 122,
</span><span class="comment">*</span><span class="comment">           .    1997). It has better mathematical foundation and
</span><span class="comment">*</span><span class="comment">           .    improves accuracy in some examples.  ====
</span>            IF( ABS( REAL( H( K, K-1 ) ) ).LE.ULP*TST ) THEN
               AB = MAX( CABS1( H( K, K-1 ) ), CABS1( H( K-1, K ) ) )
               BA = MIN( CABS1( H( K, K-1 ) ), CABS1( H( K-1, K ) ) )
               AA = MAX( CABS1( H( K, K ) ),
     $              CABS1( H( K-1, K-1 )-H( K, K ) ) )
               BB = MIN( CABS1( H( K, K ) ),
     $              CABS1( H( K-1, K-1 )-H( K, K ) ) )
               S = AA + AB
               IF( BA*( AB / S ).LE.MAX( SMLNUM,
     $             ULP*( BB*( AA / S ) ) ) )GO TO 50
            END IF
   40    CONTINUE
   50    CONTINUE
         L = K
         IF( L.GT.ILO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           H(L,L-1) is negligible
</span><span class="comment">*</span><span class="comment">
</span>            H( L, L-1 ) = ZERO
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Exit from loop if a submatrix of order 1 has split off.
</span><span class="comment">*</span><span class="comment">
</span>         IF( L.GE.I )
     $      GO TO 140
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Now the active submatrix is in rows and columns L to I. If
</span><span class="comment">*</span><span class="comment">        eigenvalues only are being computed, only the active submatrix
</span><span class="comment">*</span><span class="comment">        need be transformed.
</span><span class="comment">*</span><span class="comment">
</span>         IF( .NOT.WANTT ) THEN
            I1 = L
            I2 = I
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( ITS.EQ.10 .OR. ITS.EQ.20 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Exceptional shift.
</span><span class="comment">*</span><span class="comment">
</span>            S = DAT1*ABS( REAL( H( I, I-1 ) ) )
            T = S + H( I, I )
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Wilkinson's shift.
</span><span class="comment">*</span><span class="comment">
</span>            T = H( I, I )
            U = SQRT( H( I-1, I ) )*SQRT( H( I, I-1 ) )
            S = CABS1( U )
            IF( S.NE.RZERO ) THEN
               X = HALF*( H( I-1, I-1 )-T )
               SX = CABS1( X )
               S = MAX( S, CABS1( X ) )
               Y = S*SQRT( ( X / S )**2+( U / S )**2 )
               IF( SX.GT.RZERO ) THEN
                  IF( REAL( X / SX )*REAL( Y )+AIMAG( X / SX )*