TESTING/EIG/dstect.f

   1 *> \brief \b DSTECT
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       SUBROUTINE DSTECT( N, A, B, SHIFT, NUM )
  12 *
  13 *       .. Scalar Arguments ..
  14 *       INTEGER            N, NUM
  15 *       DOUBLE PRECISION   SHIFT
  16 *       ..
  17 *       .. Array Arguments ..
  18 *       DOUBLE PRECISION   A( * ), B( * )
  19 *       ..
  20 *
  21 *
  22 *> \par Purpose:
  23 *  =============
  24 *>
  25 *> \verbatim
  26 *>
  27 *>    DSTECT counts the number NUM of eigenvalues of a tridiagonal
  28 *>    matrix T which are less than or equal to SHIFT. T has
  29 *>    diagonal entries A(1), ... , A(N), and offdiagonal entries
  30 *>    B(1), ..., B(N-1).
  31 *>    See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
  32 *>    Matrix", Report CS41, Computer Science Dept., Stanford
  33 *>    University, July 21, 1966
  34 *> \endverbatim
  35 *
  36 *  Arguments:
  37 *  ==========
  38 *
  39 *> \param[in] N
  40 *> \verbatim
  41 *>          N is INTEGER
  42 *>          The dimension of the tridiagonal matrix T.
  43 *> \endverbatim
  44 *>
  45 *> \param[in] A
  46 *> \verbatim
  47 *>          A is DOUBLE PRECISION array, dimension (N)
  48 *>          The diagonal entries of the tridiagonal matrix T.
  49 *> \endverbatim
  50 *>
  51 *> \param[in] B
  52 *> \verbatim
  53 *>          B is DOUBLE PRECISION array, dimension (N-1)
  54 *>          The offdiagonal entries of the tridiagonal matrix T.
  55 *> \endverbatim
  56 *>
  57 *> \param[in] SHIFT
  58 *> \verbatim
  59 *>          SHIFT is DOUBLE PRECISION
  60 *>          The shift, used as described under Purpose.
  61 *> \endverbatim
  62 *>
  63 *> \param[out] NUM
  64 *> \verbatim
  65 *>          NUM is INTEGER
  66 *>          The number of eigenvalues of T less than or equal
  67 *>          to SHIFT.
  68 *> \endverbatim
  69 *
  70 *  Authors:
  71 *  ========
  72 *
  73 *> \author Univ. of Tennessee
  74 *> \author Univ. of California Berkeley
  75 *> \author Univ. of Colorado Denver
  76 *> \author NAG Ltd.
  77 *
  78 *> \date November 2011
  79 *
  80 *> \ingroup double_eig
  81 *
  82 *  =====================================================================
  83       SUBROUTINE DSTECT( N, A, B, SHIFT, NUM )
  84 *
  85 *  -- LAPACK test routine (version 3.4.0) --
  86 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  87 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  88 *     November 2011
  89 *
  90 *     .. Scalar Arguments ..
  91       INTEGER            N, NUM
  92       DOUBLE PRECISION   SHIFT
  93 *     ..
  94 *     .. Array Arguments ..
  95       DOUBLE PRECISION   A( * ), B( * )
  96 *     ..
  97 *
  98 *  =====================================================================
  99 *
 100 *     .. Parameters ..
 101       DOUBLE PRECISION   ZERO, ONE, THREE
 102       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, THREE = 3.0D0 )
 103 *     ..
 104 *     .. Local Scalars ..
 105       INTEGER            I
 106       DOUBLE PRECISION   M1, M2, MX, OVFL, SOV, SSHIFT, SSUN, SUN, TMP,
 107      $                   TOM, U, UNFL
 108 *     ..
 109 *     .. External Functions ..
 110       DOUBLE PRECISION   DLAMCH
 111       EXTERNAL           DLAMCH
 112 *     ..
 113 *     .. Intrinsic Functions ..
 114       INTRINSIC          ABS, MAX, SQRT
 115 *     ..
 116 *     .. Executable Statements ..
 117 *
 118 *     Get machine constants
 119 *
 120       UNFL = DLAMCH( 'Safe minimum' )
 121       OVFL = DLAMCH( 'Overflow' )
 122 *
 123 *     Find largest entry
 124 *
 125       MX = ABS( A( 1 ) )
 126       DO 10 I = 1, N - 1
 127          MX = MAX( MX, ABS( A( I+1 ) ), ABS( B( I ) ) )
 128    10 CONTINUE
 129 *
 130 *     Handle easy cases, including zero matrix
 131 *
 132       IF( SHIFT.GE.THREE*MX ) THEN
 133          NUM = N
 134          RETURN
 135       END IF
 136       IF( SHIFT.LT.-THREE*MX ) THEN
 137          NUM = 0
 138          RETURN
 139       END IF
 140 *
 141 *     Compute scale factors as in Kahan's report
 142 *     At this point, MX .NE. 0 so we can divide by it
 143 *
 144       SUN = SQRT( UNFL )
 145       SSUN = SQRT( SUN )
 146       SOV = SQRT( OVFL )
 147       TOM = SSUN*SOV
 148       IF( MX.LE.ONE ) THEN
 149          M1 = ONE / MX
 150          M2 = TOM
 151       ELSE
 152          M1 = ONE
 153          M2 = TOM / MX
 154       END IF
 155 *
 156 *     Begin counting
 157 *
 158       NUM = 0
 159       SSHIFT = ( SHIFT*M1 )*M2
 160       U = ( A( 1 )*M1 )*M2 - SSHIFT
 161       IF( U.LE.SUN ) THEN
 162          IF( U.LE.ZERO ) THEN
 163             NUM = NUM + 1
 164             IF( U.GT.-SUN )
 165      $         U = -SUN
 166          ELSE
 167             U = SUN
 168          END IF
 169       END IF
 170       DO 20 I = 2, N
 171          TMP = ( B( I-1 )*M1 )*M2
 172          U = ( ( A( I )*M1 )*M2-TMP*( TMP / U ) ) - SSHIFT
 173          IF( U.LE.SUN ) THEN
 174             IF( U.LE.ZERO ) THEN
 175                NUM = NUM + 1
 176                IF( U.GT.-SUN )
 177      $            U = -SUN
 178             ELSE
 179                U = SUN
 180             END IF
 181          END IF
 182    20 CONTINUE
 183       RETURN
 184 *
 185 *     End of DSTECT
 186 *
 187       END