dstect.f

Go to the documentation of this file.

*> \brief \b DSTECT
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*  Definition:
*  ===========
*
*       SUBROUTINE DSTECT( N, A, B, SHIFT, NUM )
* 
*       .. Scalar Arguments ..
*       INTEGER            N, NUM
*       DOUBLE PRECISION   SHIFT
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION   A( * ), B( * )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*>    DSTECT counts the number NUM of eigenvalues of a tridiagonal
*>    matrix T which are less than or equal to SHIFT. T has
*>    diagonal entries A(1), ... , A(N), and offdiagonal entries
*>    B(1), ..., B(N-1).
*>    See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
*>    Matrix", Report CS41, Computer Science Dept., Stanford
*>    University, July 21, 1966
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The dimension of the tridiagonal matrix T.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*>          A is DOUBLE PRECISION array, dimension (N)
*>          The diagonal entries of the tridiagonal matrix T.
*> \endverbatim
*>
*> \param[in] B
*> \verbatim
*>          B is DOUBLE PRECISION array, dimension (N-1)
*>          The offdiagonal entries of the tridiagonal matrix T.
*> \endverbatim
*>
*> \param[in] SHIFT
*> \verbatim
*>          SHIFT is DOUBLE PRECISION
*>          The shift, used as described under Purpose.
*> \endverbatim
*>
*> \param[out] NUM
*> \verbatim
*>          NUM is INTEGER
*>          The number of eigenvalues of T less than or equal
*>          to SHIFT.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date November 2011
*
*> \ingroup double_eig
*
*  =====================================================================
      SUBROUTINE dstect( N, A, B, SHIFT, NUM )
*
*  -- LAPACK test routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      INTEGER            N, NUM
      DOUBLE PRECISION   SHIFT
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   A( * ), B( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, THREE
      parameter                ( zero = 0.0d0, one = 1.0d0, three = 3.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 = dlamch( 'Safe minimum' )
      ovfl = dlamch( 'Overflow' )
*
*     Find largest entry
*
      mx = abs( a( 1 ) )
      DO 10 i = 1, n - 1
         mx = max( mx, abs( a( i+1 ) ), abs( b( i ) ) )
   10 CONTINUE
*
*     Handle easy cases, including zero matrix
*
      IF( shift.GE.three*mx ) THEN
         num = n
         RETURN
      END IF
      IF( shift.LT.-three*mx ) THEN
         num = 0
         RETURN
      END IF
*
*     Compute scale factors as in Kahan's report
*     At this point, MX .NE. 0 so we can divide by it
*
      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
*
      num = 0
      sshift = ( shift*m1 )*m2
      u = ( a( 1 )*m1 )*m2 - 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 = 2, n
         tmp = ( b( i-1 )*m1 )*m2
         u = ( ( a( i )*m1 )*m2-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 DSTECT
*
      END