LAPACK  3.4.2 LAPACK: Linear Algebra PACKage
dstect.f
Go to the documentation of this file.
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