dsvdct.f

来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 162 行

      SUBROUTINE DSVDCT( N, S, E, SHIFT, NUM )
*
*  -- LAPACK test routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      INTEGER            N, NUM
      DOUBLE PRECISION   SHIFT
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   E( * ), S( * )
*     ..
*
*  Purpose
*  =======
*
*  DSVDCT counts the number NUM of eigenvalues of a 2*N by 2*N
*  tridiagonal matrix T which are less than or equal to SHIFT.  T is
*  formed by putting zeros on the diagonal and making the off-diagonals
*  equal to S(1), E(1), S(2), E(2), ... , E(N-1), S(N).  If SHIFT is
*  positive, NUM is equal to N plus the number of singular values of a
*  bidiagonal matrix B less than or equal to SHIFT.  Here B has diagonal
*  entries S(1), ..., S(N) and superdiagonal entries E(1), ... E(N-1).
*  If SHIFT is negative, NUM is equal to the number of singular values
*  of B greater than or equal to -SHIFT.
*
*  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
*  Matrix", Report CS41, Computer Science Dept., Stanford University,
*  July 21, 1966
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The dimension of the bidiagonal matrix B.
*
*  S       (input) DOUBLE PRECISION array, dimension (N)
*          The diagonal entries of the bidiagonal matrix B.
*
*  E       (input) DOUBLE PRECISION array of dimension (N-1)
*          The superdiagonal entries of the bidiagonal matrix B.
*
*  SHIFT   (input) DOUBLE PRECISION
*          The shift, used as described under Purpose.
*
*  NUM     (output) INTEGER
*          The number of eigenvalues of T less than or equal to SHIFT.
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ONE
      PARAMETER          ( ONE = 1.0D0 )
      DOUBLE PRECISION   ZERO
      PARAMETER          ( ZERO = 0.0D0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I
      DOUBLE PRECISION   M1, M2, MX, OVFL, SOV, SSHIFT, SSUN, SUN, TMP,
     $                   TOM, U, UNFL
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           DLAMCH
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, MAX, SQRT
*     ..
*     .. Executable Statements ..
*
*     Get machine constants
*
      UNFL = 2*DLAMCH( 'Safe minimum' )
      OVFL = ONE / UNFL
*
*     Find largest entry
*
      MX = ABS( S( 1 ) )
      DO 10 I = 1, N - 1
         MX = MAX( MX, ABS( S( I+1 ) ), ABS( E( I ) ) )
   10 CONTINUE
*
      IF( MX.EQ.ZERO ) THEN
         IF( SHIFT.LT.ZERO ) THEN
            NUM = 0
         ELSE
            NUM = 2*N
         END IF
         RETURN
      END IF
*
*     Compute scale factors as in Kahan's report
*
      SUN = SQRT( UNFL )
      SSUN = SQRT( SUN )
      SOV = SQRT( OVFL )
      TOM = SSUN*SOV
      IF( MX.LE.ONE ) THEN
         M1 = ONE / MX
         M2 = TOM
      ELSE
         M1 = ONE
         M2 = TOM / MX
      END IF
*
*     Begin counting
*
      U = ONE
      NUM = 0
      SSHIFT = ( SHIFT*M1 )*M2
      U = -SSHIFT
      IF( U.LE.SUN ) THEN
         IF( U.LE.ZERO ) THEN
            NUM = NUM + 1
            IF( U.GT.-SUN )
     $         U = -SUN
         ELSE
            U = SUN
         END IF
      END IF
      TMP = ( S( 1 )*M1 )*M2
      U = -TMP*( TMP / U ) - SSHIFT
      IF( U.LE.SUN ) THEN
         IF( U.LE.ZERO ) THEN
            NUM = NUM + 1
            IF( U.GT.-SUN )
     $         U = -SUN
         ELSE
            U = SUN
         END IF
      END IF
      DO 20 I = 1, N - 1
         TMP = ( E( I )*M1 )*M2
         U = -TMP*( TMP / U ) - SSHIFT
         IF( U.LE.SUN ) THEN
            IF( U.LE.ZERO ) THEN
               NUM = NUM + 1
               IF( U.GT.-SUN )
     $            U = -SUN
            ELSE
               U = SUN
            END IF
         END IF
         TMP = ( S( I+1 )*M1 )*M2
         U = -TMP*( TMP / U ) - SSHIFT
         IF( U.LE.SUN ) THEN
            IF( U.LE.ZERO ) THEN
               NUM = NUM + 1
               IF( U.GT.-SUN )
     $            U = -SUN
            ELSE
               U = SUN
            END IF
         END IF
   20 CONTINUE
      RETURN
*
*     End of DSVDCT
*
      END