zhseqr.f

计算矩阵的经典开源库．全世界都在用它．相信你也不能例外．

     $                LDZ, INFO )
         RETURN
      END IF
      MAXB = MAX( 2, MAXB )
      NS = MIN( NS, MAXB, NSMAX )
*
*     Now 1 < NS <= MAXB < NH.
*
*     Set machine-dependent constants for the stopping criterion.
*     If norm(H) <= sqrt(OVFL), overflow should not occur.
*
      UNFL = DLAMCH( 'Safe minimum' )
      OVFL = RONE / UNFL
      CALL DLABAD( UNFL, OVFL )
      ULP = DLAMCH( 'Precision' )
      SMLNUM = UNFL*( NH / ULP )
*
*     ITN is the total number of multiple-shift QR iterations allowed.
*
      ITN = 30*NH
*
*     The main loop begins here. I is the loop index and decreases from
*     IHI to ILO in steps of at most MAXB. Each iteration of the loop
*     works with the active submatrix in rows and columns L to I.
*     Eigenvalues I+1 to IHI have already converged. Either L = ILO, or
*     H(L,L-1) is negligible so that the matrix splits.
*
      I = IHI
   60 CONTINUE
      IF( I.LT.ILO )
     $   GO TO 180
*
*     Perform multiple-shift QR iterations on rows and columns ILO to I
*     until a submatrix of order at most MAXB splits off at the bottom
*     because a subdiagonal element has become negligible.
*
      L = ILO
      DO 160 ITS = 0, ITN
*
*        Look for a single small subdiagonal element.
*
         DO 70 K = I, L + 1, -1
            TST1 = CABS1( H( K-1, K-1 ) ) + CABS1( H( K, K ) )
            IF( TST1.EQ.RZERO ) THEN
               TST1 = ZLANHS( '1', I-L+1, H( L, L ), LDH, RWORK )
***
*              Increment op count
               OPS = OPS + 5*( I-L+1 )*( I-L ) / 2
***
            END IF
            IF( ABS( DBLE( H( K, K-1 ) ) ).LE.MAX( ULP*TST1, SMLNUM ) )
     $         GO TO 80
   70    CONTINUE
   80    CONTINUE
         L = K
***
*        Increment op count
         OPST = OPST + 5*( I-L+1 )
***
         IF( L.GT.ILO ) THEN
*
*           H(L,L-1) is negligible.
*
            H( L, L-1 ) = ZERO
         END IF
*
*        Exit from loop if a submatrix of order <= MAXB has split off.
*
         IF( L.GE.I-MAXB+1 )
     $      GO TO 170
*
*        Now the active submatrix is in rows and columns L to I. If
*        eigenvalues only are being computed, only the active submatrix
*        need be transformed.
*
         IF( .NOT.WANTT ) THEN
            I1 = L
            I2 = I
         END IF
*
         IF( ITS.EQ.20 .OR. ITS.EQ.30 ) THEN
*
*           Exceptional shifts.
*
            DO 90 II = I - NS + 1, I
               W( II ) = CONST*( ABS( DBLE( H( II, II-1 ) ) )+
     $                   ABS( DBLE( H( II, II ) ) ) )
   90       CONTINUE
***
*           Increment op count
            OPST = OPST + 2*NS
***
         ELSE
*
*           Use eigenvalues of trailing submatrix of order NS as shifts.
*
            CALL ZLACPY( 'Full', NS, NS, H( I-NS+1, I-NS+1 ), LDH, S,
     $                   LDS )
            CALL ZLAHQR( .FALSE., .FALSE., NS, 1, NS, S, LDS,
     $                   W( I-NS+1 ), 1, NS, Z, LDZ, IERR )
            IF( IERR.GT.0 ) THEN
*
*              If ZLAHQR failed to compute all NS eigenvalues, use the
*              unconverged diagonal elements as the remaining shifts.
*
               DO 100 II = 1, IERR
                  W( I-NS+II ) = S( II, II )
  100          CONTINUE
            END IF
         END IF
*
*        Form the first column of (G-w(1)) (G-w(2)) . . . (G-w(ns))
*        where G is the Hessenberg submatrix H(L:I,L:I) and w is
*        the vector of shifts (stored in W). The result is
*        stored in the local array V.
*
         V( 1 ) = ONE
         DO 110 II = 2, NS + 1
            V( II ) = ZERO
  110    CONTINUE
         NV = 1
         DO 130 J = I - NS + 1, I
            CALL ZCOPY( NV+1, V, 1, VV, 1 )
            CALL ZGEMV( 'No transpose', NV+1, NV, ONE, H( L, L ), LDH,
     $                  VV, 1, -W( J ), V, 1 )
            NV = NV + 1
