LAPACK  3.4.2 LAPACK: Linear Algebra PACKage
dsvdct.f
Go to the documentation of this file.
1 *> \brief \b DSVDCT
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 DSVDCT( N, S, E, SHIFT, NUM )
12 *
13 * .. Scalar Arguments ..
14 * INTEGER N, NUM
15 * DOUBLE PRECISION SHIFT
16 * ..
17 * .. Array Arguments ..
18 * DOUBLE PRECISION E( * ), S( * )
19 * ..
20 *
21 *
22 *> \par Purpose:
23 * =============
24 *>
25 *> \verbatim
26 *>
27 *> DSVDCT counts the number NUM of eigenvalues of a 2*N by 2*N
28 *> tridiagonal matrix T which are less than or equal to SHIFT. T is
29 *> formed by putting zeros on the diagonal and making the off-diagonals
30 *> equal to S(1), E(1), S(2), E(2), ... , E(N-1), S(N). If SHIFT is
31 *> positive, NUM is equal to N plus the number of singular values of a
32 *> bidiagonal matrix B less than or equal to SHIFT. Here B has diagonal
33 *> entries S(1), ..., S(N) and superdiagonal entries E(1), ... E(N-1).
34 *> If SHIFT is negative, NUM is equal to the number of singular values
35 *> of B greater than or equal to -SHIFT.
36 *>
37 *> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
38 *> Matrix", Report CS41, Computer Science Dept., Stanford University,
39 *> July 21, 1966
40 *> \endverbatim
41 *
42 * Arguments:
43 * ==========
44 *
45 *> \param[in] N
46 *> \verbatim
47 *> N is INTEGER
48 *> The dimension of the bidiagonal matrix B.
49 *> \endverbatim
50 *>
51 *> \param[in] S
52 *> \verbatim
53 *> S is DOUBLE PRECISION array, dimension (N)
54 *> The diagonal entries of the bidiagonal matrix B.
55 *> \endverbatim
56 *>
57 *> \param[in] E
58 *> \verbatim
59 *> E is DOUBLE PRECISION array of dimension (N-1)
60 *> The superdiagonal entries of the bidiagonal matrix B.
61 *> \endverbatim
62 *>
63 *> \param[in] SHIFT
64 *> \verbatim
65 *> SHIFT is DOUBLE PRECISION
66 *> The shift, used as described under Purpose.
67 *> \endverbatim
68 *>
69 *> \param[out] NUM
70 *> \verbatim
71 *> NUM is INTEGER
72 *> The number of eigenvalues of T less than or equal to SHIFT.
73 *> \endverbatim
74 *
75 * Authors:
76 * ========
77 *
78 *> \author Univ. of Tennessee
79 *> \author Univ. of California Berkeley
80 *> \author Univ. of Colorado Denver
81 *> \author NAG Ltd.
82 *
83 *> \date November 2011
84 *
85 *> \ingroup double_eig
86 *
87 * =====================================================================
88  SUBROUTINE dsvdct( N, S, E, SHIFT, NUM )
89 *
90 * -- LAPACK test routine (version 3.4.0) --
91 * -- LAPACK is a software package provided by Univ. of Tennessee, --
92 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
93 * November 2011
94 *
95 * .. Scalar Arguments ..
96  INTEGER n, num
97  DOUBLE PRECISION shift
98 * ..
99 * .. Array Arguments ..
100  DOUBLE PRECISION e( * ), s( * )
101 * ..
102 *
103 * =====================================================================
104 *
105 * .. Parameters ..
106  DOUBLE PRECISION one
107  parameter( one = 1.0d0 )
108  DOUBLE PRECISION zero
109  parameter( zero = 0.0d0 )
110 * ..
111 * .. Local Scalars ..
112  INTEGER i
113  DOUBLE PRECISION m1, m2, mx, ovfl, sov, sshift, ssun, sun, tmp,
114  \$ tom, u, unfl
115 * ..
116 * .. External Functions ..
117  DOUBLE PRECISION dlamch
118  EXTERNAL dlamch
119 * ..
120 * .. Intrinsic Functions ..
121  INTRINSIC abs, max, sqrt
122 * ..
123 * .. Executable Statements ..
124 *
125 * Get machine constants
126 *
127  unfl = 2*dlamch( 'Safe minimum' )
128  ovfl = one / unfl
129 *
130 * Find largest entry
131 *
132  mx = abs( s( 1 ) )
133  DO 10 i = 1, n - 1
134  mx = max( mx, abs( s( i+1 ) ), abs( e( i ) ) )
135  10 continue
136 *
137  IF( mx.EQ.zero ) THEN
138  IF( shift.LT.zero ) THEN
139  num = 0
140  ELSE
141  num = 2*n
142  END IF
143  return
144  END IF
145 *
146 * Compute scale factors as in Kahan's report
147 *
148  sun = sqrt( unfl )
149  ssun = sqrt( sun )
150  sov = sqrt( ovfl )
151  tom = ssun*sov
152  IF( mx.LE.one ) THEN
153  m1 = one / mx
154  m2 = tom
155  ELSE
156  m1 = one
157  m2 = tom / mx
158  END IF
159 *
160 * Begin counting
161 *
162  u = one
163  num = 0
164  sshift = ( shift*m1 )*m2
165  u = -sshift
166  IF( u.LE.sun ) THEN
167  IF( u.LE.zero ) THEN
168  num = num + 1
169  IF( u.GT.-sun )
170  \$ u = -sun
171  ELSE
172  u = sun
173  END IF
174  END IF
175  tmp = ( s( 1 )*m1 )*m2
176  u = -tmp*( tmp / u ) - sshift
177  IF( u.LE.sun ) THEN
178  IF( u.LE.zero ) THEN
179  num = num + 1
180  IF( u.GT.-sun )
181  \$ u = -sun
182  ELSE
183  u = sun
184  END IF
185  END IF
186  DO 20 i = 1, n - 1
187  tmp = ( e( i )*m1 )*m2
188  u = -tmp*( tmp / u ) - sshift
189  IF( u.LE.sun ) THEN
190  IF( u.LE.zero ) THEN
191  num = num + 1
192  IF( u.GT.-sun )
193  \$ u = -sun
194  ELSE
195  u = sun
196  END IF
197  END IF
198  tmp = ( s( i+1 )*m1 )*m2
199  u = -tmp*( tmp / u ) - sshift
200  IF( u.LE.sun ) THEN
201  IF( u.LE.zero ) THEN
202  num = num + 1
203  IF( u.GT.-sun )
204  \$ u = -sun
205  ELSE
206  u = sun
207  END IF
208  END IF
209  20 continue
210  return
211 *
212 * End of DSVDCT
213 *
214  END