LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
ssvdct.f
Go to the documentation of this file.
1 *> \brief \b SSVDCT
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 SSVDCT( N, S, E, SHIFT, NUM )
12 *
13 * .. Scalar Arguments ..
14 * INTEGER N, NUM
15 * REAL SHIFT
16 * ..
17 * .. Array Arguments ..
18 * REAL E( * ), S( * )
19 * ..
20 *
21 *
22 *> \par Purpose:
23 * =============
24 *>
25 *> \verbatim
26 *>
27 *> SSVDCT 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 REAL array, dimension (N)
54 *> The diagonal entries of the bidiagonal matrix B.
55 *> \endverbatim
56 *>
57 *> \param[in] E
58 *> \verbatim
59 *> E is REAL 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 REAL
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 *> \ingroup single_eig
84 *
85 * =====================================================================
86  SUBROUTINE ssvdct( N, S, E, SHIFT, NUM )
87 *
88 * -- LAPACK test routine --
89 * -- LAPACK is a software package provided by Univ. of Tennessee, --
90 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
91 *
92 * .. Scalar Arguments ..
93  INTEGER N, NUM
94  REAL SHIFT
95 * ..
96 * .. Array Arguments ..
97  REAL E( * ), S( * )
98 * ..
99 *
100 * =====================================================================
101 *
102 * .. Parameters ..
103  REAL ONE
104  parameter( one = 1.0e0 )
105  REAL ZERO
106  parameter( zero = 0.0e0 )
107 * ..
108 * .. Local Scalars ..
109  INTEGER I
110  REAL M1, M2, MX, OVFL, SOV, SSHIFT, SSUN, SUN, TMP,
111  $ TOM, U, UNFL
112 * ..
113 * .. External Functions ..
114  REAL SLAMCH
115  EXTERNAL slamch
116 * ..
117 * .. Intrinsic Functions ..
118  INTRINSIC abs, max, sqrt
119 * ..
120 * .. Executable Statements ..
121 *
122 * Get machine constants
123 *
124  unfl = 2*slamch( 'Safe minimum' )
125  ovfl = one / unfl
126 *
127 * Find largest entry
128 *
129  mx = abs( s( 1 ) )
130  DO 10 i = 1, n - 1
131  mx = max( mx, abs( s( i+1 ) ), abs( e( i ) ) )
132  10 CONTINUE
133 *
134  IF( mx.EQ.zero ) THEN
135  IF( shift.LT.zero ) THEN
136  num = 0
137  ELSE
138  num = 2*n
139  END IF
140  RETURN
141  END IF
142 *
143 * Compute scale factors as in Kahan's report
144 *
145  sun = sqrt( unfl )
146  ssun = sqrt( sun )
147  sov = sqrt( ovfl )
148  tom = ssun*sov
149  IF( mx.LE.one ) THEN
150  m1 = one / mx
151  m2 = tom
152  ELSE
153  m1 = one
154  m2 = tom / mx
155  END IF
156 *
157 * Begin counting
158 *
159  u = one
160  num = 0
161  sshift = ( shift*m1 )*m2
162  u = -sshift
163  IF( u.LE.sun ) THEN
164  IF( u.LE.zero ) THEN
165  num = num + 1
166  IF( u.GT.-sun )
167  $ u = -sun
168  ELSE
169  u = sun
170  END IF
171  END IF
172  tmp = ( s( 1 )*m1 )*m2
173  u = -tmp*( tmp / u ) - sshift
174  IF( u.LE.sun ) THEN
175  IF( u.LE.zero ) THEN
176  num = num + 1
177  IF( u.GT.-sun )
178  $ u = -sun
179  ELSE
180  u = sun
181  END IF
182  END IF
183  DO 20 i = 1, n - 1
184  tmp = ( e( i )*m1 )*m2
185  u = -tmp*( tmp / u ) - sshift
186  IF( u.LE.sun ) THEN
187  IF( u.LE.zero ) THEN
188  num = num + 1
189  IF( u.GT.-sun )
190  $ u = -sun
191  ELSE
192  u = sun
193  END IF
194  END IF
195  tmp = ( s( i+1 )*m1 )*m2
196  u = -tmp*( tmp / u ) - sshift
197  IF( u.LE.sun ) THEN
198  IF( u.LE.zero ) THEN
199  num = num + 1
200  IF( u.GT.-sun )
201  $ u = -sun
202  ELSE
203  u = sun
204  END IF
205  END IF
206  20 CONTINUE
207  RETURN
208 *
209 * End of SSVDCT
210 *
211  END
subroutine ssvdct(N, S, E, SHIFT, NUM)
SSVDCT
Definition: ssvdct.f:87