***
*           Increment op count
            OPST = OPST + 8*NV*( N+1 ) + 6*( NV+1 )
***
*
*           Scale V(1:NV) so that max(abs(V(i))) = 1. If V is zero,
*           reset it to the unit vector.
*
            ITEMP = IZAMAX( NV, V, 1 )
***
*           Increment op count
            OPST = OPST + 2*NV
***
            RTEMP = CABS1( V( ITEMP ) )
            IF( RTEMP.EQ.RZERO ) THEN
               V( 1 ) = ONE
               DO 120 II = 2, NV
                  V( II ) = ZERO
  120          CONTINUE
            ELSE
               RTEMP = MAX( RTEMP, SMLNUM )
               CALL ZDSCAL( NV, RONE / RTEMP, V, 1 )
***
*              Increment op count
               OPST = OPST + 2*NV
***
            END IF
  130    CONTINUE
*
*        Multiple-shift QR step
*
         DO 150 K = L, I - 1
*
*           The first iteration of this loop determines a reflection G
*           from the vector V and applies it from left and right to H,
*           thus creating a nonzero bulge below the subdiagonal.
*
*           Each subsequent iteration determines a reflection G to
*           restore the Hessenberg form in the (K-1)th column, and thus
*           chases the bulge one step toward the bottom of the active
*           submatrix. NR is the order of G.
*
            NR = MIN( NS+1, I-K+1 )
            IF( K.GT.L )
     $         CALL ZCOPY( NR, H( K, K-1 ), 1, V, 1 )
            CALL ZLARFG( NR, V( 1 ), V( 2 ), 1, TAU )
***
*           Increment op count
            OPST = OPST + 10*NR + 12
***
            IF( K.GT.L ) THEN
               H( K, K-1 ) = V( 1 )
               DO 140 II = K + 1, I
                  H( II, K-1 ) = ZERO
  140          CONTINUE
            END IF
            V( 1 ) = ONE
*
*           Apply G' from the left to transform the rows of the matrix
*           in columns K to I2.
*
            CALL ZLARFX( 'Left', NR, I2-K+1, V, DCONJG( TAU ),
     $                   H( K, K ), LDH, WORK )
*
*           Apply G from the right to transform the columns of the
*           matrix in rows I1 to min(K+NR,I).
*
            CALL ZLARFX( 'Right', MIN( K+NR, I )-I1+1, NR, V, TAU,
     $                   H( I1, K ), LDH, WORK )
***
*           Increment op count
            OPS = OPS + 4*( 4*NR-2 )*( I2-I1+2+MIN( NR, I-K ) )
***
*
            IF( WANTZ ) THEN
*
*              Accumulate transformations in the matrix Z
*
               CALL ZLARFX( 'Right', NH, NR, V, TAU, Z( ILO, K ), LDZ,
     $                      WORK )
***
*              Increment op count
               OPS = OPS + 4*( 4*NR-2 )*NH
***
            END IF
  150    CONTINUE
*
*        Ensure that H(I,I-1) is real.
*
         TEMP = H( I, I-1 )
         IF( DIMAG( TEMP ).NE.RZERO ) THEN
            RTEMP = DLAPY2( DBLE( TEMP ), DIMAG( TEMP ) )
            H( I, I-1 ) = RTEMP
            TEMP = TEMP / RTEMP
            IF( I2.GT.I )
     $         CALL ZSCAL( I2-I, DCONJG( TEMP ), H( I, I+1 ), LDH )
            CALL ZSCAL( I-I1, TEMP, H( I1, I ), 1 )
***
*           Increment op count
            OPST = OPST + 6*( I2-I1+1 )
***
            IF( WANTZ ) THEN
               CALL ZSCAL( NH, TEMP, Z( ILO, I ), 1 )
***
*              Increment op count
               OPST = OPST + 6*NH
***
            END IF
         END IF
*
  160 CONTINUE
*
*     Failure to converge in remaining number of iterations
*
      INFO = I
      RETURN
*
  170 CONTINUE
*
*     A submatrix of order <= MAXB in rows and columns L to I has split
*     off. Use the double-shift QR algorithm to handle it.
*
      CALL ZLAHQR( WANTT, WANTZ, N, L, I, H, LDH, W, ILO, IHI, Z, LDZ,
     $             INFO )
      IF( INFO.GT.0 )
     $   RETURN
*
*     Decrement number of remaining iterations, and return to start of
*     the main loop with a new value of I.
*
      ITN = ITN - ITS
      I = L - 1
      GO TO 60
*
  180 CONTINUE
***
*     Compute final op count
      OPS = OPS + OPST
***
      WORK( 1 ) = MAX( 1, N )
      RETURN
*
*     End of ZHSEQR
*
      END