LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*> \ingroup double_eig
84*
85* =====================================================================
86 SUBROUTINE dsvdct( 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 DOUBLE PRECISION SHIFT
95* ..
96* .. Array Arguments ..
97 DOUBLE PRECISION E( * ), S( * )
98* ..
99*
100* =====================================================================
101*
102* .. Parameters ..
103 DOUBLE PRECISION ONE
104 parameter( one = 1.0d0 )
105 DOUBLE PRECISION ZERO
106 parameter( zero = 0.0d0 )
107* ..
108* .. Local Scalars ..
109 INTEGER I
110 DOUBLE PRECISION M1, M2, MX, OVFL, SOV, SSHIFT, SSUN, SUN, TMP,
111 $ TOM, U, UNFL
112* ..
113* .. External Functions ..
114 DOUBLE PRECISION DLAMCH
115 EXTERNAL dlamch
116* ..
117* .. Intrinsic Functions ..
118 INTRINSIC abs, max, sqrt
119* ..
120* .. Executable Statements ..
121*
122* Get machine constants
123*
124 unfl = 2*dlamch( '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 DSVDCT
210*
211 END
subroutine dsvdct(n, s, e, shift, num)
DSVDCT
Definition dsvdct.f:87